Christian Tanguy
email: christian.tanguy@orange-ft.com

Exact solutions for the two-and all-terminal reliabilities of a simple ladder network

Keywords: network reliability, star-triangle transformation, transfer matrix, zeros of the reliability polynomial, sensitivity, algebraic structures

The exact calculation of network reliability in a probabilistic context has been a long-standing issue of practical importance, but a difficult one, even for planar graphs, with perfect nodes and with edges of identical reliability p. Many approaches (determination of bounds, sums of disjoint products algorithms, Monte Carlo evaluations, studies of the reliability polynomials, etc.) can only provide approximations when the network's size increases.

We consider here a ladder graph of arbitrary size corresponding to real-life network configurations, and give the exact, analytical solutions for the all-and two-terminal reliabilities. These solutions use transfer matrices, in which individual reliabilities of edges and nodes are taken into account. The special case of identical edge and node reliabilities -p and ρ, respectively -is solved. We show that the zeros of the two-terminal reliability polynomial exhibit structures which differ substantially for seemingly similar networks, and we compare the sensitivity of various edges. We discuss how the present work may be further extended to lead to a catalog of exactly solvable networks in terms of reliability, which could be useful as elementary bricks for a new and improved set of bounds or benchmarks in the general case.

Introduction

Network reliability has long been a practical issue, and will remain so for years, since networks have now entered an era of Quality of Service (QoS). IP networks, mobile phone networks, transportation networks, electrical power networks, etc., have become "commodities." Connection availability rates of 99.999% are now an objective for telecommunication network operators, and premium services may only be deployed -and correspondingly billed -if the connection reliability is close enough to unity. Reliability is therefore a crucial parameter in the design and analysis of the various kinds of networks we daily use.

Not surprisingly, the study of network reliability has led to a huge body of literature, starting with the work of Moore and Shannon [START_REF] Moore | Reliable circuits using less reliable relays[END_REF], and including excellent textbooks and surveys [START_REF] Ball | Handbooks in operations research and management science[END_REF][START_REF] Barlow | Mathematical theory of reliability[END_REF][START_REF] Colbourn | The combinatorics of network reliability[END_REF][START_REF] Shier | Network reliability and algebraic structures[END_REF][START_REF] Shooman | Probabilistic reliability: an engineering approach[END_REF]. In what follows, we consider a probabilistic approach, in which the network is represented by an undirected graph G = (V, E), where V is a set of nodes (also called vertices) and E is a set of undirected edges (or links), each of which having a probability p n or p e to operate correctly. Failures of the different constituents are assumed to occur at random, and to be statistically independent events. Among the different measures of reliability, one often considers the k-terminal reliability, namely the probability that a given subset K of k nodes (K ⊂ E) are connected. The most common instances are the all-terminal reliability Rel A (K ≡ E) and the two-terminal reliability Rel 2 (s → t), which deals with a particular connection between a source s and a destination t. Both of them are affine functions of each p n and p e .

The sheer number of possible system states, namely 2 |E|+|V | , clearly precludes the use of an "enumeration of states" strategy for realistic networks, and shows that the final expression may be extremely cumbersome. Consequently, most studies have considered graphs with perfect nodes (p n ≡ 1) and edges of identical reliability p; radio broadcast networks have also been described by networks with perfectly reliable edges but imperfect nodes [START_REF] Aboelfotoh | Computing 2-terminal reliability for radio-broadcast networks[END_REF][START_REF] Graver | You may rely on the reliability polynomial for much more than you might think[END_REF]. It was shown early on -see for instance the discussion in [START_REF] Colbourn | The combinatorics of network reliability[END_REF][START_REF] Welsh | Complexity: Knots, Colourings and Counting[END_REF] -that the calculation of k-terminal reliability is #P-hard in the general case, even with the following simplifying and restricting assumptions that (i) the graph is planar (ii) all nodes are perfectly reliable (iii) all edges have the same reliability p. All reliabilities are then expressed as a polynomial in p, called the reliability polynomial.

The difficulty of the problem has stimulated many approaches: partitioning techniques [START_REF] Dotson | A new analysis technique for probabilistic graphs[END_REF][START_REF] Yoo | A comparison of algorithms for terminal-pair reliability[END_REF], sum of disjoint products [START_REF] Abraham | An improved algorithm for network reliability[END_REF][START_REF] Balan | Preprocessing minpaths for sum of disjoint products[END_REF][START_REF] Heidtmann | Smaller sums of disjoint products by subproduct inversion[END_REF][START_REF] Rai | A survey of efficient reliability computation using disjoint products approach[END_REF][START_REF] Rauzy | A practical comparison of methods to assess sum-of-products[END_REF][START_REF] Soh | Experimental results on preprocessing of path/cut terms in sum of disjoint products technique[END_REF], graph simplifications (series-parallel reductions [START_REF] Moore | Reliable circuits using less reliable relays[END_REF], delta-wye transformations [START_REF] Chari | The delta-wye approximation procedure for two-terminal reliability[END_REF][START_REF] Egeland | On dependence and reliability computations[END_REF][START_REF] Gadani | System effectiveness evaluation using star and delta transformations[END_REF][START_REF] Rosenthal | Transformations for simplifying network reliability calculations[END_REF][START_REF] Wang | Transformations of star-delta and delta-star reliability networks[END_REF], factoring [START_REF] Wood | A factoring algorithm using polygon-to-chain reductions for computing K-terminal network reliability[END_REF]), determination of various lower and upper bounds to the reliability polynomial [START_REF] Ball | Handbooks in operations research and management science[END_REF][START_REF] Beichelt | Bounds on the reliability of binary coherent systems[END_REF][START_REF] Brecht | Improving reliability bounds in computer networks[END_REF][START_REF] Brown | Non-Stanley bounds for network reliability[END_REF][START_REF] Chen | Bounds on the reliability of distributed systems with unreliable nodes & links[END_REF][START_REF] Colbourn | The combinatorics of network reliability[END_REF][START_REF] Scott Provan | Bounds on the reliability of networks[END_REF], Monte-Carlo simulations [START_REF] Fishman | A comparison of four Monte Carlo methods for estimating the probability of st connectedness[END_REF][START_REF] Karger | A randomized fully polynomial time approximation scheme for the all terminal reliability problem[END_REF][START_REF] Nel | Combining Monte Carlo estimates and bounds for network reliability[END_REF], genetic [START_REF] Coit | Reliability optimization of series-parallel systems using a genetic algorithm[END_REF] and ordered binary decision diagram (OBDD) algorithms [START_REF] Kuo | Determining terminal pair reliability based on edge expansion diagrams using OBDD[END_REF][START_REF] Rauzy | A new methodology to handle Boolean models with loops[END_REF][START_REF] Yeh | OBDD-based evaluation of k-terminal network reliability[END_REF][START_REF] Yeh | Analyzing network reliability with imperfect nodes using OBDD[END_REF]. The reliability polynomial has also been extensively studied [START_REF] Chari | Reliability polynomials: a survey[END_REF][START_REF] Colbourn | Some open problems on reliability polynomials[END_REF][START_REF] Oxley | Chromatic, flow and reliability polynomials: the complexity of their coefficients[END_REF], with the aim of deriving some useful and hopefully general information from the structure of its coefficients [START_REF] Chari | Reliability polynomials: a survey[END_REF][START_REF] Colbourn | Some open problems on reliability polynomials[END_REF] or the location of its zeros in the complex plane [START_REF] Brown | Roots of the reliability polynomial[END_REF].

A recent breakthrough has been obtained by statistical physicists, who observed that the Tutte polynomial T (G, x, y) of a graph G is actually equivalent to the Potts model partition function of the q-state Potts model [START_REF] Shrock | Exact Potts model partition functions on ladder graphs[END_REF][START_REF] Welsh | The Potts model and the Tutte polynomial[END_REF], which they were able to calculate for various recursive families of graphs [START_REF] Chang | Reliability polynomials and their asymptotic limits for families of graphs[END_REF]. The all-terminal reliability Rel A (G, p) of such graphs is then deduced from T (G, 1, 1 1-p ). Royle and Sokal [START_REF] Royle | The Brown-Colbourn conjecture on zeros of reliability polynomials is false[END_REF] proved that the Brown-Colbourn conjecture [START_REF] Brown | Roots of the reliability polynomial[END_REF] on the location of the zeros of the all-terminal reliability polynomial, while valid for series-parallel reducible graphs, does not hold for some families of graphs. While these results are extremely valuable to understand a few properties of graphs and all-terminal reliability polynomials, they still assume that nodes are perfect and that edges have the same reliability.

In recent years, the tremendous growth of Internet traffic has called for a better evaluation of the reliability of connections in optical networks. This, of course, strongly depends on the connection under consideration. Actual failure rates and maintenance data show that a proper evaluation of twoterminal reliabilities must put node and edge equipments on an equal footing, i.e., both edge (fiber links, optical amplifiers) and node (optical cross-connects, routers) failures must be taken into account. The possibility of node failure has been considered in early papers [START_REF] Aboelfotoh | Computing 2-terminal reliability for radio-broadcast networks[END_REF][START_REF] Evans | Optimally reliable graphs for both edge and vertex failures[END_REF][START_REF] Hänsler | Exact Calculation of Computer Network Reliability[END_REF], to quote but a few. Adaptation of algorithms to include imperfect nodes has been -sometimes controversially -addressed [START_REF] Ke | Reliability evaluation for distributed computing networks with imperfect nodes[END_REF][START_REF] Netes | Consideration of node failures in network-reliability calculation[END_REF][START_REF] Theologou | Factoring & reductions for networks with imperfect vertices[END_REF][START_REF] Torrieri | Calculation of node-pair reliability in large networks with unreliable nodes[END_REF][START_REF] Yeh | Analyzing network reliability with imperfect nodes using OBDD[END_REF]. The two-variable approach for bounds to the reliability polynomial, by Bulka and Dugan [START_REF] Bulka | Network s-t reliability bounds using a 2-dimensional reliability polynomial[END_REF] and Chen and He [START_REF] Chen | Bounds on the reliability of distributed systems with unreliable nodes & links[END_REF], is also worth mentioning. In order to be realistic, different edge reliabilities should be used too: for instance, the failure rate of optical fiber links is often assumed to increase with their length.

In this work, we give the exact solution to the two-terminal reliability of a simple ladder network, displayed in Fig. 1, where successive nodes are labelled S i or T j , and where the larger black dots mark the source s and terminal t we consider in our two-terminal reliability calculations (we choose these special nodes to lie at the graph extremities, since we can always reduce to this case by a series-parallel simplification). This network is a simplified description of a standard (nominal + backup paths) architecture, with additional connections between transit nodes for enabling the so-called "local protection" policy, which bypasses faulty intermediate nodes and/or edges. Such an architecture of "absolutely reliable nodes and unreliable edges," with up to 25 edges, was chosen as Example 5 in [START_REF] Heidtmann | Smaller sums of disjoint products by subproduct inversion[END_REF] for a comparison of different "sum of disjoint products" minimizing algorithms, or by Rauzy [START_REF] Rauzy | A new methodology to handle Boolean models with loops[END_REF] as well as Kuo and collaborators [START_REF] Kuo | Determining terminal pair reliability based on edge expansion diagrams using OBDD[END_REF][START_REF] Yeh | OBDD-based evaluation of k-terminal network reliability[END_REF][START_REF] Yeh | Analyzing network reliability with imperfect nodes using OBDD[END_REF]. By letting the individual node and edge reliabilities take arbitrary values, we actually do not add to the complexity of the problem but make the internal structure of the problem more discernible (a similar approach has been fruitful in the context of graph coloring [START_REF] Biggs | Recursive families of graphs[END_REF][START_REF] Biggs | Matrix method for chromatic polynomials[END_REF]). Indeed, it is then easier to exploit to the full the recursive nature of the ladder graph while using the delta-wye transformation for graphs with unreliable nodes [START_REF] Gadani | System effectiveness evaluation using star and delta transformations[END_REF].

We show that the two-terminal reliability has a beautiful algebraic structure, as its exact expression is given by a product of 3 × 3 transfer matrices (see eqs. ( 13) and ( 16) below for the configuration of Fig. 1(a), or eqs. ( 17) and [START_REF] Chari | The delta-wye approximation procedure for two-terminal reliability[END_REF] for Fig. 1(b)). Consequently, it can also be determined for an arbitrary size (length) of the network. Our aim is (i) to give a detailed derivation of the final results, so that researchers or engineers involved in reliability studies can readily use an easyto-implement formula, since worksheet applications are well up to the task (ii) evaluate the implications of these findings, which suggest a new and potentially fruitful approach to reliability in the case of more general graphs or in combinatorial problems such as the enumeration of self-avoiding walks in lattices of restricted width (iii) emphasize anew the importance of algebraic structures of the underlying graphs in the determination of network reliability [START_REF] Biggs | Algebraic graph theory[END_REF][START_REF] Shier | Network reliability and algebraic structures[END_REF].

Our paper is organized as follows: In Section 2, we briefly recall the formulae for the delta-wye transformation for unreliable nodes. In Section 3, we define the notations for the different edge and node reliabilities, used in the detailed derivation of the main results. We then consider in Section 4 the case where all edges and nodes have identical reliabilities p and ρ, respectively, the size of the network appearing simply as an integer n. We give there the analytic solution of the two-terminal reliability Rel 2 (p, ρ; n) for different configurations, and their asymptotic form when n → ∞. In Section 5, we expand the preceding expressions for p and ρ close to unity, and determine when it is possible to approximate the two-terminal reliability in terms of multiple failures. We also provide in Section 6 the associated generating functions, a very useful tool in combinatorics perfectly suited to reliability studies, since they encode all the necessary information in a beautifully simple form. Prompted by the nearly universal character of the Brown-Colbourn conjecture [START_REF] Brown | Roots of the reliability polynomial[END_REF], we then address the location of zeros of the two-terminal reliability polynomials in Section 7, and show that their structures may well be quite distinct even for seemingly similar networks. We then give in Section 8 a short glimpse of the sensitivity issue [START_REF] Henley | Probabilistic Risk Assessment: Reliability Engineering, Design and Analysis[END_REF][START_REF] Rubino | Sensitivity computation in network reliability analysis[END_REF], namely the influence of a given component to the overall reliability, by comparing the contributions of particular edges (the rungs of the ladder). For the sake of completeness, we derive in Section 9 the all-terminal reliability for the ladder networks under consideration, for arbitrary values of edge reliabilities. Finally, we conclude by indicating several directions in which the present results may be further extended, so that, for instance, a catalog of exactly solvable networks -in terms of reliability -may be given rapidly [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of the Brecht-Colbourn ladder and the generalized fan[END_REF][START_REF] Tanguy | Transfer matrix method for two-and all-terminal reliabilities in ladder[END_REF], which could be useful as elementary bricks for a new and improved set of bounds or benchmarks for alternative methods in the general case.

Triangle-star transformation for unreliable nodes

The triangle-star -also called delta-star, delta-wye, and ∆-Y -transformation has been used many times to simplify calculations of network reliability [START_REF] Chari | The delta-wye approximation procedure for two-terminal reliability[END_REF][START_REF] Colbourn | The combinatorics of network reliability[END_REF][START_REF] Egeland | On dependence and reliability computations[END_REF][START_REF] Gadani | System effectiveness evaluation using star and delta transformations[END_REF][START_REF] Wood | A factoring algorithm using polygon-to-chain reductions for computing K-terminal network reliability[END_REF][START_REF] Rosenthal | Transformations for simplifying network reliability calculations[END_REF][START_REF] Wang | Transformations of star-delta and delta-star reliability networks[END_REF], even though it has mainly been applied in a perfect nodes context, to provide upper and lower bounds to the exact reliability. Here, we fully exploit this transformation in the case of imperfect nodes in order to obtain exact results. Since it plays a crucial part in the derivation, we recall the formulae first derived by Gadani [START_REF] Gadani | System effectiveness evaluation using star and delta transformations[END_REF].

Consider three particular nodes A, B, and C of a network represented in the left part of Fig. 2 (the reliability of the nodes are given by the same uppercase variables, in order to avoid an unnecessary multiplication of notations). The reliability of the edge connecting A and B is given by (lowercase) c, with similar notation for the remaining edges of the triangle. The aim of the triangle-star transformation is to replace the triangle by a star, which is possible by the addition of a new unreliable vertex O and three new edges connecting O to A, B and C, with reliabilities p A , p B , and p C , respectively. Both networks are equivalent provided that the following compatibility relations hold [START_REF] Gadani | System effectiveness evaluation using star and delta transformations[END_REF] Note that the first three equalities correspond to the probability that the two nodes under consideration are connected, while the last gives the probability that all three nodes are connected. A word of caution -given by Gadani -is worth mentioning in the case of successive triangle-star transformations: the triangles should have no common edge or node.

p A O p C = b + a c B -a b c B, (1)
p A O p B = c + a b C -a b c C, (2) p B O p C = a + b c A -a b c A, (3) p A O p B p C = a b + b c + a c -2 a b c. (4)

Derivation of the main results

We first name the different edge and node reliabilities of the ladder diagrams represented in Fig. 1(a), and detail how we can use the triangle-star transformation of the preceding section. Assuming that S 0 is always the source node, and using node reliabilities S i or T j , we note a n the reliability of the edge (S n-1 , S n ), b n for (S n , T n ), and finally c n for (T n-1 , T n ), as shown in Fig. 3(a).

Calculation of Rel

2 (S 0 → T n )
The two-terminal reliability can be written

Rel 2 (S 0 → T n ) = T n F n (b 0 , a 1 , b 1 , c 1 , • • • , a n , b n , c n ). ( 5
)
because the end-node reliabilities obviously factor themselves out of the expression. The transformation from Fig. 3(a) to Fig. 3(b) is a mere series simplification, where a single edge of reliability a n b n S n is now connecting S n-1 to T n . We can then apply the triangle-star transformation to the triangle constituted by {S n-1 , T n-1 , T n }, and define new edges p 1 , p 2 , p 3 , and a node O as represented in Fig. 3(c). From eqs. ( 1)-( 4), we have

p 1 O p 2 = a n b n S n + b n-1 c n T n-1 -a n b n-1 b n c n S n T n-1 , (6)
p 1 O p 3 = c n + b n-1 a n b n S n-1 S n -a n b n-1 b n c n S n-1 S n , (7)
p 2 O p 3 = b n-1 + a n b n c n S n T n -a n b n-1 b n c n S n T n , (8) p 1 O p 2 p 3 = a n b n c n S n + a n b n-1 b n S n + b n-1 c n -2 a n b n-1 b n c n S n . (9)
It is then obvious that the two-terminal reliability of the new network displayed in Fig. 3(c) must be formally similar to the one we started from. The product T n p 1 O (a dangling link) is easily factored, and we are left with the calculation of Rel 2 (S 0 → O). This is nothing but the calculation of Rel 2 (S 0 → T n-1 ), where O has replaced T n-1 and where adjustments are due: b n-1 and c n-1 must be replaced by p 2 and c n-1 p 3 T n-1 , respectively. In short, we have

F n = p 1 O F n-1 (b n-1 → p 2 , c n-1 → c n-1 p 3 T n-1 ). ( 10
)
This relation between F n and F n-1 , in which the recursive nature of the ladder graph is now fully apparent, can be further simplified. Indeed, assuming for F n the simple form

F n = β n b n S n + γ n c n + δ n b n c n S n (11)
in eqs. ( 5) and [START_REF] Biggs | Recursive families of graphs[END_REF], and replacing all occurrences of p 1 , p 2 , p 3 , and O, we easily find

       β n γ n δ n        = M n •        β n-1 γ n-1 δ n-1        , (12)
with the transfer matrix M n given by

M n =        a n S n-1 a n b n-1 c n-1 S n-1 T n-1 a n b n-1 c n-1 S n-1 T n-1 b n-1 S n-1 T n-1 c n-1 T n-1 b n-1 c n-1 S n-1 T n-1 -a n b n-1 S n-1 T n-1 -a n b n-1 c n-1 S n-1 T n-1 a n (1 -2 b n-1 ) c n-1 S n-1 T n-1        .( 13
)
We now need to calculate F 1 . This is straightforward, since we have a simple series-parallel graph:

F 1 = S 0 (a 1 b 1 S 1 + b 0 c 1 T 0 -b 0 a 1 b 1 c 1 S 1 T 0 ), (14)
from which we deduce β 1 = a 1 S 0 , γ 1 = b 0 S 0 T 0 , and δ 1 = -b 0 a 1 S 0 T 0 , or even more simply

       β 1 γ 1 δ 1        = M 1 •        1 0 0        . (15)
Note that eq. ( 15) is correct, even if c 0 is undefined; c 0 may be set to 0 for the sake of simplicity. We finally obtain the beautifully simple and compact result for the two-terminal reliability between nodes S 0 and T n :

Rel 2 (S 0 → T n ; n ≥ 1) = T n • (b n S n , c n , b n c n S n ) • M n M n-1 • • • M 1 •        1 0 0        , (16)
with M n given by eq. ( 13). We can also write

T n • (b n S n , c n , b n c n S n ) = (0, 1, 0) • M n+1 . Note that Rel 2 (S 0 → T 0 ) = b 0 S 0 T 0 may be obtained by replacing M n • • • M 1 ,
which is undefined for n = 0, by the 3 × 3 identity matrix -or no matrix at all. It is also worth remarking that each edge or node reliability -which may be arbitrary -appears once and once only in the different matrices or vectors, respecting the general property that Rel 2 (S 0 → T n ) is an affine function of each component reliability. We have independently checked the above expressions with a sum of disjoint products method for the first values of n. Numerically speaking, the formula given in eq. ( 16) is very easy to implement, even in a worksheet application where multiplication of 3 × 3 matrices is routinely performed; it also easily applies for ladders of arbitrary size, which is also quite important. This is the final answer to the two-terminal reliability of Example 5 of [START_REF] Heidtmann | Smaller sums of disjoint products by subproduct inversion[END_REF], which only gave the number of disjoint products for up to 25 unreliable edges, and of benchmark networks #29-#30 of [START_REF] Kuo | Determining terminal pair reliability based on edge expansion diagrams using OBDD[END_REF][START_REF] Rauzy | A new methodology to handle Boolean models with loops[END_REF][START_REF] Yeh | OBDD-based evaluation of k-terminal network reliability[END_REF][START_REF] Yeh | Analyzing network reliability with imperfect nodes using OBDD[END_REF], where Rel 2 (S 0 → T n ) was calculated for n = 19 and n = 99 using OBDD algorithms.

Calculation of Rel

2 (S 0 → S n )
Let us now turn to the calculation of the two-terminal reliability between two nodes located on the same side of the ladder, as S 0 and S n in Fig. 1(b). The recipe is the same as above: the edges b n and c n are merged with the node T n to give a new equivalent edge between T n-1 and S n , of reliability b n c n T n .

Here again, we can use the triangle-star transformation to derive new recursion relations, and transfer matrices appear. The calculations being straightforward and similar to those in the preceding subsection, we give only the final result

Rel 2 (S 0 → S n ; n ≥ 1) = S n • (a n , b n T n , a n b n T n ) • M n M n-1 • • • M 1 •        1 0 0        , (17)
with

M n =        a n-1 S n-1 b n-1 S n-1 T n-1 a n-1 b n-1 S n-1 T n-1 a n-1 b n-1 c n S n-1 T n-1 c n T n-1 a n-1 b n-1 c n S n-1 T n-1 -a n-1 b n-1 c n S n-1 T n-1 -b n-1 c n S n-1 T n-1 a n-1 (1 -2 b n-1 ) c n S n-1 T n-1        (18)
and the additional convention 17) also applies in the case n = 0, which gives Rel 2 (S 0 → S 0 ) = S 0 , as it should.

a 0 = 1. Since S n • (a n , b n T n , a n b n T n ) = (1, 0, 0)• M n+1 , eq. (
4 Identical reliabilities p and ρ While eqs. ( 16) and ( 17) give the two-terminal reliability in the general case, we consider here the special case where edges and nodes have reliabilities p and ρ, respectively. Our exact results may prove useful because identical edge reliabilities are usually taken for granted, and combinatorial aspects play an important part in the corresponding calculations [START_REF] Colbourn | The combinatorics of network reliability[END_REF]. We show in this section that the two-terminal reliability of the ladder graph has a simple analytic expression for arbitrary n, which we derive from the above expansions.

Calculation of Rel

(S 0 →Tn) 2

(p, ρ)

Taking a i = b i = c i ≡ p, and S i = T i ≡ ρ, we get Rel (S 0 →Tn) 2 (p, ρ) = (p ρ 2 , p ρ, p 2 ρ 2 ) •        p ρ p 3 ρ 2 p 3 ρ 2 p ρ 2 p ρ p 2 ρ 2 -p 2 ρ 2 -p 3 ρ 2 p 2 (1 -2 p) ρ 2        n •               .( 19
)
We could of course replace (p ρ 2 , p ρ, p 2 ρ 2 ) by (0, 1, 0), provided that the matrix exponent is n + 1 instead of n. The next step is to find the eigenvalues and eigenvectors of the matrix appearing in eq. ( 19), so that the final result may be further simplified. Obviously, a p n ρ n factor can be factored out of the matrix contribution. We are left with the determination of the eigenvalues and possibly eigenvectors of

M T =        1 p 2 ρ p 2 ρ ρ 1 p ρ -p ρ -p 2 ρ p (1 -2 p) ρ        . ( 20
)
From the first two rows of M T , 1 -p ρ is obviously an eigenvalue of the matrix. By using standard linear algebra, the characteristic polynomial of M T can be factored, the three eigenvalues x 0 and x ± of M T being simply

x 0 = 1 -p ρ, (21)
x ± = 1 + 2 p (1 -p) ρ ± √ 1 + 4 p 2 ρ -8 p 3 ρ 2 + 4 p 4 ρ 2 2 . ( 22
)
While the determination of the eigenvectors may be used to finish the calculation, this is by no means necessary. Indeed, M T being diagonalizable -the three roots are distinct -Rel

(S 0 →Tn) 2 (p, ρ) must be of the form Rel (S 0 →Tn) 2 (p, ρ) = p n ρ n α 0 x n 0 + α + x n + + α -x n -. (23)
From the first values of Rel

(S 0 →Tn) 2 (p, ρ) for n = 0, 1, 2 Rel (S 0 →T 0 ) 2 (p, ρ) = p ρ 2 , (24)
Rel

(S 0 →T 1 ) 2 (p, ρ) = p 2 ρ 3 2 -p 2 ρ , (25)
Rel

(S 0 →T 2 ) 2 (p, ρ) = p 3 ρ 4 3 -2 p 2 ρ + p 2 ρ 2 (1 -p) (1 -2 p) , (26)
which have been independently checked by a sum of disjoint products procedure, and the inversion of a 3 × 3 Vandermonde matrix based on x 0 and x ± , we obtain α 0 and α ± . Finally Rel

(S 0 →Tn) 2 (p, ρ) = p n ρ n+1 2 -(1 -p ρ) n+1 + (1 + p ρ) x n+1 + -x n+1 - x + -x - -p ρ (1 -2 p + p ρ) x n + -x n - x + -x - . (27)
Setting p = 0.9 and ρ = 0.9 or 1 for n = 19 or 99 confirms the numerical results found for the 2 × 20 and 2 × 100 ladder networks [START_REF] Kuo | Determining terminal pair reliability based on edge expansion diagrams using OBDD[END_REF][START_REF] Rauzy | A new methodology to handle Boolean models with loops[END_REF][START_REF] Yeh | OBDD-based evaluation of k-terminal network reliability[END_REF][START_REF] Yeh | Analyzing network reliability with imperfect nodes using OBDD[END_REF].

We discuss the implications of eq. ( 27) in section 4.4 below. We will see in Section 6 that the simplicity of the result stems from our choosing a unique reliability p for all edges. Before that, we turn to other two-terminal configurations.

Calculation of Rel

(S 0 →Sn) 2 (p, ρ)
In the S 0 → S n configuration, eqs. ( 17) and ( 18) lead to Rel

(S 0 →Sn) 2 (p, ρ) = ρ (p, p ρ, p 2 ρ) •        p ρ p ρ 2 p 2 ρ 2 p 3 ρ 2 p ρ p 3 ρ 2 -p 3 ρ 2 -p 2 ρ 2 p 2 (1 -2 p) ρ 2        n •        1/p 0 0        .( 28
)
Note that while a 0 ≡ 1 is the only a i which should not be set to p, we can compensate this missing p term in M 1 by introducing a 1/p term in the rightmost vector in eq. ( 17). The procedure is the same as above: a (p ρ) n prefactor can safely be extracted from the matrix on the right-hand side of eq. ( 28), so that we have to find the eigenvalues of another matrix, namely

M S M S =        1 ρ p ρ p 2 ρ 1 p 2 ρ -p 2 ρ -p ρ p (1 -2 p) ρ        . (29)
Unsurprisingly, the eigenvalues of M S are again given by x 0 and x ± . Finding the prefactors associated with x n 0 and x n ± from the first three values of Rel (S 0 →Sn) 2 (p, ρ) is straightforward, and the final expression is Rel

(S 0 →Sn) 2 (p, ρ) = p n ρ n+1 2 +(1 -p ρ) n+1 + (1 + p ρ) x n+1 + -x n+1 - x + -x - -p ρ (1 -2 p + p ρ) x n + -x n - x + -x - . (30)
It is worth noting that the expressions given in eqs. ( 27) and ( 30) only differ by a '±' sign.

Calculation of Rel

(S 0 →Un) 2 (p, ρ)
Let us introduce another configuration -the symmetrical ladder network of Fig. 4 -which looks quite similar to the two previous ones because (i) it has been considered in the literature as a way to describe a "gamma multistage interconnection network" [START_REF] Rai | A survey of efficient reliability computation using disjoint products approach[END_REF], an architecture demonstrating the potential usefulness of recursive algorithms by Hänsler [START_REF] Hänsler | Comments on "A fast recursive algorithm to calculate the reliability of a communication network[END_REF] or Fratta and Montanari [START_REF] Fratta | A recursive method based on case analysis for computing network reliability[END_REF] (ii) it will shed some light on what can be expected from the study of the two-terminal reliability polynomial, as regards the location of its zeros. We do not give the exact expression of Rel (S 0 →Un) 2

for arbitrary edge and node reliabilities, since it can be derived from eqs. ( 13) and ( 16) by merely taking T 0 = S n = b 0 = b n = 1. We limit ourselves to the (p, ρ) configuration and find, after a treatment similar to that of the previous sections,

Rel (S 0 →Un) 2 (p, ρ) = p n ρ n+1 x + -x - (2 -p) (x n + -x n -) +p (1 -2 ρ + p ρ) (x n-1 + -x n-1 -) . (31)
Note that x 0 is absent here: only two eigenvalues appear, even though the network architecture differs from the previous ones by two edges only.

Variation of the eigenvalues with p and ρ

From eqs. ( 27), [START_REF] Fratta | A recursive method based on case analysis for computing network reliability[END_REF], and ( 31), it appears that the two-terminal reliabilities for the ladders are given by sums of terms such as λ n 0 , λ n + , and λ n -, with

λ 0 = p ρ (1 -p ρ), (32)
λ ± = 1 2 p ρ 1 + 2 p (1 -p) ρ ± 1 + 4 p 2 ρ -8 p 3 ρ 2 + 4 p 4 ρ 2 , (33)
where n appears only in the exponents. A natural first step is to assess these quantities as p and ρ vary, in order to assess the relative contributions of the different terms. For the sake of simplicity, we have displayed on Fig. 5 the variation of λ 0 , λ + , and λ -with p, when ρ = 1. All three quantities are positive and vanish when p goes to zero, as expected. However, only λ + tends to 1 when p → 1 (remember that for p → 1, Rel 2 → 1, implying that one of the eigenvalues tends to 1). In this particular case, we have λ + > λ 0 > λ -in the entire range 0 < p < 1.

An important issue is the behavior of Rel (S 0 →Tn,Sn,Un) 2

(p, ρ) as the size of the network -or n -increases. We see on Fig. 5 that for p = 0.9, λ + is about ten times larger than λ 0 and λ -. Consequently, when n is sufficiently largethis may be valid for n ∼ 10 already -we have to a very good approximation Rel

(S 0 →Tn,Sn,Un) 2 (p, ρ) ≈ α + λ n + , (34)
which is observed in the numerical evaluations of Kuo et al. and Rauzy [START_REF] Kuo | Determining terminal pair reliability based on edge expansion diagrams using OBDD[END_REF][START_REF] Rauzy | A new methodology to handle Boolean models with loops[END_REF][START_REF] Yeh | OBDD-based evaluation of k-terminal network reliability[END_REF][START_REF] Yeh | Analyzing network reliability with imperfect nodes using OBDD[END_REF]. Basically, adding one "cell" to the ladder network will merely result in another λ + factor in the two-terminal reliability, so that λ + represents a scaling of the reliability; it is as if we had n equipments of reliability λ + in series. When ρ = 1, this general behavior is the same, even if λ -is negative for 1 2-ρ ≤ p ≤ 1, λ + still prevails, and the power-law behavior of eq. ( 34) holds.

Failures decomposition

In some instances, when p and ρ are assumed to be very close to 1, calculations often involve the unreliability of the connection U 2 = 1 -Rel 2 . Let q = 1p and η = 1ρ. The network unreliability is then expanded in q and η, even though one should be careful about such a practice, when q and η differ by orders of magnitude. For instance, keeping only three terms in U 2 , we find from eq. ( 27) that the following formula is valid for n greater than 3 (but not too large, of course)

U 2 (S 0 → T n ) = 2 η + (3 n -4) η 2 -(12 n -16) η 3 + 4 n η q -(20 n -6) η 2 q + (n + 2) q 2 -(4 n + 14) η q 2 +2 (n -1) q 3 + • • • (35)
We must keep in mind, however, that q and η must be sufficiently small for such an expansion to make sense. Limiting ourselves to η = 0 (the perfect nodes configuration), we get

U 2 (S 0 → T n ) = (n + 2) q 2 + 2 (n -1) q 3 - n 2 + 7 n + 10 2 q 4 + • • • (36)
with the same result for U 2 (S 0 → S n ), whereas

U 2 (S 0 → U n ) = n q 2 + 2 (n -1) q 3 - n (n + 3) 2 q 4 + • • • (37)
The "small parameter" in the expansion is therefore of the order of n q. A rule of thumb is to ensure this quantity is small indeed when considering expansions such as eq. ( 35).

Generating functions

There is yet another way to express all the abovementioned results is a simple and compact form, namely that of the so-called "ordinary generating functions", which are a fundamental tool in combinatorics [START_REF] Stanley | Enumerative combinatorics[END_REF] and extensively used in [START_REF] Chang | Reliability polynomials and their asymptotic limits for families of graphs[END_REF]. The matrix transfer formalism and results obtained in the preceding sections clearly indicate the existence of recursion relations between successive reliability polynomials (see [START_REF] Biggs | Recursive families of graphs[END_REF] for an interesting counterpoint for the chromatic polynomials of graph ladders), which we can write in an even more compact form than eqs. ( 27) and [START_REF] Fratta | A recursive method based on case analysis for computing network reliability[END_REF]. If we can write -assuming all the eigenvalues are distinct -

Rel (n) 2 (p, ρ) = 3 i∈{0,+,-} α i λ n i , (38)
the ordinary generating function G(x) defined by

G(x) = ∞ n=0 Rel (n) 2 (p, ρ) x n (39)
is easily calculated, since we have to sum geometric series :

G(x) = 3 i∈{0,+,-} α i 1 -λ i x = N (x) D(x) . ( 40
)
Of course, the summation can be performed too when α i is a polynomial function of n. G(x) is a rational fraction of x (see [START_REF] Stanley | Enumerative combinatorics[END_REF] for a thorough discussion of the topics). Its denominator D(x) = i (1λ i x) may be deduced very simply from the characteristic polynomial P carac (x) of the unfactored transfer matrix by D(x) = (-x) 3 P carac (1/x). The determination of G(x), which encodes all the two-terminal reliabilities for the recursive family of graphs indexed by n, is actually quite simple: we only need to multiply the known leading terms of the expansion of G(x) by D(x) to see N (x) emerge.

For instance, in the case of Rel (S 0 →Tn) 2

(p, ρ), we easily find D T (x) from eq. ( 19)

D T (x) = 1 -p ρ x (2 + p ρ (1 -2 p)) + p 2 ρ 2 x 2 (1 + p ρ (1 -2 p) (2 -p ρ)) -p 4 ρ 4 x 3 (1 -p ρ) (1 -2 p + p ρ) , (41)
and deduce from eqs. ( 24)-( 26) the simple result (mathematical softwares allowing symbolic calculus make such computations an easy task)

N T (x) = p ρ 2 1 -(1 -p) p 2 ρ 2 x . (42)
The preceding treatment works well too if all the a i , b i , and c i are replaced by a, b, and c, respectively, instead of by a unique value p. The final result reads, after simplification

G T (a, b, c, ρ; x) = b ρ 2 N T (a, b, c, ρ; x) D T (a, b, c, ρ; x) , (43)
N T (a, b, c, ρ; x) = 1 -a (1 -b) c ρ 2 x, (44)
D T (a, b, c, ρ; x) = 1 -[a + c + a c ρ -2 a b c ρ] ρ x +a c [1 + (1 -2 b) (a + c) ρ -(1 -a -c) b 2 ρ 2 ] ρ 2 x 2 -(1 -b ρ) (1 -2 b + b ρ) a 2 c 2 ρ 4 x 3 . (45)
When a, b and c are not identical, no simple factorization of D T occurs, so that we must use the roots of a polynomial of degree three in x, the literal expressions of which are not concise. However, the expression of the generating function makes it possible to calculate the two-terminal reliability by a mere partial fraction decomposition, which can be handled quite routinely.

In the case of Rel (S 0 →Sn) 2

, G S (x) = N S (x)/D S (x), with D S (x) = D T (x) and

N S (x) = ρ -p ρ 2 1 + p ρ (1 -2 p) -p 2 ρ 2 (1 -p) x +p 3 ρ 4 1 -2 p + p 2 ρ (2 -ρ) x 2 . ( 46
)
While this expression seems more complicated than eq. ( 42), we can check that

G S (x) -G T (x) = p ρ 2 (1 -p ρ) 2 x 1 -p ρ (1 -p ρ) x , (47)
which could be directly deduced from eqs. ( 27) and [START_REF] Fratta | A recursive method based on case analysis for computing network reliability[END_REF].

The generating function G U for S 0 → U n is given by

G U (x) = p ρ 2 x (2 -p) + p 2 ρ (1 -2 ρ + p ρ) x 1 -p ρ (1 + 2 p (1 -p) ρ) x + p 3 ρ 3 (1 -2 p + p ρ) x 2 . (48)
Here, the denominator is only a polynomial in x of degree two, leading to the decomposition in eq. ( 31). We have included the case n = 1, for which the reliability is that of two links in parallel between S 0 and U 1 , namely ρ 2 (2 pp 2 ); this has the advantage of simplifying the rational G U (x), and making transparent the derivation of eq. ( 31).

Zeros of the two-terminal reliability polynomials

One way to understand the structure of the different reliability polynomials is to study the locations of their zeros in the complex plane. Such a study has been fruitfully performed in the case of the chromatic polynomial [START_REF] Biggs | Recursive families of graphs[END_REF][START_REF] Biggs | Matrix method for chromatic polynomials[END_REF][START_REF] Salas | Transfer matrices and partition-function zeros for antiferromagnetic Potts models. I. General theory and square-lattice chromatic polynomial[END_REF], most notably in the context of the four-color theorem. In the reliability context, some effort has been done to discover general properties for the allterminal reliability Rel A (p) [START_REF] Chari | Reliability polynomials: a survey[END_REF][START_REF] Colbourn | Some open problems on reliability polynomials[END_REF][START_REF] Oxley | Chromatic, flow and reliability polynomials: the complexity of their coefficients[END_REF], its main byproduct being the Brown-Colbourn conjecture [START_REF] Brown | Roots of the reliability polynomial[END_REF], according to which all the zeros are to be found in the region |1 -p| < 1. Although valid for series-parallel graphs, this remarkable conjecture does not strictly hold in the general case (but not by far) [START_REF] Royle | The Brown-Colbourn conjecture on zeros of reliability polynomials is false[END_REF]. As mentioned in the introduction, the all-reliability polynomial is linked to the Tutte polynomial, an invariant of the graph. It has also been studied extensively by Chang and Shrock for various recursive families of graphs [START_REF] Chang | Reliability polynomials and their asymptotic limits for families of graphs[END_REF], who give the limiting curves where all zeros of the polynomials converge.

In this section, we want to locate the zeros of the two-terminal reliability polynomial, in order to see whether some insight may also be found in this case. Truly, this polynomial depends on the couple (source, terminal), but some structures could still be expected. We show indeed that the zeros tend to aggregate along sets of curves, which can substantially differ even for two seemingly similar ladder graphs.

Calculation of the limiting curves: General results

As n grows, the number of zeros of the reliability polynomial in the complex plane increases. Because of the matrix transfer property, we have recursion relations between polynomials corresponding to successive values of n, and the solution looks like eq. ( 23). The general treatment of the problem has been done by Beraha, Kahane, and Weiss [START_REF] Beraha | Limits of zeros of recursively defined families of polynomials[END_REF], but may be understood in the following, simplifying way: if the reliability polynomial has the form i α i λ i (p) n (where λ i are the eigenvalues of the transfer matrix), then at large n, only the two eigenvalues of greater modules, say λ 1 and λ 2 , will prevail, so that the reliability polynomial will vanish when |λ 1 (p)| = |λ 2 (p)| (of course, it might be three or more eigenvalues of equal moduli; our simplification is sufficient here). This defines a set of curves in the complex plane, where all zeros should accumulate in the n → ∞ limit. The interested reader should refer to the work of Salas and Sokal [START_REF] Salas | Transfer matrices and partition-function zeros for antiferromagnetic Potts models. I. General theory and square-lattice chromatic polynomial[END_REF] for a very detailed discussion of the convergence to the limiting curves. This behavior would not be modified when one or more of the α i 's is a polynomial in n [START_REF] Beraha | Limits of zeros of recursively defined families of polynomials[END_REF].

In the following we limit ourselves to the perfect-nodes case (ρ = 1) and look for the limiting curves for the architectures considered in Section 4, because it is illustrative enough, and because all expressions are then simple and analytical. The trivial solution p = 0 will not always be considered in the following.

Calculation for the symmetrical ladder (S 0 → U n )

Let us first consider the simplest case for Rel 2 (S 0 → U n ) given in eq. ( 31) because only two eigenvalues are present (see Fig. 4). We have to find p such that |x + /x -| = 1, or

x + x - = 1 = 1 + 2 p (1 -p) + 1 + 4 p 2 (1 -p) 2 1 + 2 p (1 -p) -1 + 4 p 2 (1 -p) 2 . ( 49
)
Introducing t = 2 p (1p), we must find t such that

1 + t + √ 1 + t 2 1 + t - √ 1 + t 2 = e i θ (-π ≤ θ ≤ π), (50)
which can be rewritten as

√ 1 + t 2 = i (1 + t) tan θ 2 . ( 51
)
Squaring both sides of eq. ( 51) implies

t 2 + 2 sin 2 θ 2 + 1 = 0, (52)
which yields two possible solutions, t ± , of modulus 1:

t ± = -sin 2 θ 2 ± i 1 -sin 4 θ 2 . ( 53
)
Only the solutions that verify eq. ( 51) must be kept ; they may be rewritten in condensed form

t = -e i ϕ - π 2 ≤ ϕ ≤ π 2 , ( 54
)
which represents a half-circle (the real part of t ± must be negative), so that the searched values of p are finally

p = 1 ± √ 1 + 2 e i ϕ 2 - π 2 ≤ ϕ ≤ π 2 , (55)
It is then straightforward to give the coordinates of some particular points of the curves. On the real axis, the solutions for ϕ = 0 are (1 ± √ 3)/2, i.e., +1.366 and -0.366. The endpoints of the two curves are actually 1 2 (1

± 1+ √ 5
2 ) ± i

) 56
We see on Fig. 6 that for n = 10 already, the zeros are located very close to the two well-separated arcs of quasi-circular shape defined by eq. ( 55). Note that with respect to the diagrams by Chang and Shrock for the all-terminal reliability polynomials calculated for several recursive families of graphs [START_REF] Chang | Reliability polynomials and their asymptotic limits for families of graphs[END_REF], we now have accumulation points in the region Re(p) < 0. When ρ = 1, the global appearance is unchanged, even if the zeros move away from the origin. (p): the dots corresponds to n = 10 and eq. ( 56), the continuous curves to the limit n → ∞ given by eq. ( 55). The trivial solution p = 0 is not displayed.

When ρ approaches zero, the structure distorts itself to look like a circle of radius of the order of (2 ρ) -1/2 centered at ( 1 2 , 0). The convergence of the roots to the limiting curve follows the general expression of Salas and Sokal [START_REF] Salas | Transfer matrices and partition-function zeros for antiferromagnetic Potts models. I. General theory and square-lattice chromatic polynomial[END_REF].

Calculation for

Rel 2 (S 0 → S n ) and Rel 2 (S 0 → T n )
In these configurations, we now have three eigenvalues. We must therefore consider |x ± /x 0 | in addition to |x + /x -|, calculated in the preceding section. We must find all the p's such that

x ± x 0 = 1 + 2 p (1 -p) ± 1 + 4 p 2 (1 -p) 2 2 (1 -p) = e i θ (-π ≤ θ ≤ π), (57)
so that

Z = ± 1 + 4 p 2 (1 -p) 2 = 2 (1 -p) e i θ -(1 + 2 p (1 -p)) . ( 58
)
Here again, we square both sides of eq. ( 58), keeping in mind we will have to check that the real part of the right-hand-side of eq. ( 58) must comply with our determination of the square root (our angle argument being taken between -π and +π, we choose Re √ • • • ≥ 0). After simplification, we get

1 + 2 p (1 -p) = (1 -p) e i θ + p e -i θ . ( 59
)
This equation has two possible solutions, namely

p ± = 1 2 1 + i sin θ ± 1 + (1 -cos θ) 2 . ( 60
)
Which values of θ satisfy with eq. ( 58) ? Using eqs. ( 58) and ( 59), we find

Z = (1 -p ± ) e i θ -p ± e -i θ = ∓ cos θ 1 + (1 -cos θ) 2 + i sin θ (1 -cos θ), (61)
so that Re(Z) and Re(∓ cos θ) have the same sign. Since for x + /x 0 , we must have Re(Z) ≥ 0 and thus Re(∓ cos θ) ≥ 0, we find

p[|x + /x 0 | = 1] =      p + ( π 2 ≤ |θ| ≤ π) p - (0 ≤ |θ| ≤ π 2 ) , (62)
p[|x -/x 0 | = 1] =      p - ( π 2 ≤ |θ| ≤ π) p + (0 ≤ |θ| ≤ π 2 ) . (63)
We also have to check, when considering p[|x + /x 0 | = 1] for instance, that the third eigenvalue x -has a smaller modulus. Consequently, using eqs. ( 59) and (61),

x + x - = (1 -p) e i θ + p e -i θ + ((1 -p) e i θ -p e -i θ ) (1 -p) e i θ + p e -i θ -((1 -p) e i θ -p e -i θ ) = 1 -p p ≥ 1, (64)
and therefore that Re(p) ≤ 1 2 , where p stands for p ± . We deduce from eq. ( 61) that the p + solution is to be discarded. The same procedure for p[|x -/x 0 | = 1] shows that only p -should be considered.

Finally, we must assess the modulus of x 0 along the curves given by eq. ( 55). Here, the situation is slightly more complicated, and we checked numerically that the p -solution of eq. ( 55) should be discarded. The final result

p =        1 + √ 1 + 2 e i ϕ 2 (-π 2 ≤ ϕ ≤ π 2 ) 1 2 1 -1 + (1 -cos θ) 2 + i sin θ (-π ≤ θ ≤ π) (65)
is displayed in Fig. 7. (p): the dots corresponds to n = 10 (eq. ( 66)), the continuous curves to the limit n → ∞ given in eq. ( 65).

When n = 10, we have

We have plotted in Fig. 7 the zeros for n = 10 and the n → ∞ limit. Note that the asymptotic limit seems to be reached faster in the Re(p) > 0 half-plane.

The new curve endpoints of Fig. 7 are 1 2 (1 -√ 2 ± i) ≈ -0.207 ± 0.5 i, and

1 2 (1 - √ 5
) ≈ -0.618. For ρ = 1, the structure of the zeros spreads out from the origin, as in the preceding subsection.

We want to emphasize that while the structure of the zeros in Figs. 6 and7 is similar in the Re(p) > 0 half-plane -the n → ∞ limit is identicalthis is certainly not the case in the Re(p) < 0 half-plane, even though the two graphs differ by two links only. Inferences from these structures of zeros should therefore be considered cautiously.

Sensitivity

An important parameter appearing in reliability studies is the sensitivity s i [START_REF] Rubino | Sensitivity computation in network reliability analysis[END_REF], which measures the influence of a particular equipment i (node or edge) on the global reliability; the sensitivity analysis is akin to the various importance criteria compiled by Henley et Kumamoto [START_REF] Henley | Probabilistic Risk Assessment: Reliability Engineering, Design and Analysis[END_REF]. It is usually defined as the derivative of the total -be it two-or all-terminal -reliability with the individual reliability p i

s i = ∂Rel 2,A ∂p i = Rel 2,A (p i = 1) -Rel 2,A (p i = 0). (67)
The two definitions are equivalent, because Rel 2,A is an affine function of all the p i 's. The calculation can be performed quite easily, since the desired p i appears in only one of the transfer matrices of Section 3. For instance, let us consider the sensitivity of the successive "rungs" of the ladder -a similar work could be done for the nodes, without difficulty. In the following, we restrict ourselves to the case Rel 2 (S 0 → T n ), with perfect nodes and identical edge reliabilities p, and want to assess the influence of the edge location to the overall performance of the network.

From eq. ( 13), taking the derivative of a matrix with respect to a given b n 0 is not difficult, since each of its elements is an affine function of b n 0 . A new matrix appears in the problem, and when all edge reliabilities are set to p (we take ρ = 1 for the sake of simplicity), we have Figure 8 shows that these sensitivities vary with p and that they uniformly decrease from the ladder's ends to its center. Note that only a few coefficients of the polynomials differ. If the length of the ladder increases, the maximum sensitivity will actually decrease in magnitude, and be located at higher values of p. Besides, if n is odd, then the minima are shared by the two central rungs, i.e., n 0 = (n ± 1)/2.

s bn 0 = (0, 1, 0) •        p p 3 p 3 p p p 2 -p 2 -p 3 p 2 (1 -2 p)        n-n 0 •        0 p 2 p 2 1 0 p -p -p 2 -2 p 2        •        p p 3 p 3 p p p 2 -p 2 -p 3 p 2 (1 -2 p)        n 0 •        1 0 0        . (68)
Because of the form of eq. ( 68), the final result is expected to be a sum of different power-law contributions, with products of two of the different eigenvalues already observed for the reliability. The easier way is then to calculate the generating functions in a S 0 → T 2 n ladder, for both b 0 and b n ≡ b central . We finally get G sens (b i ) = n=0 s b i (S 0 → T 2 n ) x n with G sens = N sens /D sens :

N sens (b central ) = 1 -p 2 (1 + p + 3 p 2 -6 p 3 + 2 p 4 ) x + (1 -p) p 5 (1 + 3 p -3 p 2 -2 p 3 + 2 p 4 ) x 2 -p 9 (1 -p) 3 x 3 , (73)
D sens (b central ) = 1 -p 2 (1 -p) 2 x 1 -p 3 (1 -p) x × 1 -p 2 (1 + 2 p + 2 p 2 -8 p 3 + 4 p 4 ) x + p 6 (1 -p) 2 x 2 . ( 74
)
whereas

G sens (b 0 ) = 1 -p 2 (1 + 3 p 2 -5 p 3 + 2 p 4 ) x + (2 -p) p 6 (1 -p) 2 x 2 (1 -p 2 (1 -p) 2 x) 1 -p 2 (1 + 2 p + 2 p 2 -8 p 3 + 4 p 4 ) x + p 6 (1 -p) 2 x 2 .( 75
)
The partial fraction decompositions of G sens (b 0 ) and G sens (b central ) lead to

s b 0 = 1 2 p 2 n    (1 -p) 2 n + x 2 n + 2   1 + 1 -2 p 2 1 + 4 p 2 (1 -p) 2   + x 2 n - 2   1 - 1 -2 p 2 1 + 4 p 2 (1 -p) 2      (76)
and

s b central = 1 2 p 2 n (1 -p) 2 n + 4 p 2 (1 -p) 1 + 4 p 2 (1 -p) 2 p n (1 -p) n + x 2 n + 4   1 + 1 -2 p 2 1 + 4 p 2 (1 -p) 2   2 + x 2 n - 4   1 - 1 -2 p 2 1 + 4 p 2 (1 -p) 2   2      , (77)
so that the ratio s b central /s b 0 , in the large n limit, is given by

s b central s b 0 -→ 1 2   1 + 1 -2 p 2 1 + 4 p 2 (1 -p) 2   , (78)
as displayed in Fig. 9 (for p = 1q close to 1, this ratio expands as 2 q -6 q 3 + 8 q 4 + • • •). Actually, this limit is reached very quickly, since for the S 0 → T 6 ladder, the ratio s b 3 /s b 0 is barely distinguishable from the asymptotic limit in the entire range 0 ≤ p ≤ 1, which is displayed in Fig. 9. It also shows that the overall reliability is more sensitive to "access" rungs located at the edges of the ladder network than to those in the middle, even if all edge reliabilities are identical. We expect the all-terminal reliability R n for the simple ladder diagram of Fig. 1 should exhibit the same behavior, with a transfer matrix depending on a n , b n , and c n , and a generating function that is a rational fraction of p (the common reliability of links) and x. Actually, the calculation of the all-terminal reliability may be somewhat easier because all nodes can be considered perfect without loss of generality [START_REF] Colbourn | The combinatorics of network reliability[END_REF]. We can then use the usual factoring decomposition to establish a relationship between R n and R n-1 . From Fig. 10, we easily find that

R n = [(1 -a n ) b n + a n (1 -b n )] c n R n-1 + a n b n R n-1 (b n-1 → b n-1 // c n ).( 79
)
Using the ansatz

R n = α n (a n + b n ) + β n a n b n , we get    α n β n    =    (a n-1 + b n-1 ) c n a n-1 b n-1 c n (a n-1 + b n-1 ) (1 -2 c n ) + c n (1 -2 b n-1 ) a n-1 (c n + b n-1 (1 -3 c n ))       α n-1 β n-1    , ( 8
where actually α 0 = 1 and β 0 = 0. Calling M n the matrix appearing in eq. ( 80), we obtain (with the condition a 0 = 0 in M 1 )

R n = (a n + b n , a n b n ) • M n • M n-1 • • • M 1 •    1 0    . (81)
When all the reliabilities are equal to p, the transfer matrix is simply

M R =    2 p 2 p 3 p (3 -5 p) p 2 (2 -3 p)    . (82)
The eigenvalues of M R ,

ζ ± = p 2 4 -3 p ± √ 12 -20 p + 9 p 2 2
, are easily obtained, and finally

R n = 1 p √ 12 -20 p + 9 p 2 ζ n+1 + -ζ n+1 - ( 83
)
and

G R (x) = p 1 -p 2 (4 -3 p) x + p 4 (1 -p) x 2 . ( 84
)
The last results were obtained in [START_REF] Chang | Reliability polynomials and their asymptotic limits for families of graphs[END_REF]: the polynomials aggregate at the curve defined by p = 1 + 1 3 e i ϕ , with cos ϕ ≥ 1 3 . It is worth stressing that the new scale ζ + for the all-terminal reliability is different from the eigenvalues found in the two-terminal study, but that a power-law behavior is still the rule.

Conclusion and perspectives

We have found the solution of the two-and all-terminal reliabilities for a simple ladder graph, which corresponds nonetheless to realistic network architectures, especially in telecommunication networks for IP transport, but not only. Node and edge failures are put on an equal footing, and the simple formulae relying on transfer matrices may be directly implemented, even in a worksheet application. When identical reliabilities p and ρ are assumed for edges and nodes, respectively, we have given the analytical solution of the two-terminal reliability, for which only numerical determinations were previously available. We have shown that, while it may be useful to study the location of the zeros of the two-terminal reliability polynomial, it may strongly depend on the existence of a limited number of particular edges, especially in the Re(p) < 0 region. We have also given the generating functions of the reliability polynomials, and provided a glimpse of sensitivity studies, which could be performed straightforwardly.

Although we used a delta-star transformation to solve a particular family of graphs, the present work may clearly be extended in several directions, which we outline in the following, and which we shall develop elsewhere.

Brecht-Colbourn ladder

Another simple ladder graph is the Brecht-Colbourn ladder [START_REF] Brecht | Improving reliability bounds in computer networks[END_REF][START_REF] Prékopa | The use of binomial moments for bounding network reliability[END_REF], which has been considered as a case study, both for the two-and all-terminal reliability, in order to evaluate the quality of bounds for the reliability polynomial in the case of perfect nodes. It has been also brought forward in a special case of broadcasting network [START_REF] Graver | You may rely on the reliability polynomial for much more than you might think[END_REF]. Using the same delta-star transformation, the application of which is slightly more tricky, we have been able to solve this problem exactly, where the intermediate nodes have a connectivity degree of 4 instead of 3. The complete results -the transfer matrix is then of dimension four for imperfect nodes, but still of dimension three for perfect nodes; a few eigenvalues may be complex -will be given elsewhere [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of the Brecht-Colbourn ladder and the generalized fan[END_REF], along with results for a generalized fan [START_REF] Aggarwal | A simple method for reliability evaluation of a communication system[END_REF][START_REF] Neufeld | The most reliable series-parallel networks[END_REF].

More general recursive families of graphs

Even though the delta-star transformation may not apply successfully to all networks, it seems quite clear that a similar decomposition through transfer matrices of two-and all-terminal reliabilities should exist too for a ladder of K 4 graphs, or ladders of greater width. In order to make such calculations useful for applications, imperfect nodes as well as imperfect edges should be considered. All is needed is a recursion relation between successive graphs, when one "elementary cell" is added. This implies a new expression for the deletion-contraction theorem, in which the -most useful -linearity with respect to all individual edge or node reliabilities must be preserved, whereas to our knowledge the edge reliabilities have been renormalized to account for the unreliability of the nodes they connect [START_REF] Theologou | Factoring & reductions for networks with imperfect vertices[END_REF][START_REF] Torrieri | Calculation of node-pair reliability in large networks with unreliable nodes[END_REF]. This new expression will be given elsewhere [START_REF] Tanguy | Transfer matrix method for two-and all-terminal reliabilities in ladder[END_REF], along with an application to other recursive families of graphs (K 4 ladders and K 3 cylinders).

What should we expect ? Basically the same kind of behavior as detailed in the present work, with a factorization of the reliability in terms of transfer matrices, the dimension of which may substantially increase to reflect the interplay of different edges/nodes in the overall reliability, and the underlying algebraic structure of the graph. It should be clear, however, that all graphs might not behave that way [START_REF] Biggs | Recursive families of graphs[END_REF]. For instance, the generating function of the two-terminal reliability for the complete graphs K n , for which we give a recursion relation in the case of imperfect nodes and edges in the appendix, might not be a rational fraction.

Failure frequency of systems

Another performance index for networks, mostly investigated in the context of power distribution systems [START_REF] Hayashi | System failure-frequency analysis using a differential operator[END_REF][START_REF] Shi | General formulas for calculating the steady-state frequency of system failure[END_REF][START_REF] Singh | A new method to determine the failure-frequency of a complex system[END_REF], is the failure frequency of a connection. Calculations can readily be performed when the assumption is made that all equipments fail at a constant failure rate λ e and brought back to operation via a constant repair rate µ e (and p e = µ e µ e + λ e

). The failure frequency is then the sum over all equipments of the products of the individual failure rate, reliability and sensitivity (as defined in Section 8). The present results clearly indicate that such calculations are made very simple by the transfer matrix formulation, even for extended networks. We shall present various examples in a forthcoming paper [START_REF] Tanguy | Exact expressions for the failure frequency of some ladder networks[END_REF].

Reliability incertitude

In recent years, Coit and collaborators [START_REF] Coit | System-reliability confidence-intervals for complex-systems with estimated component-reliability[END_REF][START_REF] Coit | System optimization with component reliability estimation uncertainty: a multi-criteria approach[END_REF][START_REF] Jin | Variance of system-reliability estimates with arbitrarily repeated components[END_REF] have considered the possible influence of the reliability uncertainty to the determination of an "optimal network architecture", and have mostly dealt with series-parallel architectures.

In the context of meshed networks, and recalling our assumption of statistically independent failures for the network constituents, we could safely replace in the final expression of the reliability each p e or p n by the random variable P e or P n of which it is the mean value. We can thus straightforwardly assess the influence of the full statistical properties of P e or P n (variance, skewness, kurtosis, etc.). This issue will be reported elsewhere. We can also replace, for non-reparable systems, each p e by the often used exp(-λ e t).

Bounds, combinatorics, etc.

There are obviously many directions in which this work may be useful in the determination of bounds. The first one is to use our present results on simple ladder graphs, or their future extensions [START_REF] Tanguy | Exact solutions for the two-and all-terminal reliabilities of the Brecht-Colbourn ladder and the generalized fan[END_REF][START_REF] Tanguy | Transfer matrix method for two-and all-terminal reliabilities in ladder[END_REF], as possible upper or lower bounds to more complex graphs. If a particular graph looks like a special instance of a recursive family of graphs, we expect the generating function to be a rational fraction again. The dimension of the corresponding transfer matrix may be probed by trying to find recursion relations between successive reliability polynomials. Of course, if the dimension of the corresponding transfer matrix is large, the degree of the numerator and denominator of the fraction may be too large for a complete solution to be obtained easily. Even so, the knowledge that the generating function is rational may be an indication that an approximate generating function might be quite useful, since Padé approximants [START_REF] Baker | Padé approximants[END_REF] are known for their devilish knack of getting very close to the exact function [START_REF] Press | Numerical recipes in C: The art of scientific computing[END_REF]. This suggests that we may profitably consider rational fractions -deduced from approximate generating functions -for bounds, instead of exclusively relying on polynomials. This will be demonstrated in detail elsewhere [START_REF] Tanguy | Padé approximants for network reliability[END_REF].

Our calculations may also provide some information on some combinatorial issues, such as the enumeration of self-avoiding walks on lattices of restricted width. Finally, the exact results found for classes of arbitrarily large networks may prove useful for testing different algorithms (Monte Carlo, genetic, OBDD, etc.) in numerically exacting configurations, where edge and node unreliabilities have to be taken into account.

A Two-terminal reliability for a complete graph with imperfect nodes

The exact all-terminal reliability polynomial A n = Rel A (K n ) for the complete graph K n with perfect nodes and identical edges has been known for decades, as well as the corresponding two-terminal reliability polynomial T n = Rel 2 (K n ) ( [START_REF] Colbourn | The combinatorics of network reliability[END_REF], pp. [START_REF] Hänsler | Exact Calculation of Computer Network Reliability[END_REF][START_REF] Hänsler | Comments on "A fast recursive algorithm to calculate the reliability of a communication network[END_REF]. They may be found by recurrence, starting with A 1 = 1:

A n = 1 - n-1 j=1    n -1 j -1    A j (1 -p) j (n-j) (A.1)
and

T n = n j=2    n -2 j -2    A j (1 -p) j (n-j) . (A.2)
Taking imperfect nodes with identical reliability ρ into account does not substantially change the all-terminal reliability, since one merely has to multiply the result of eq. (A.1) by ρ to the power of the number of operating nodes [START_REF] Colbourn | The combinatorics of network reliability[END_REF]. The situation is quite different for K-terminal reliability, and in particular for two-terminal reliability, since all the possible paths between the source and the destination may visit different numbers of nodes. This can be very useful since the complete graph configuration represents the highest possible upper bound, because all its nodes are connected.

Defining T n (p, ρ) the two-terminal reliability polynomial with imperfect edges and nodes -so as to keep the notation of Colbourn [START_REF] Colbourn | The combinatorics of network reliability[END_REF], we have

Fig. 1 .

 1 Fig. 1. Different source-terminal connections for lattice graphs.

Fig. 2 .

 2 Fig. 2. Triangle-star transformation for unreliable nodes. A, B, and C are the node reliabilities, with a, b, and c the edge reliabilities of the initial network, and p A , p B , p C , and O those of the transformed network.

Fig. 3 .

 3 Fig. 3. Two-step decomposition of the ladder network from (a) to (c), using the triangle-star transformation.

Fig. 4 .

 4 Fig. 4. Symmetrical ladder network.

Fig. 5 .

 5 Fig. 5. Dependence of λ 0 , λ + , and λ -on p, when ρ = 1.

Fig. 6 .

 6 Fig. 6. Location of the zeros of Rel (S 0 →Un) 2

Fig. 7 .

 7 Fig. 7. Location of the zeros of Rel (S 0 →Tn) 2

ForFig. 8 .

 8 Fig. 8. Variation of the sensitivity Sens 2 of the ladder's rungs, depending on their location, in the Rel 2 (S 0 → T 6 ) case.

Fig. 9 .

 9 Fig. 9. Ratio s b central /s b 0 for Rel 2 (S 0 → T 2n ) in the large n limit.

Fig. 10 .

 10 Fig. 10. Successive deleting/contracting operations lead to a recursion relation between R n and R n-1 (see text).

T 2 2 + 2 p 2 ρ 3 -7 p 4 ρ 4 + 7 p 5 ρ 4 -2 p 6 ρ 4 + p 3 -2 ρ 3 + 2 ρ 4 ,P 1 

 2234443341 (p, ρ) = p ρ 2 , T 3 (p, ρ) = p ρ 2 + p 2 ρ 3p 3 ρ 3 , (A.3) T 4 (p, ρ) = p ρ (p) = p, P 2 (p) = p 2p 3, and P 3 (p) = 2 p 3 -7 p 4 + 7 p 5 -2 p 6 , etc. We can then invert eq. (A.4) given for T n (p, ρ = 1) to obtain all the P j (p)'s, using the inversion formula P j (p) = T k+1 (p, ρ = 1).(A.5)Combining Eqs. (A.4) and (A.5) gives an easily computable expression for T n (p, ρ).

  Rel 2 (S 0 → T 10 ) = p 11 (11 + 155 p 2 -99 p 3 + 40 p 4 -907 p 5 -296 p 6 +1448 p 7 + 3121 p 8 -1102 p 9 -7989 p 10 -1747 p 11 +14806 p 12 + 4776 p 13 -24168 p 14 + 176 p 15 + 35072 p 16 -38016 p 17 + 19072 p 18 -4864 p 19 + 512 p 20 ).

All-terminal reliabilityAs mentioned in the introduction, the all-terminal reliability is another useful measure of the network availability, namely the probability that all nodes are connected. We refer the reader to the work of Chang and Shrock[START_REF] Chang | Reliability polynomials and their asymptotic limits for families of graphs[END_REF] for the explicit expressions of Rel A for various recursive families of graphs, among which our simple ladder graph, when all edges have the same reliability p. For the sake of completeness, we show that their result may be slightly generalized for edges with distinct reliabilities. Here again, the final, analytical expression can be written in a concise form using transfer matrices. In the context of graph theory, this can be viewed as the factorization of a particular value of the multi-variate Tutte polynomial considered by Wu[START_REF] Wu | Graph theory in statistical physics[END_REF] and Sokal[START_REF] Sokal | The multivariate Tutte polynomial (alias Potts model) for graphs and matroids[END_REF].

Acknowledgments

I am indebted to Élisabeth Didelet for giving me the opportunity to study reliability issues in depth, and to Joël Le Meur, François Gallant, Veluppillai Chandrakumar, Adam Ouorou, and Éric Gourdin for encouragement and support. Gérard Cohen provided important information. I would also like to thank Christine Leroy, Wojtek Bigos, Annalisa Morea, Antoine Pardigon, and Annie Druault-Vicard for useful discussions.


David Uystepruyst
email: david.uystepruyst@univ-valenciennes.fr

François Monnoyer

A numerical study of the evolution of the blast wave shape in rectangular tunnels

When the explosion of condensed materials occurs in square or circular cross-section tunnel, the subsequent blast wave reveals two patterns: threedimensional close to the explosive charge and one-dimensional far from the explosion. Pressure decays for these two patterns have been thoroughly studied. However, when the explosion occurs in rectangular cross-section tunnel, which is the most regular geometry for underground networks, the blast wave exhibits a third, two-dimensional, patterns. In order to assess the range of these three patterns, several numerical simulation of blast waves were carried out varying the width and the height of the rectangular cross-section as well as the mass of the charge. Laws are presented to localize the transition zones between the 3D and the 2D patterns, and between the 2D and the 1D patterns, as functions of non-dimensional width and height. The numerical results of the overpressure are compared to existing 3D and 1D

laws. An overpressure decay law is proposed to represent the 2D pattern.

Knowing the two transition zones and the overpressure decays within these

zones, an algorithm is presented to efficiently predict an overpressure map. This algorithm is validated by comparison with experimental data. Blast wave is a topical research subject as shown by the several recent works [START_REF] Igra | Review of methods to attenuate shock/blast waves[END_REF][START_REF] Langlet | Transient response of a plate-liquid system under an aerial detonation: Simulations and experiments[END_REF][START_REF] Kim | Effect of propagation behaviour of expanding spherical flames on the blast wave generated during unconfined gas explosions[END_REF], and in particular confined blast wave [START_REF] Song | Dynamic response of composite shell under axial explosion impact load in tunnel[END_REF][START_REF] Buonsanti | 3-D simulation of tunnel structures under blast loading[END_REF]. One of the application of the confined domain is the underground network. Underground The first pattern cited above, ie. the free-field pattern, is indubitably the most studied. From these studies, scaling laws were derived, as the laws from Baker, Cox, Westine, Kulesz and Strehlow [START_REF] Baker | Fundamental studies in engineering[END_REF], [START_REF] Mills | The design of concrete structures to resist explosions and weapons effects[END_REF], [START_REF] Brode | Blast waves from a spherical charge[END_REF] or [START_REF] Henrych | The dynamics of explosion and its use[END_REF] relating the maximum overpressure peak to the distance from the explosive charge.

In fact, [START_REF] Henrych | The dynamics of explosion and its use[END_REF] proposed one of the most common free-field decay law, which is expressed as:

             ∆P max P ref Henrych = 14.072 Z + 5.54 Z 2 - 0.357 Z 3 + 0.00625 Z 4 if 0.05 ≤ Z ≤ 0.3 ∆P max P ref Henrych = 6.194 Z - 0.326 Z 2 + 2.132 Z 3 if 0.3 ≤ Z ≤ 1 ∆P max P ref Henrych = 0.662 Z + 4.05 Z 2 + 3.288 Z 3 if 1 ≤ Z . (1)
The third propagation pattern is the case of a confined explosion, which considers that the blast wave propagates inside a confined space that is strong enough to withstand the explosive charge impulse (e.g., a tunnel). Among the few reported experiments investigating air detonation in underground environments, some gave overpressure decays laws during this third pattern. [START_REF] Curran | Underground storage of ammunition -experiments concerning accidental detonation in an underground chamber[END_REF] proposed the following pressure-distance law for various explosives weights:

∆p p 0 = M Φ 2 x 0.8 . (2)
[17] determined the following overpressure-distance decay relationship:

∆p p 0 = 7.028 M A r 0.514 . (3)
Applying the energy concentration concept (ECF), which is detailed in section 2.3, [START_REF] Silvestrini | Energy concentration factor. A simple concept for the prediction of blast propagation in partially confined geometries[END_REF] the very similar law: The objective of this paper is to determine the position of both the 3D-2D and the 2D-1D transition zones and to predict the overpressure occurring in these three zones. To accomplish this, the detonation of different quantities of TNT explosives inside a perfectly rigid tunnel was simulated.

∆p p 0 = 7.43538 M A r 0.51 . (4)
The following section presents the geometrical configuration, the numerical methodology and the ECF method. Moreover, scaling laws that eliminates the solution's parametric dependence on the explosive energy, the weight of the explosives and the real tunnel cross-sectional size are provided. These scaling laws transform the infinite number of solutions into a single solution that demonstrates a monotonic transition from one wave pattern to another pattern. Afterwards, the numerical results are brought in section 3. These numerical results consisted of overpressure history for different widths and heights of the cross-section, and of different detonation mass. In section 4, fitting power law are proposed to estimate the locations of both transition zones. Furthermore, an algorithm giving the complete map of the overpressure pattern is introduced. This algorithm is validated by comparison with previous experimental data.

Configurations, numerical details and ECF method

Calculation domains and non-dimensional numbers

The calculation domains are presented in figure 1. They consist on rectangular cross-sectional tunnel of length L = 25 m. The cross-sectional widths and heights ranged from 3 to 7 meters and 2 to 5 meters respectively. Especially, a square cross-sectional tunnel of 3 m 2 area was considered to verify results from [START_REF] Benselama | A numerical study of the evolution of the blast wave shape in tunnels[END_REF]. In order to study the effect of the width, four additional configurations with 3 m height and width ranging from 4 to 7 m are investigated. A 6 m width and 2 m height domain is considered because it yields the same hydraulic diameter as the 3 m side square section. Then, two additional configurations with 6 m width and respectively 4 and 5 m height are studied to evaluate the effect of the height. Furthermore, TNT charge ranging from 1.000 to 10.000 kg are investigated.

It can be anticipated that the first reflection modifying the 3D shape of the blast wave in a 2D shape depends on the height of the domain H.

Therefore a first parameter α H , inspired by the work of [START_REF] Benselama | A numerical study of the evolution of the blast wave shape in tunnels[END_REF] is considered: α H = 100

d c H . (6)
In the same way, the difference between the second and the first reflection modifying the 2D shape of the bast wave in a 1D shape depends on difference between the width, W , and the height of the domain. Thus, a second parameter is:

α W -H = 100 d c W -H . (7)
The parameter α W based on the width of the domain will also be used for convenience.

The numerical results are highlighted by the pressure coefficient defined as:

C p = p -p 0 p 0 . (8)

Numerical details

1D calculation was performed as long as the incident blast wave did not encounter any obstacle. This strategy obviously yield a faster and more accurate computations. During this 1D calculation, the TNT gases, composed of hot detonation products, obey the Jones-Wilkins-Lee law, which is expressed by [START_REF] Mader | Numerical modeling of explosives and propellants[END_REF]:

p = A 1 - ω R 1 ρ ρ c exp(-R 1 ρ c /ρ) + B 1 - ω R 2 ρ ρ c exp(-R 2 ρ c /ρ) + ωρe, (9)
where the parameters A, B, R 1 , R 2 , ω and ρ c depend on the explosive material. Table [START_REF] Baker | Fundamental studies in engineering[END_REF] provides these parameters. Details of the 1D solver, as well as the 1D-3D transfer method are available in [START_REF] Benselama | A 1D-3D mixed method for the numerical simulation of blast waves in confined geometries[END_REF]. The same 1D calculation was used for all configurations.

The initial condition for the 3D calculation, thus the result of the 1D calculation, was located in the center of the cross-section, see figure 1, using a symmetry condition in the longitudinal axis. which were solved by a proprietary software [START_REF] Benselama | A 1D-3D mixed method for the numerical simulation of blast waves in confined geometries[END_REF]. The numerical method on which this software's solver is based is a unstructured finite-volume cellcentered approach using the traditional upwind scheme and a two-stage explicit time integration technique, yielding a second-order accuracy in both space and time. In order to prevent numerical oscillations, which may occur in regions with strong gradients, the Total Variation Diminishing (TVD) minmod scheme is used [START_REF] Van Leer | Towards the ultimate conservative scheme, V. A second-order sequel to Godunov's method[END_REF]. The spatial discretization is performed with an automatic Cartesian grid generator [START_REF] Deister | Selforganizing hybrid cartesian grid generation and application to external and internal flows, numerical flow simulation III[END_REF]. In addition, the Courant-Friedrichs-Lewy (CFL) condition has to be satisfied in order to guarantee the stability of the time integration technique. The successive gauges, also centered into the cross-section, were located 0.2 m from each other knowing that the first one was 1 m from the explosive charge.

The initial computational domains are split into 0.5 m cube elements and a dynamic adaptive mesh method was applied to finely catch the incident wave: the refinement criterion was based on the density value. A fine mesh was kept after the incident wave passage to guarantee an accurate prediction of the complex waves pattern following the incident wave. The averaged number of elements was ranged from 4 × 10 6 for the smallest domain to 10 7

for the largest one. The smallest refined elements were 1.5 × 10 -2 m ensuring a mesh wave number λ of 67, 106, 128 and 144 for the TNT charges of 1.000, 4.000, 7.000 and 10.000 kg, respectively. The mesh wave number has to satisfy in the ideal case λ > 100 [START_REF] Benselama | A 1D-3D mixed method for the numerical simulation of blast waves in confined geometries[END_REF].

The concept of Energy Concentration Factor

The concept of energy concentration factor introduced by Silvestrini et al.

[16] is adapted in the present paper. The idea is to scale the abscissa of freefield decay law, for instance Henrych's law [START_REF] Baker | Fundamental studies in engineering[END_REF], by geometrical consideration.

As shown in figure 2, for the same volume the cylinder shape blast wave (2D) has a radius R ′ defined as R ′ = 4R 3 /3W where R is the radius of the equivalent spherical shape blast wave (3D). Then, the new scale distance is:

Z ′ = R ′ M 1/3 c = 4Z 3 M 1/3 c 3W .
It has to be noted that the ECF method was essentially developed and validated for 3D to 1D transition. Therefore, one objective of the present paper is to see whether the ECF method is available for 3D to 2D transition.

Results and discussion

3.1. Effects of the width, α H = 3.55, M c = 1.000 kg However, the 2D pattern of the configuration α H = 5.33 seems to go towards the 2D pattern of the configuration α H = 3.55 before Z = 7 m/kg 1/3 . In the same way, the 2D pattern of the configuration α H = 3.55 seems to go towards the 2D pattern of the configuration α H = 2.66. These two observations yield the conclusion that the overpressure peak does not depend on the height. The 1D law developed from ECF method well predicts the 1D pattern of α H = 11.15. However, figure 5 shows the disability of the ECF method to predict the overpressure peak during the 2D behavior of the wave when the first reflection is occurring quickly, e.g. for a low value of Z.

Effects of the charge of the mass

Z (m/kg The results of the previous sections have shown the ECF method is particularly accurate for the 1D pattern of blast waves. It has been also shown that this method is not suitable for the 2D pattern especially for high value of α H . However, it has been proved that the 2D pattern of the incident blast wave is the same for every configurations. Therefore, a polynomial interpolation, as the one of equation ( 1) can be proposed. Finally, knowing the transition locations, Z H and Z W , the complete incident wave pattern can be deduced.

Figure 6 shows the distributions of the discontinuity locations Z H and Z W -H = Z W -Z H versus α H and α W -H , respectively. It is shown both discontinuity locations depend linearly on the parameter α. Using the leastsquare algorithm, best correlations are gave by:

Z H = 0.07992 (α H /100) 1.317 , (10)
and

Z W -H = 0.07994 (α W -H /100) 1.351 . (11)
The first fit has a correlation coefficient of 0.999, while the second one has a correlation factor of 0.971. Actually, it even sounds possible to use a common power fit for both discontinuity locations. Finally, it has to be said that these laws are close to the law giving the transitional zone 3D-1D (5)

found in [START_REF] Benselama | A numerical study of the evolution of the blast wave shape in tunnels[END_REF] for square sections.

Interpolating the clearest points of 2D patterns superimposition gives the 2D law:

(C p ) 2D = 33.78 Z 3 + 3.74 Z 2 + 1.52 Z . (12)
Finally the algorithm of calculation of the overpressure peak can, then, be written as:

1. calculate Z H = 0.07992 (α H /100) 2. for Z ≤ Z H , the pressure is given by (1), 3. calculate Z W = Z H + 0.07994 (α W -H /100) 1.351 , 4. for Z H < Z ≤ Z W , the pressure is given by (12),

for

Z > Z W , Z ′ = 3ZW H 2π 1/3 , 6
. for Z ′ , the pressure is given by (1).

Validation of the algorithm

Validation is done by comparison with experimental data from a series Based on the maximum pressure effect [START_REF] Formby | Blast characteristic and TNT equivalence values for some commercial explosives detonated at ground level[END_REF][START_REF] Wharton | Air blast TNT equivalence for a range of commercial blasting explosives[END_REF], the TNT equivalence factor was found to be 1.18. In the experimental setup, 5 sensors located every 2 m from the charge were considered.

Figure 7 shows the peak of overpressure of the incident wave for the experimental data and the predicted map done by the algorithm. It is shown the predicted pressure is in good convenience with experimental data. The transition zones are correctly determined by the algorithm. Moreover, the values of experimental pressure peak are well fitted by the predicted pressure.

Conclusion

Numerical computations are carried out to predict the locations of the 3D-2D transition and the 2D-1D transition of a blast wave in a rectangular cross-sectional tunnel. In order to avoid any lack of generality, different width and height domains as well as different charge of mass were investigated. The numerical results were compared to free-field law and to the ECF method for the 2D and the 1D patterns.

The main results of this work are:

• the 3D-2D transition is only a function of α H

• as soon as the wave became 2D, the transition to the 1D pattern is only a function of α W -H

• the ECF method is suitable to predict the overpressure of the 2D pattern for the lower value of α H

• a polynomial function is proposed to predict the overpressure peak of the 2D pattern to provide the lacks of the ECF method for high value of α H

• the ECF method is suitable to predict the overpressure of the 1D pattern An algorithm has been proposed to efficiently predict the overpressure pattern. This algorithm has been validated by comparison with experimental data. This algorithm may be very helpful to efficiently drawing up the blast wave damages.

Keywords: 3 c 2 (u 2 + v 2 + w 2 )

 32222 Damage assessment for explosions, three-dimensional simulation, blast wave, confined domain, adaptive grid Nomenclature Greek letters α x = 100d c /x ratio of the explosive diameter d c to the tunnel hydraulic diameter d H , height H, width H or W -H ∆x mesh cell size (m) λ mesh wavenumber, λ = m 1//∆x (kg 1/3 /m) ω constant for Jones-Wilkins-Lee (JWL) equation of state Φ tunnel cross-section diameter (m) ρ density (kg/m 3 ) Latin letters A tunnel cross-sectional area (m 2 ) A, B, C, R 1 , R 2 constants for the JWL equation of state C p pressure coefficient C p = (pp 0 )/p 0 d diameter (m) E total specific internal energy, E = e + 1 zones between 3D and 2D patterns (x = H), and between 2D and 1D patterns (x = W ). Z W -H = Z W -Z H (m/kg 1/3

  networks such as subway station have rectangular cross-section. Typically, the height is smaller than the width. Therefore, blast waves occurring in such a domain exhibit three patterns: (a) a free-field pattern, known to yield fast overpressure decay, while the blast wave does not reach any obstacle; (b) a two-dimensional (2D) pattern after the first reflection (vertical: in the direction of the height) and (c) a one-dimensional (1D) pattern after the second reflection (transversal: in the direction of the width). 2D and 1D patterns obviously induced lower overpressure decays involving more dramatical damages not only for the structures but also for the peoples. The knowledge of the global behavior of blast waves in such a confined domain is thus decisive for safety reasons.

[ 3 ]

 3 have investigated the position of the transition zone between the three-dimensional (3D) and the 1D patterns for a square cross-sectional tunnel. They performed detonations of charges ranging between 0.150 kg and 15.000 kg of trinitrotoluene (TNT) in a 5 m 2 cross-sectional tunnel. Defining the ratio size α = 100d c /d h , where d c was the charge diameter and d h was the hydraulic diameter of the tunnel. They shown the transition zone is located at : square cross-sectional tunnel the four Mach reflections (upper and lower in both transversal and vertical coordinates) occur at the same time producing a large pressure increase when the Mach reflections catch up the incident wave. The shape of the blast wave looks totally different in rectangular cross-sectional tunnel where transversal and vertical reflections do not occur simultaneously.

Figure 1 :

 1 Figure 1: Calculation domains.

Figure 3

 3 Figure 3 shows the distribution of the overpressure peak of the incident wave for α H = 3.55 (typically H = 3 m and M c =1.000 kg). The square

Figure 2 :

 2 Figure 2: Scaling used for the ECF method.

Figure 3 :

 3 Figure 3: Distribution of the maximum pressure over the axial gauges for different domain height, α H = 3.55, M c =1.000 kg. 2D and 1D laws were calculated on α W -H = ∞.

Figure 5

 5 Figure5shows the distribution of the overpressure peak of the incident wave for mass of the charge ranging from 1.000 kg to 10.000 kg. Every calculation were done on the same domain, W = 6 m and H = 2 m. The numerical results are compared with 3D law, and 2D and 1D laws calculated with the ECF method on α H = 11.15 and α W -H = 5.73, ie. for a charge of 10.000 kg.The main information in figure5is that 2D patterns are the same for every mass. Configurations α H = 11.15, α H = 10.19 and α H = 8.45 clearly have overpressure peaks superimposed between both jumps. Furthermore, the 2D pattern of α H = 5.33 seems to be the continuity of the three configurations

Figure 4 :

 4 Figure 4: Distribution of the maximum pressure over the axial gauges for different domain height, α W = 1.77, M c = 11.000 kg. 2D and 1D laws were calculated on α H = 5.33.

Figure 5 :

 5 Figure 5: Distribution of the maximum pressure over the axial gauges for different mass of the charge. 2D and 1D laws were calculated with α H = 11.15 and α W -H = 5.73.

4 .

 4 Determination of Z H , Z W and the overpressure peak value in 2D and 1D

  of trials carried out by INERIS (French Institute of Industrial Environment and Risks) and already presented in[START_REF] Benselama | A numerical study of the evolution of the blast wave shape in tunnels[END_REF]. The trial used in the present pa-

Figure 7 :

 7 Figure 7: Comparison between experimental data and predicted value of pressure peak.

Table 1 :

 1 Parameters used for the JWL law describing the T.N.T. explosive material

	Specific	Density,	A	B	R 1	R 2	ω
	energy, e c	ρ c
	kJ/kg	kg/m 3	10 11 P a 10 9 P a	-	-	-
	4870	1580	3.73	3.74	4.15	0.90	0.35
	Blast wave propagation is governed by the unsteady Euler equations,

  1.317 ,

		Z H =f( H ) numerical results Z H =0.0799/( H /100) 1.317
		Z W-H =f( W-H ) numerical results Z W-H =0.0799/( W-H /100) 1.351
	10	1
	Z H , Z W-H
	10 0
		5	10	15 20 25
		w , h-w

Figure 6: Discontinuity locations Z H and Z W -H as functions of α H and α W -H .

Acknowledgments

This research was carried out for the TRISTAN project, supported by the French National Research Agency (ANR).

List of Tables

Parameters used for the JWL law describing the T.N.