Benoit Vereecke
email: vereecke@lmt.ens-cachan.fr

Henri Bavestrello
email: bavestre@lmt.ens-cachan.fr

D Dureisseix
email: david.dureisseix@lmgc.univ-montp2.fr

An extension of the FETI domain decomposition method for incompressible and nearly incompressible problems

Keywords: Stokes, Herrmann, multilevel domain decomposition, scalability, mixed formulation

come

An extension of the FETI domain decomposition method for incompressible and nearly incompressible problems

Introduction

Incompressible or nearly incompressible problems arise, for instance, from fluid flows, polymer injection, casting, extruding, spinning, or from the simulation of elastomer structures, etc. Because low order divergence-free displacement discretizations with finite elements may lead to severe locking problems (see [START_REF] Babuska | Locking effects in the finite element approximation of elasticity problems[END_REF] for instance), the treatment of such problems is usually performed with an approximation of the incompressibility (or divergence-free) condition. For instance, one may use the penalty method, or equivalently for the case of elasticity, a Poisson's ratio close to 0.5; but these methods lead to ill-conditioned systems of equations. When the size of the problem increases, this can prevent the use 1 of iterative methods, though they are often less CPU expensive than direct approaches, [START_REF] Bramley | Iterative methods and formulations for incompressible FEM flow solvers[END_REF][START_REF] Bauer | Performance of parallel preconditioners for adaptive hp FEM discretization of incompressible flows[END_REF]. Furthermore, the precision obtained for the pressure estimation is poor. Other techniques use a conforming or non-conforming mixed finite element (displacement -pressure formulation), with or without under-integration [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Stenberg | Mixed hp finite element methods for problems in elasticity and Stokes flow[END_REF].

In this paper, we are dealing with mixed formulation and pressure -displacement (or celerity for Stokes problem) coupled problems, within the context of large-scale finite element simulations. In this situation, domain decomposition methods are efficient computational strategies well suited to parallel architecture computers, especially multilevel approaches, see [START_REF] Tallec | Domain decomposition methods in computational mechanics[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF]Recent Developments in Domain Decomposition Methods[END_REF] for instance. Several multilevel domain decomposition methods have already been applied to incompressible problems. Concerning overlapping Schwarz methods, the reader can refer to [START_REF] Klawonn | A comparison of overlapping Schwarz methods and block preconditioners for saddle point problems[END_REF]. For non-overlapping Schur methods, the primal BDD method [START_REF] Mandel | Balancing domain decomposition[END_REF] (for Balancing Domain Decomposition) has been extended: see for instance [START_REF] Cowsar | Balancing domain decomposition for mixed finite elements[END_REF][START_REF] Tallec | Non-overlapping domain decomposition methods for adaptive hp approximations of the Stokes problem with discontinuous pressure fields[END_REF][START_REF] Pavarino | Balancing Neumann-Neumann methods for incompressible Stokes equations[END_REF], but less has been done for dual methods. Usually, the incompressible condition is treated iteratively with a Uzawa iteration (i.e. with a physical partitioning of the unknowns), or a pressure correction algorithm is used. In these cases, the inner-loop linear systems are solved with classical domain decomposition; see [START_REF] Léné | Some strategies to compute elastomeric lamified composite structures[END_REF] and [START_REF] Grant | Coarse grain parallel finite element simulations for incompressible flows[END_REF][START_REF] Axelsson | Solving the Stokes problem on a massively parallel computer[END_REF] for an example of each of these approaches. A monolevel dual Schur method has been used in [START_REF] Calgaro | On the domain decomposition method for the generalized Stokes problem with continuous pressure[END_REF], but the extension to a multilevel scheme has not been addressed. Concerning other monolevel Schur methods for Stokes problem, one may refer to the early work [START_REF] Bramble | A domain decomposition technique for Stokes problems[END_REF] or [START_REF] Ainsworth | Domain decomposition preconditioners for p and hp finite element approximations of Stokes equations[END_REF]. In all of these cases, a monolevel scheme does not possess the numerical scalability property.

We are using herein a dual and multilevel domain decomposition methodthe FETI method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] (for Finite Element Tearing and Interconnecting) -and we take into account the flexibility of such a method to extend the original algorithm to incompressible or nearly incompressible linear elasticity problems, while keeping the coupling between displacement and pressure inside each subdomain. We first discuss the treatment of incompressibility with a mixed formulation, for the cases where the pressure discretization is discontinuous, and propose a first extension of the method that preserve the original performances obtained on compressible cases. The test cases concern 2D plane strain problems and illustrates the influence of the size of the problem, and of the number of subdomains. The resulting algorithms are still optimal and numerically scalable. Finally, an augmentation of the method dedicated to (nearly) incompressible problems is proposed to improve the previous performances. Both the original FETI and the dual-primal FETI-DP [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part I: a faster alternative to the two-level FETI method[END_REF] methods are extended to (nearly) incompressible problems and lead to the FETI-I and FETI-DPI approaches.

Treatment of incompressibility

For compressible linear elasticity, the material behavior is modelled with the Hooke's law that relates the stress σ to the strain ε through a bijective linear operator D. With an isotropic material, only 2 independent coefficients are needed, for instance, the Young's modulus E and the Poisson's ratio ν, or the Lamé's coefficients µ (shear modulus) and λ:

σ = Dε = 2µε + λ(Tr ε)1 = E 1 + ν ε + νE (1 + ν)(1 -2ν) (Tr ε)1
For nearly incompressible materials, when U is the displacement field, div U = Tr ε(U ) → 0, ν → 0.5 and λ → ∞, such that λ(Tr ε) remains finite. Without a particular care, this coefficient is expected to be poorly evaluated [START_REF] Ladevèze | Accuracy and optimal meshes in finite element computation for nearly incompressible materials[END_REF]; moreover, when a purely displacement-oriented finite element formulation is used, the condition number of the rigidity matrix increases as ν goes close to 0.5. This phenomenon is amplified when the size of the finite element model increases, and can prevent a classical iterative algorithm to converge, if applied to such an ill-conditioned problem.

Mixed formulation

To circumvent the previous difficulties, a mixed formulation is classically used, see [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF][START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Pironneau | Finite Element Methods for Fluids[END_REF][START_REF] Herrmann | Elasticity equations for incompressible and nearly incompressible materials by a variational theorem[END_REF] for instance. Both the displacement and the pressure (related to the term H = λ Tr ε) are discretized. The drawbacks are a larger problem with an higher fill-in, and a rigidity matrix which is no more symmetric positive definite.

In order to simplify, consider a structure Ω made with a linear isotropic material. It is subjected to a body force f d and an external force field F d on a part ∂ 2 Ω of its boundary. A classical variational formulation (also called perturbated Lagrangian) consists in finding the saddle point (U , H) of:

Ω µ Tr[ε(U )ε(U )] - 1 2λ H 2 + H Tr ε(U ) dΩ - Ω U • f d dΩ - ∂2Ω U • F d dS
Once U and H are correctly discretized, the corresponding linear system of equations is:

K C T C -N u h = f 0 (1)
u and h are the vectors of unknowns describing the fields U and H; K is a standard rigidity matrix associated to a material with a shear modulus µ and a null second Lamé coefficient; C is a coupling term and N becomes null when the physical material is incompressible:

u T C T h = h T Cu = Ω H Tr ε(U )dΩ ; h T 1 N h 2 = Ω H 1 1 λ H 2 dΩ
K and N are both symmetric definite positive (SPD), but this is not the case for the matrix of the system (1). For compressible cases, N is regular and the elimination of the unknown h in (1) can be performed with a condensation and leads to: (K + C T N -1 C)u = f . When ν → 0.5, N → 0 and this formulation is equivalent to a penalization procedure, intended to prescribe Tr ε(U ) ≈ 0.

For incompressible cases, N = 0 and h is interpreted as a Lagrange multiplier in (1) that enforces the discrete form of the divergence-free condition:

Cu = 0.
A trivial generalization of the previous formulation can be obtained with a slightly different splitting of spherical terms in the constitutive relation:

σ = 2µε + λ(Tr ε)1 + p1
with p = (λ -λ) Tr ε; λ is a parameter to be chosen ( λ = 0 recovers the previous formulation with p = H). The corresponding variational formulation is to find the saddle point of:

Ω 1 2 Tr[ε(U ){2µε(U ) + λ Tr ε(U )1}] + p Tr ε(U ) - 1 2(λ -λ) p 2 dΩ+ - Ω U • f d dΩ - ∂2Ω U • F d dS
and the discrete linear system of equations is:

K C T C -M u p = f 0 (2)
K is a standard rigidity matrix associated to a material with a shear modulus µ and a second Lamé coefficient equal to λ; M still becomes null when the physical material is incompressible:

p T 1 M p 2 = Ω p 1 1 λ - λ p 2 dΩ

Field discretization

To avoid spurious terms in the solution and locking phenomena, the discretization of both fields U and p must be chosen accordingly to the LBB condition, or must pass successfully the patch tests, [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Zienkiewicz | The patch test for mixed formulations[END_REF][START_REF] Dvorkin | On the convergence of incompressible finite element formulations[END_REF]. In this paper, we use a discontinuous discretization of the pressure p, a continuous one for the displacement U , and the finite element in all of the forecoming bidimensional examples is the straight-edge triangle with a linear interpolation of the pressure, a quadratic interpolation of the displacement augmented with a third order bubble (P2 bubble -P1 discontinuous). This element is LBB stable and 2d order.

Domain decomposition

In order to solve the systems (1) or (2), several strategies can be used. A direct monolithic factorisation requires a suitable degree of freedom (dof) renumbering scheme and/or a suitable dof pivoting due to the lack of SPD property for the global problem matrix. Such a direct approach suffers from a high computational cost when the size of the problem increases. An iterative scheme can then be of interest to keep these costs affordable. Among them, partitioning techniques can be applied: for instance, a physical splitting between displacement and pressure dof, or a geometric splitting of the domain in subdomains. We are interested herein in such methods: the domain is partitioned into non overlapping subdomains, the displacement-pressure coupled problem remains coupled on each of them. If the pressure is not subjected to continuity conditions, "gluing" requirements on the interface between subdomains is only concerned with a null displacement jump.

The algorithms we are developing are based on the dual multilevel domain decomposition FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. More precisely, two versions of this method, the socalled FETI and FETI-DP formulations, are extended to get the FETI-I and The proposed test problems used throughout this paper concern an idealized 2D incompressible polymer flow through an extrusion die, treated with plane strain assumption (remark that with plane stress, the incompressibility condition does not prescribe any constraint on the in-plane kinematics). Due to symmetry, only half of the problem is modelled on Figure 1. A parabolic displacement field is prescribed on the left, while there is a free surface on the right. The structure is decomposed into n subdomains Ω s , s = 1, . . . , n, and a global interface Γ between them. h denotes a characteristic size of the finite elements, and H, the size of the subdomains. With an overall length of the problem set to 1, the subdomain size H varies from 1/6 to 1/15 when the number of subdomains increases from 12 to 75, on Figure 2. Figure 3 presents the solution obtained with n = 27 subdomains and a mesh refinement h/H = 1/6.

H h Ω s Γ Ω U.y = 0 U = 0 F = 0 U = U d 1

Extension of the FETI method to (nearly) incompressible problems

The FETI method is a dual domain decomposition method: the discrete kinematic "gluing" condition between subdomains is dualized via Lagrange multipliers [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. When using a conforming finite element mesh for all the domain Ω, the Lagrange multipliers λ are defined at each node of the global interface Γ between subdomains. They represent the corresponding nodal force acting from a subdomain s on a subdomain s ′ (when using the convention s < s ′ to uniquely define the sign of λ). Therefore, each subdomain s needs a signed boolean matrix B (s) to map the global force vector λ onto the forces acting on its boundary: -B (s) T λ. The kinematic constraint on the jump of displacement between neighboring subdomains on the global interface is therefore:

✂✁ ✄ ✂☎ p ✆ ✝ ✟✞ ✁ ✂✠ ✝ ✂✁ ✂✡ ✝ ✂✁ ✡ ✂✄
s B (s) u (s) = 0 (3)
The finite element equilibrium of a subdomain s, submitted to the external generalized forces f (s) and the action of the multipliers, is:

K (s) H u (s) p (s) = f (s) 0 -B (s) T 0 λ with K (s) H = K (s) C (s) T C (s) -M (s) (4) K (s)
H is symmetric, but not positive and not obviously definite: for floating subdomains, i.e. subdomains with no prescribed displacement, or with internal mechanisms, the kernel of K (s) H contains rigid body motions. The set of such rigid body motions and null pressure of such a subdomain s is denoted by R (s) .

The first step to build the FETI-I method is to condense the unknowns u (s) and p (s) onto the multiplier λ: provided that the right hand side of ( 4) is orthogonal to R (s) , i.e.

R (s) T f (s) 0 -B (s) T 0 λ = 0 ( 5
)
the solution is obtained with an indefinite rigid body motion field using one of the pseudo-inverses

K (s) H + of K (s)
H : (4) is equivalent to

u (s) p (s) = K (s) H + f (s) 0 -B (s) T 0 λ + R (s) α (s)
α (s) are the arbitrary coordinates of the rigid body motion basis R (s) . Using this expression in the constraint (3) leads to

F λ -Gα -d = 0 with F = s B (s) 0 K (s) H + B (s) T 0 d = s B (s) 0 K (s) H + f (s) 0 G = B (1) 0 R (1) . . . B (n) 0 R (n) α T = α (1) T . . . α (n) T
In order to complete the problem, the equations ( 5) are needed for each subdomain. Assembling them leads to:

G T λ = e with e = R (1) T f (1) 0 . . . R (n) T f (n) 0 (6)
Finally, the dual problem is to find (λ, α) such that:

F -G -G T 0 λ α = d -e (7)
Table 1: The projected conjugate gradient algorithm

Initialize λ 0 = G(G T G) -1 e r 0 = d -F λ 0 Iterate k = 1, 2, • • • until convergence Project w k-1 = P T r k-1 Precondition z k-1 = F -1 w k-1
Re-project

y k-1 = P z k-1 Conjugate ζ k = y k-1 T w k-1 y k-2 T w k-2 (ζ 1 = 0) Search direction p k = y k-1 + ζ k p k-1 (p 1 = y 0 ) Matrix-vector product q k = F p k Line search η k = p k-1 T w k-1 p k T q k Update λ k = λ k-1 + η k p k Residual r k = r k-1 -η k q k
A key point is that, though we use a mixed formulation, F is still symmetric positive (not definite, especially when there are redundancies in kinematic constraints (3)).

To prove it, consider first the initial formulation (1), and an arbitrary multiplier field λ; for the subdomain s, define

v (s) π (s) = K (s) H + B (s) T λ 0 Then, λ T F (s) λ = B (s) λ T 0 K (s) H + B (s) T λ 0 = B (s) λ T 0 v (s) π (s) = B (s) λ T 0 v (s) -π (s) = v (s) T π (s) T K (s) H v (s) -π (s) = v (s) T K (s) v (s) + π (s) T N (s) π (s)
As K (s) and N (s) are positive, F (s) and F = s F (s) are positive too. The next step is to recall that the generalized formulation ( 2) is strictly equivalent to the previous one for the local solution u (s) , and p (s) = λh (s) /(λ -λ), so the positivity is ensured for this formulation also, as soon as λ = λ. Therefore, a projected conjugate gradient is applicable to solve [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]. A projection is required to enforce the admissibility condition (6) for λ. The projector is for instance

P = 1 -G(G T G) -1 G T .
The corresponding overall algorithm does not differ from the classical FETI method; in order to be self-contained, it is recalled in Table 1.

One has to notice that

• at each iteration, the algorithm requires 2 projections; • each projection is a coarse problem on the entire domain with the rigid body motion coefficients α as unknowns;

• the computational kernel is the matrix-vector product F v. F is never explicitly assembled; the product is performed in parallel as it requires independent Neumann problems on each subdomains;

• the preconditioner F -1 is intended to be easily parallelizable. Its design is recalled in a following section.

This approach is a straightforward extension of the FETI method to mixed problems with a discontinuous discretized pressure field. Using a continuous approximation of the pressure would require a "gluing" in the pressure field at the interface. The operator F would be no more positive, and so, an other iteration scheme is required (GMRes for instance); this feature is out of the scope of the present paper.

An other situation is very similar to the one we consider herein: the case where the pressure is chosen to be continuous over any subdomain, but discontinuous throughout the interface as discussed in [START_REF] Calgaro | On the domain decomposition method for the generalized Stokes problem with continuous pressure[END_REF]. In such a case, the discretized solution obviously depends on the domain decomposition.

Extension of the dual-primal version to (nearly) incompressible problems

Recently, a modified version of the FETI method have been designed under the FETI-DP acronym (for FETI Dual-Primal) [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part I: a faster alternative to the two-level FETI method[END_REF][START_REF] Farhat | A scalable dual-primal domain decomposition method[END_REF]. Using an upgraded gluing kinematical condition, it allows the local rigidities of the subdomains to be regular, and needs only one coarse problem resolution per iteration. The key point is the definition of "corner" nodes that will constitute the coarse nodes of the problem. A simple way to define them is to use the geometric corners of the subdomains in 2D, but other possibilities exist, as soon as clamping these nodes leads to a regular Neumann problem on every subdomain, and to a regular coarse problem [START_REF] Lesoinne | A FETI-DP corner selection algorithm for three-dimensional problems[END_REF].

For the example of the Figure 1, and n = 27 subdomains, the coarse nodes are represented on Figure 4.

For (nearly) incompressible problems, it begins also with the splitting of the local dof (per subdomain) into local corner dof u

∀ u (s) r ⋆ , p (s) ⋆ , u (s) c ⋆ = B (s) c u ⋆ c ,    u (s) c ⋆ u (s) r ⋆ p (s) ⋆    T K (s) H    u (s) c u (s) r p (s)    =    u (s) c ⋆ u (s) r ⋆ p (s) ⋆    T    f (s) c f (s) r 0    -    u (s) c ⋆ u (s) r ⋆ p (s) ⋆    T    0 B (s) r T 0    λ r with K (s) H =     K (s) cc K (s) cr C (s) c T K (s) rc K (s) rr C (s) r T C (s) c C (s) r -M (s)    
It can be split in the local part per subdomain s:

K (s) Hrr u (s) r p (s) = f (s) r 0 -B (s) r T 0 λ r - K (s) rc C (s) c B (s) c u c with K (s) Hrr = K (s) rr C (s) r T C (s) r -M (s) (8)
and the remaining coarse part which is assembled on the entire domain:

s B (s) c T K (s) cc B (s) c u c + s B (s) c T K (s) cr C (s) c T u (s) r p (s) = s B (s) c T f (s) c With a correct choice of local corner dof, K (s)
Hrr is regular and ( 8) is equivalent to:

u (s) r p (s) = K (s) Hrr -1 f (s) r 0 -B (s) r T 0 λ r - K (s) rc C (s) c B (s) c u c (9)
The kinematic constraints on the interface depends only on u (s)

r dof: s B (s) r 0 u (s) r p (s) = 0
With the previous expression of the local solution, (3.2) and ( 9) lead to:

F rr F rc F cr -K ⋆ cc λ r u c = d r -f ⋆ c (10)
with

F rr = s B (s) r 0 K (s) Hrr -1 B (s) r T 0 d r = s B (s) r 0 K (s) Hrr -1 f (s) r 0 F rc = s B (s) r 0 K (s) Hrr -1 K (s) rc C (s) c B (s) c F cr = F T rc K ⋆ cc = s B (s) c T K (s) cc ⋆ B (s) c K (s) cc ⋆ = K (s) cc -K (s) cr C (s) c T K (s) Hrr -1 K (s) rc C (s) c f ⋆ c = s B (s) c T f (s) c ⋆ B (s) c f (s) c ⋆ = f (s) c -K (s) cr C (s) c T K (s) Hrr -1 f (s) r 0
The last step is the condensation of the corner dof u c on the Lagrange multipliers λ r in [START_REF] Mandel | Balancing domain decomposition[END_REF],

u c = K ⋆ cc -1 (f ⋆ c + F cr λ r )
, to get:

F ⋆ λ r = d ⋆ with F ⋆ = F rr + F rc K ⋆ cc -1 F cr and d ⋆ = d r -F rc K ⋆ cc -1 f ⋆ c
(11) Notice that the matrix problem F ⋆ is still symmetric positive. The proof that F rr is positive is the same as the one previously used for FETI-I algorithm. To prove that

K ⋆ cc is positive, let us consider    v (s) c v (s) r π (s)    = K (s) H +   f (s) c 0 0   ⇒ K (s) H    v (s) c v (s) r π (s)    =   f (s) c 0 0   ⇒ K (s) cc ⋆ v (s) c = f (s) c
As it has already been proved that c must be the restriction to local c dof of a rigid body motion u (s) ; and (ii) progressing from a subdomain to its neighbors, if sufficiently connected by some local c dof, one can conclude that all u (s) are the restriction on the subdomains of the same global rigid body motion u. If the initial problem is regular, the only admissible rigid body motion is u = 0, so u c = 0 too. This concludes to the regularity of K ⋆ cc . There are therefore less severe conditions on the choice of corner modes to ensure regularity of K ⋆ cc in such a way that less corner nodes are required. They are discussed in [START_REF] Lesoinne | A FETI-DP corner selection algorithm for three-dimensional problems[END_REF] and are directly applicable to the present case. Notice that, in 2D, when each subdomain shares with every neighbor that has a common edge at least two corner nodes, the regularity is guaranteed, though the number of corner nodes is not minimal. This is the case for the following examples in this paper.

f (s) c T 0 0 K (s) H +   f (s) c 0 0   ≥ 0 one gets f (s) c T 0 0    v (s) c v (s) r π (s)    = f (s) c T v (s) c = v (s) c T K (s) cc ⋆ v (s) c ≥ 0 So, K (
To discuss the definiteness of F ⋆ , let us now consider

λ T r F ⋆ λ r = λ T r F rr λ r + (λ T r F rc )K ⋆ cc -1 (F cr λ r ) = 0
Due to the positivity of F rr and K ⋆ cc , this leads to λ T r F rr λ r = 0. With

v (s) r π (s) = K (s) Hrr -1 B (s) r T λ r 0 10 as for FETI-I, it leads to s v (s) r T K (s) rr v (s) r +π (s) T N (s) π (s) = 0. With K (s)
rr and N (s) SPD, we can conclude to v (s) r = 0 and π (s) = 0. So, B (s) r T λ r = 0 for any subdomain s. If there is no redundant multiplier (i.e. if each multiplier enforces an independant gluing condition), this allows us to conclude to λ r = 0. So, F ⋆ is definite, and therefore, is SPD. Notice that with redundant multipliers, F ⋆ is not definite any more, but this does not prevent a conjugate gradient algorithm to converge.

The FETI-DPI algorithm is a classical conjugate gradient applied to the problem [START_REF] Cowsar | Balancing domain decomposition for mixed finite elements[END_REF]. The algorithm is then similar with the one in Table 1, but no more projection is needed.

One has to notice that

• at each iteration, this algorithm requires the solution of only one global coarse problem consisting in finding the displacement of corner dof, during matrix-vector product F ⋆ v;

• as previously, the choice of the preconditioner is discussed in the next section.

A first choice of preconditioner and the corresponding comparisons

The first idea for preconditioning is to re-use the classical preconditioners in the case of a compressible problem. The so-called Dirichlet preconditioner is fully parallel and solves a Dirichlet-like problem per subdomain (with a prescribed displacement on the boundary):

F -1 D = s W B (s) S (s) bb 0 0 0 B (s) T W (12)
F -1 D applies to a kinematics jump on the interface and returns a generalized force on the same interface. W is a diagonal matrix storing in each of its entries the inverse of the multiplicity of an interfaced dof, and the subscript b denotes the dof that lie on the boundary of the subdomain, and the subscript i, the remaining ones. In particular, the standard rigidity matrix can be split in:

K (s) = K (s) bb K (s) bi K (s) ib K (s) ii
and a compressible preconditioner is:

S (s) bb = K (s) bb -K (s) bi K (s) ii -1 K (s)
ib . This preconditioner is a function of the coefficient λ, or equivalently of the coefficient ν = 1 2 λ/( λ + µ), used to compute the stiffness K (s) . This preconditioner can be applied to each of the FETI-I and FETI-DPI methods (for the second one, the matrix B (s) r is used and b dof are a partition of the r dof; there is no residual belonging to c nodes because there is no kinematic jump on these nodes all along iterations).

To compare the proposed methods to the standard FETI and FETI-DP ones, the test case proposed in Figure 1 is used. The number of subdomains is n = 48 (which corresponds to H = 1/12) and a refinement of the mesh such that h/H = 1/6. The different algorithms have been implemented with Matlab TM . 2 report the number of iterations required to reach convergence. The convergence criteria is w k / d ≤ 10 -7 . This number of iterations is plotted versus the Poisson's ratio of the real material.

For compressible cases (ν < 0.5) both classical methods and adapted methods (with a mixed formulation) have similar performances; the methods suited to incompressible problems are then not competitive with the classical ones for this regime, as their finite element problems are larger, due to the pressure unknowns. For the mixed formulations, the results are plotted for different values of the parameter ν. When ν approaches 0.5, the bad conditioning of the preconditioner decreases global performance, and the number of iterations is higher than for the other cases. Taking ν = 0.29 or even ν = 0 (classical mixed formulation) leads to good results.

For nearly incompressible cases (ν ≈ 0.5), the classical methods exhibit difficulties to converge because the problem to solve is equivalent to a penalized one, and so, is ill-conditioned. For FETI-I and FETI-DPI methods, the convergence is preserved, even for the purely incompressible case (ν = 0.5).

The characteristics of the original FETI methods, for compressible linear elasticity problems, concern the condition number κ of the iterations on the interface problem: it is asymptotically bounded as κ = O(1 + log H h ) m . Usually m = 3, or m = 2 with an appropriate projector, see [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: connections and new results[END_REF].

Optimality test

The optimality test consists in keeping the number n of subdomains constant while increasing the size of the problem, i.e. decreasing the ratio h/H. The algorithm is optimal if the number of iteration weakly depends on h/H. Figure 6 and Table 3 report this number of iterations versus the ratio h/H. The different curves correspond to an incompressible problem (ν = 0.5) and to ν = 0.3 for FETI-I and FETI-DPI methods. In order to compare to a reference, the number of iterations for a compressible case (ν = 0.3) and standard FETI is also provided. FETI-IA and FETI-DPIA entries correspond to augmented versions of the algorithms that will be discussed in a later section.

When using n = 75 subdomains, and h/H = 1/10, the problem possesses 15 000 elements, 45 000 pressure dof and about 61 000 displacement dof (without the bubble dof which has been condensed at the element level).

All the tested approaches are optimal; nevertheless, for the incompressible case, the convergence is not as efficient as for the compressible case: the number of iterations is about 35% higher. Dual-primal versions of the methods need systematically less iterations, but their costs are not necessarily the same. The section 5 will discuss this last point.

Numerical scalability test

The numerical scalability test consists in keeping a local problem of constant size (h/H constant) while increasing the number of subdomains (decreasing H), and so, increasing the size of the global problem. Figure 7 and Table 4 report the results for the same kind of situations than for the previous test. The scalability is similar for all the methods. The previous remarks concerning the number of iterations are still applicable.

A second Dirichlet-like preconditioner

The next idea for a preconditioner is to use the same mixed formulation as for the direct problem. The form of the preconditioner is still the same, but now:

S (s) bb = K (s) bb -K (s) bi C (s) b K (s) Hii -1 K (s) ib C (s) b T with K (s) Hii = K (s) ii C (s) i T C (s) i -M (s) (13) S (s)
bb does not depend any more on ν so, for compressible cases, the results are expected to be similar to those of the previous preconditioner.

For incompressible case, this preconditioner cannot be applied since

K (s)
Hii is singular; the next section will discuss this topic.

For the (nearly) incompressible case, using this preconditioner is not obviously efficient: the residual (kinematic jump on interface) can have components with a significant volume variation per subdomain. Due to its nearly incompressibility nature, the preconditioner will return a force field on the boundary with a high hydrostatic component, and so, will impair convergence. To exemplify this point, let us consider again the test case of the die problem with a nearly incompressible material (ν = 0.4995), n = 48 subdomains and h/H = 1/6. We use FETI-I method with ν = 0.3 and different preconditioners: F -1 C for compressible preconditioner of section 3. 3, andF -1 H for this last preconditioner using [START_REF] Pavarino | Balancing Neumann-Neumann methods for incompressible Stokes equations[END_REF]. The convergence rate is related to the condition number κ of the iteration matrix F -1 D P T F P . Once we get rid of null eigenvalues corresponding to the unbalanced eigenvectors, and to the redundant gluing conditions, we obtain: κ(P T F P ) = 586, κ( F -1 C P T F P ) = 255 and κ( F -1 H P T F P ) = 1114. As the condition number is higher for the nearly incompressible preconditioner, the convergence is expected to be lower. This last preconditioner is therefore not useful if used alone. Nevertheless, its design is a motivation for the augmentation of the algorithms, in the next section.

Augmentation of the algorithms

Up to now, the basis of the FETI-I and FETI-DPI algorithms are very similar to the original FETI and FETI-DP ones because they do not take into account the particularities of the problem to be solved. This feature could improve the performance of the algorithms. For instance, for an incompressible case, one may use an incompressible preconditioner. This can be used only if the residual to which this preconditioner is applied, satisfies the constant volume condition per subdomain. Such a constraint is an augmentation.

This technique has been previously used to accelerate the convergence by prescribing additional constraints to the successive iterates of the algorithm [START_REF] Farhat | A unified framework for accelerating the convergence of iterative substructuring methods with Lagrange multipliers[END_REF][START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF][START_REF] Farhat | A scalable dual-primal domain decomposition method[END_REF]. The drawback is the increase of the coarse problem size. If the acceleration counterbalances this overhead of cost, the performance of the overall algorithm is improved. It is now known in the litterature that an augmentation is mandatory for FETI-DP to be scalable for second order 3D problems.

We propose in this paper an augmentation suited to the kind of problems we are dealing with. As a first step, this augmentation is built for the incompressible case.

Augmentation of FETI-I algorithm for incompressible cases

In this situation, if one wishes to use an incompressible preconditioner, the residual it is applied to must satisfy a solvability condition: the volume of the subdomain s which is completely surrounded by neighbors (or prescribed displacement boundary) must be constant when a displacement V is prescribed by the preconditioner on its boundary. Let us denote by U its prolongation inside the subdomain Ω s (U |∂Ωs = V ); the condition is:

Ωs Tr ε(U )dΩ = ∂Ωs V • ndS = 0
where n is the outward unit normal vector on the boundary. Once discretized, this condition is equivalent to p (s) 1 T C (s) u (s) = 0 with a uniform pressure field p (s) 1 on the whole subdomain (for instance equal to unity). The constraint matrix is then easily computed from the coupling matrix C (s) : c (s) = p (s) 1 T C (s) . The Dirichlet preconditioner [START_REF] Tallec | Non-overlapping domain decomposition methods for adaptive hp approximations of the Stokes problem with discontinuous pressure fields[END_REF] and ( 13) acts on the projected residual: s) . For the subdomain s, B (s) T W P T r are the displacements computed in the beginning of the preconditioning step, so the constraint must be enforced on them: (1) . . . Q (n) , which has a full column rank, the n additional constraints are

P T r = P T s B (s) u (
Q (s) T P T s B (s) u (s) = 0 with Q (s) T = c (s) B (s) T W . With Q = Q
Q T P T s B (s) u (s) = 0 ( 14
)
As this constraint is obviously satisfied for the solution of the problem (for which P T s B (s) u (s) = 0), it can be prescribed to any subdomain, even if it possesses a free boundary. Doing this avoids the detection of the subdomains that have a singular Neumann problem during preconditioning.

To enforce these constraints, an additional multiplier γ enriches the previous ones [START_REF] Farhat | A unified framework for accelerating the convergence of iterative substructuring methods with Lagrange multipliers[END_REF][START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF], and satisfies automatically the admissibility condition: λ + P Qγ replaces λ in the previous FETI-I algorithm (4), (3), [START_REF] Stenberg | Mixed hp finite element methods for problems in elasticity and Stokes flow[END_REF].

γ can be interpreted as a uniform pressure value per subdomain. Keeping it in a coarse problem allows the solution to satisfy the additional constraints at each iteration, and so, is expected to improve convergence.

With the equations ( 14), the overall problem to be solved is:

  F F P Q -G Q T P T F L 0 -G T 0 0     λ γ α   =   d Q T P T d -e   with L = Q T P T F P Q
To build the corresponding algorithm, the following ways are equivalent: using the FETI2 framework [START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF] with the interpretation of an additional coarse space correction with the regular matrix L, or, as it is done herein, condensing the new unknowns γ on the other ones to get:

F ⋆ -G -G T 0 λ α = d ⋆ -e with F ⋆ = F -F P QL -1 Q T P T F and d ⋆ = d -F P QL -1 Q T P T d
, and apply the previous FETI-I algorithm to this new problem. the overall algorithm is the same as in Table 1; the difference in the implementation is the need for the additional coarse problem (with matrix L) within the matrix-vector product F ⋆ p. The preconditioning step allows to apply the eventually singular (because incompressible) Dirichlet preconditioner to the residual: if s is a subdomain with all of its boundary subjected to a prescribed displacement v (s) = B (s) T W P T s B (s) u (s) , the Dirichlet problem to be solved on this subdomain is:

K (s) Hii v (s) i p (s) = - K (s) ib C (s) b T v (s) b If K (s)
Hii is singular, its kernel is exactly a uniform pressure 0 p (s) 1 T T . As the right hand side is orthogonal to this kernel as soon as the augmentation condition is satisfied, the solution is then defined up to an undetermined uniform pressure field per subdomain. This field is not required to be determined: first, only forces on the boundary are extracted from the preconditioned solution and, second, this correction will be automatically added in the augmented coarse problem within the multiplier γ.

Augmentation of FETI-DPI algorithm for incompressible cases

The same procedure can be applied to the FETI-DPI algorithm. Briefly, the additional constraint is this time:

Q T s B (s) r u (s) r = 0 with Q (s) T = c (s) r B (s) r T W = p (s) 1 T C (s) r B (s) r T
W because at the current iteration, the correction to the solution is performed with no residual on coarse nodes and so, the additional constraint is only dealing with r dof. The corresponding problem to be solved is:

  F rr F rr Q F rc Q T F rr Q T F rr Q Q T F rc F cr F cr Q -K ⋆ cc     λ r γ u c   =   d r Q T d r -f ⋆ c  
with a condensation of the coarse dof γ T u T c T on λ r :

F ⋆ F rr -F rr Q F rc K⋆ cc -1 Q T F rr F cr λ r = d r -F rr Q F rc K⋆ cc -1 Q T d r -f ⋆ c with K⋆ cc = Q T F rr Q Q T F rc F cr Q -K ⋆ cc
The augmented coarse space is then a mixed coarse problem whose unknowns are both the corner dof and the uniform pressure values.

To check if this coarse problem is LBB stable, in particular to avoid spurious pressure oscillations, let us consider the condensation of K ⋆ cc on coarse pressure dof γ:

Q T F rr Q+Q T F rc K⋆ cc -1 F cr Q = Q T F ⋆ Q. If F ⋆ is SPD (especially when
there is no redundant multiplier), this condensed matrix is symmetric positive. Moreover, if we consider γ T Q T F ⋆ Qγ = 0, we obtain Qγ = 0. As Q is full column rank, one gets γ = 0, so the condensed matrix is SPD too. Its minimum eigenvalue is then bounded away from zero, and the coarse matrix problem is stable.

General case and numerical results

As already mentioned, the augmentation consists in prescribing a constraint that will be satisfied at convergence, because it enforces the residual to be orthogonal to a wider subspace: Q T s B (s) u (s) = 0. The present augmentation can obviously been applied to nearly incompressible cases as well, without any modification of the algorithm. Now, the different possibilities are the following, when the material is nearly incompressible or purely incompressible: • for the non-augmented algorithms FETI-I and FETI-DPI: the preconditioner can be compressible or nearly incompressible;

• for the augmented algorithms denoted with FETI-IA and FETI-DPIA: the preconditioner can be compressible, nearly incompressible or incompressible.

The compressible Dirichlet preconditioner has already been tested for FETI-I and FETI-DPI algorithms. To compare with FETI-IA and FETI-DPIA, the results have already been given in Tables 3 and4. Obviously, the additional constraints improve convergence. The results are very similar to the compressible case, in terms of number of iterations.

As already mentioned, the treatment of compressible problems with the proposed algorithms is not efficient. For nearly incompressible problems and incompressible problems, the performance in terms of the number of iterations with different preconditioners can be found in Table 5. The tested problem is the die problem with n = 48 subdomains, h/H = 1/6 and for ν = 0.3. As it is not clearly interesting to use an incompressible preconditioner for nearly incompressible problems or nearly incompressible preconditioners for incompressible problems, the corresponding results are not reported.

Clearly, in each case, the convergence rate in increased when switching from FETI to FETI-DP, and with the augmentation of the algorithms. Concerning the influence of the preconditioner, the situation is not so clear: the compressible preconditioner (due to its regularity) is more efficient than other ones. Even with purely incompressible problems with augmentation, the use of compressible preconditioner is to be recommended. So, this is now the choice we made, in particular in the next section. of efficiency. It is not able to measure the effective parallel potentialities of an algorithm, nor its ability to take into account a particular architecture of a parallel machine to tune the algorithm (as for unlooping, vectorizing or locality of datas) and the cost of data access in memory. On the positive side, the floating point operation count does not interfere with the performance of a particular operating system and compiler; it is close to the algorithmic part and the main choices done when implementing it (full entry matrix, skyline solver. . . ). It provides also information on how the cost scales with respect to the parameters of the problem solved. For instance, we choose to consider problems with a regular decomposition into subdomains (m × m subdomains in 2D, and m × m × m subdomains in 3D); and regular meshes for each subdomain. This allows to determine the structure of the sparse matrices involved, independently of the dof renumbering scheme used.

We considered herein banded symmetric matrices, and we assumed that the renumbering scheme, the solver for local subdomain non SPD problems, and the eventual pivoting scheme, see [START_REF] Duff | The solution of augmented systems[END_REF], does not destroy the bandwidth. This assumption is optimistic for large number of dof, as the probability to detect a null pivot that does not correspond to a rigid body motion, increases. In particular, when needed, the detection of rigid body motions in subdomains will be preferably performed with a geometric inspection [START_REF] Farhat | On the general solution by a direct method of a large-scale singular system of linear equations: application to the analysis of floating structures[END_REF][START_REF] Papadrakakis | An integrated geometric-algebraic method for solving semi-definite problems in structural mechanics[END_REF].

Matrix structures

In all of the following cost estimations, we separated the initialization phase of the various algorithms (computation of elementary matrices; local factorizations; computation of coarse matrix; corresponding right hand side, and its factorization; factorization of the preconditioner) and the cost of an iteration (local solves for matrix-vector product of conjugate gradient, and for preconditioner; coarse solves; global operation on interface fields).

The assembling of matrices and vectors, or disassembling of global vectors (at element or subdomain levels) have been neglected as soon as they involve boolean matrices B (s) , B (s) c , B (s) r . . . For FETI-I, the coarse space matrix G T G is symmetric and possesses the connectivity of neighboring subdomains. For FETI-DPI, K ⋆ cc has the same connectivity as a mesh with corner nodes, and subdomains only connected with these nodes as superelements. Concerning augmented algorithms, the structure of the coarse problem is more complex because superlements are also connected by an 'internal' dof corresponding to the coarse pressure γ. Nevertheless, it still possesses a sparse pattern [START_REF] Pierson | A family of domain decomposition methods for the massively parallel solution of computational mechanics problems[END_REF], and with the same assumption as for the subdomain stiffness factorization, the costs are estimated with the corresponding banded matrix storage, and costs are estimated accordingly. For FETI-IA, the second level coarse problem matrix L = Q T P T F P Q is full and symmetric, due to the long distance coupling of the projection operator. For the tested structure of the die problem (with h/H = 1/6 and n = 48 subdomains) the corresponding morse pattern is also plotted on the previous expression.

Implementation choices

Concerning implementation efficiency, the idea is to store the information (rather than recompute it on-the-fly several times when iterating) as soon as it concerns only fields on the boundary that need only to be stored distributedly at the subdomain level.

For instance, for FETI-DPI, during initialization phase, the vectors

V (s) = K (s) Hrr -1 K (s) cr C (s) c
are precomputed with the resolution of as many local systems as there are corner nodes touching the subdomain s (and stored only on the boundary dof b for iteration phase). Then, K

s) cc ⋆ = K (s) cc -K (s) cr C (s) c T V (s) , f (s) c ⋆ = f (s) c -V (s) T f (s) (
r and for matrix vector product

F ⋆ λ r = F rr λ r +F rc K ⋆ cc -1 F cr λ r , with: F s rr λ r = B (s) r 0 K (s) Hrr -1 λ s r 0 with λ s r = B (s) r T λ r F s cr λ r = B (s) c T V (s) T λ s r u c = K ⋆ cc -1 F rc λ r = K ⋆ cc -1 s F s rc λ r F rc u c = s B (s) r V (s) u (s) c with u (s) c = B (s) c u c
For FETI-DPIA, the entries of F cr Q are computed as follows: the column j is

s B (s) c T V (s) T B (s) r T Q (j) 0
and is computed only with the contributions of the subdomains s neighbors of the subdomain j (due to the term B j) , the other contributions are null). The term on line i and column j of Q T F rr Q is therefore:

(s) r T B j r in B (s) r T Q (
s Q (i) T B (s) r 0 A sj with A sj = K (s) Hrr -1 B (s) r T Q (j) 0
Similarly, A sj is stored on a subdomain s database only for neighbors j, and only for local boundary dof b. Then

F ⋆ λ r = F rr λ r -F rr Q F rc K⋆ cc -1 Q T F rr F cr λ r with F rr Qµ = s B (s) r 0 j A sj µ j
and similar computations of F rr λ r and F cr λ r .

Numerical Results

With the previous scenario, the cost estimation can be plotted with respect to the parameters H and h/H for each of the proposed algorithms. In 2D on Figure 8, the initialization phase has been separated from the iterations, and in order to have similar costs, 20 iterations are considered. FETI-I and FETI-IA are plotted with dashed lines, while DP versions are plotted with plain lines. To give an idea of the behavior of the methods, a 3D cost estimation has been done as well and the results are plotted on Figure 9 for the initialization phase, and for 100 iterations. Both for the initialization and iteration phases, we can consider as a first step that the costs of the previous algorithms arise from different sources:

• the coarse problem: the cost is related first to the number of subdomains, and second, to the structure of the coarse problem itself;

• the subdomain problems: the cost is mainly related to the size of the local problems on subdomains; it is driven by the coefficient h/H;

• additive operations like dot products on a global interface field, which can be performed at subdomain level with few assembly or disassembly operations.

Usually, the costs of the third category are negligible when compared to the previous ones.

If the amount of operations that are done at the subdomain level is sufficiently large when compared to the cost related to the coarse problem (i.e. if h/H is sufficiently small), both the initialization and iteration costs have similar evolutions when h/H decreases. Due to the higher fill-in in the coarse problem of FETI-IA, the initialization cost grows very rapidly when the number of subdomain increases. Therefore, FETI-IA is probably not the most efficient approach to use, for the cost estimation, as well as for the complexity of implementation.

Concerning iteration phase, the cost is almost the same for all of the approaches, especially for 3D case. Concerning initialization phase, the cost increases from FETI-I, to FETI-DPI, and up to FETI-DPIA. To compare these approaches, let us recall the convergence results for 2D incompressible case, with compressible preconditioner, of Table 5. Table 6 evaluates the relative costs with the previous complexity analysis (n = 48, H = n -0.5 , h/H = 1/6). It appears that for the proposed test case, using dual-primal version is more efficient, as when an augmentation is used. Moreover, as it has been previously highlighted in original papers on FETI-DP [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part I: a faster alternative to the two-level FETI method[END_REF], it consumes less iterations and is more robust than the previous FETI and FETI2 versions.

Summary and conclusions

The framework of standard FETI and FETI-DP versions with or without augmentation has not been modified. To do so, we derived specific implementation for the mixed local problems per subdomain, without changing gluing condition between them, thanks to a discontinuous pressure discretization. This feature let us to expect an easy integration of the extension for incompressible or nearly incompressible material with other developments of the FETI method, especially for dynamics problems [START_REF] Farhat | A scalable Lagrange multiplier based domain decomposition method for implicit time-dependent problems[END_REF] and contact problems [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF].

Without augmentation of the algorithm, optimality and numerical scalability have been obtained in 2D. They are expected to pertain in 3D for the FETI-I method, but not for the FETI-DPI without augmentation, as it has been previously noticed for compressible problems [START_REF] Farhat | A scalable dual-primal domain decomposition method[END_REF]. If the proposed augmentation (only one constraint per subdomain) is not sufficient enough to recover the numerical scalability in 3D, there will be no difficulty to add the previously developed augmentations [START_REF] Farhat | A scalable dual-primal domain decomposition method[END_REF], to the specific incompressible one proposed herein. With this specific augmentation, the rigid body motion for FETI-I (and corner modes for FETI-DPI) as well as a uniform pressure per subdomain is propagated globally at each iteration. The performances obtained are then similar to the one obtained for compressible problems.

Of course, though the use of uniform pressure for the augmentation is required in the case of purely incompressible behavior, other additional stress modes may be used, for instance the generalized Trefftz basis.

The increase in size of the coarse problem leads to improved convergence rate, but increases also the cost of the coarse problem (a complexity analysis is required to predict the gain in performances) and impairs the parallel part of the implementation. With a large number of subdomains, a particular care of the coarse problem implementation is required as it is a bottleneck for the parallelization.

Figure 1 :Figure 2 :

 12 Figure 1: Decomposition in subdomains Ω s , and a global interface Γ

Figure 3 :

 3 Figure 3: Solution to the die problem with n = 27 subdomains and h/H = 1/6; left: displacement field, right: field of variable p

Figure 4 :

 4 Figure 4: Coarse nodes for the die problem with n = 27 subdomains

  corner continuity between subdomains is equivalent to extract u (s) c from a unique global vector of corner dof u c : u (s) c = B (s) c u c . The remaining kinematic continuity conditions on interfaces is still expressed with a Lagrange multiplier λ r and a boolean mapping matrix B (s) r . The variational form corresponding to the equilibrium of the subdomains is:

  s) cc ⋆ is positive, and K ⋆ cc is positive too. To check if K ⋆ cc is regular, let us look for a displacement field such that u T c K ⋆ cc u c = 0. As all K any subdomain s. As all K (s) rr are regular: (i) K (s) cc ⋆ has rigid body motions (restricted to local c dof) as a kernel, so u (s)

Figure 5 :

 5 Figure 5: Compared convergences (number of iterations) versus Poisson's ratio; left: for FETI like methods; right: for FETI-DP like methods

Figure 6 :

 6 Figure 6: Number of iterations depending of the problem size h/H (for n = 75 subdomains)

Figure 7 :

 7 Figure 7: Number of iterations depending of the number of subdomains (for h/H = 1/10)

10 Figure 8 :

 108 Figure 8: Complexity estimation for 2D problems

Figure 9 :

 9 Figure 9: Complexity estimation for 3D problems

Table 2 :

 2 Compared convergences (number of iterations) versus Poisson's ratio

	Poisson's ratio ν	0.3 0.4 0.45 0.49 0.499 0.4995 0.5
	FETI		35	39	46	83	168	183	-
	FETI-I ν = 0	35	37	40	46			50
	FETI-I ν = 0.29	35	36	38	43			47
	FETI-I ν = 0.44	42	43	44	49			52
	FETI-I ν = 0.48	56	57	59	67			72
	FETI-DP		23	25	29	50	77	89	-
	FETI-DPI ν = 0	24	27	28	33			36
	FETI-DPI ν = 0.29 22	24	26	31			34
	FETI-DPI ν = 0.44 28	27	27	33			37
	FETI-DPI ν = 0.48 39	38	37	40			47
	0.3	0.4	0.45	0.5
	Poisson'

✁

s ratio

Table 3 :

 3 Optimality test for n = 75 subdomains

	h/H			1/3 1/6 1/10
	FETI (compressible)		33	37	39
	FETI-I (incompressible)	49	50	52
	FETI-IA (incompressible)	38	39	40
	FETI-DP (compressible)	20	23	26
	FETI-DPI (incompressible)	34	36	39
	FETI-DPIA (incompressible)	21	24	26
		60
					FETI-I
	Number of iterations	20 30 40 50	FETI-DPI FETI-DPIA	(FETI) compressible case FETI-IA
			(FETI-DP) compressible case
		10
		0
			1/10	1/6	1/3
				h / H

Table 4 :

 4 Numerical scalability test for h/H = 1/10

	H(n)			1/6 (12) 1/9 (27) 1/12 (48) 1/15 (75)
	FETI (compressible)			31	35	37	39
	FETI-I (incompressible)		38	46	49	52
	FETI-IA (incompressible)		34	36	38	40
	FETI-DP (compressible)		17	22	24	26
	FETI-DPI (incompressible)		24	33	37	39
	FETI-DPIA (incompressible)	17	22	25	26
		60
					FETI-I FETI-IA
	Number of iterations	20 30 40 50	FETI-DPI FETI-DPIA	(FETI) compressible case
			(FETI-DP) compressible case
		10
		0
		1/15 (75)	1/12 (48)	1/9 (27)	1/6 (12)
					H (n)

Table 5 :

 5 Comparison of convergences for different algorithms (NI = nearly incompressible, I = incompressible, C = compressible)

				NI problems		I problem
			ν = 0.49 ν = 0.499 ν = 0.4995	ν = 0.5
	NI preconditioner FETI-I	49	77	80
	NI preconditioner FETI-IA	37	39	39
	NI preconditioner FETI-DPI	52	68	69
	NI preconditioner FETI-DPIA	22	23	23
	I preconditioner	FETI-IA				40
	I preconditioner	FETI-DPIA				25
	C preconditioner	FETI-I	43	47	47	47
	C preconditioner	FETI-IA	34	36	37	37
	C preconditioner	FETI-DPI	31	34	34	34
	C preconditioner	FETI-DPIA	21	22	22	22

Table 6 :

 6 Comparison of complexities (in Mflop) for different algorithms

		number of iterations init. cost iter. cost total cost
		nit	c1	c2	c1 + c2nit
	FETI-I	47	151	6.24	445
	FETI-IA	37	395	6.69	642
	FETI-DPI	34	180	6.51	401
	FETI-DPIA	22	225	7.19	383

Implementation and complexity analysisA complexity analysis is a way to estimate the costs of algorithms; nevertheless it should be used only for an estimation of the main trends rather than a proof


Robakowska Dagmara
email: d.robakowska@yahoo.fr

Marie Hyzorek

Pierre Mirouze

Larmande

Dagmara Robakowska Hyzorek

Marie Mirouze

Pierre Larmande

Integration and Visualization of Epigenome and Mobilome Data in Crops

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction:

In the coming years, the study of the interaction between the epigenome and the mobilome is likely to give insights on the role of TEs on genome stability and evolution.

In the present project we have created tools to collect epigenetic datasets from different laboratories and databases and translate them to a standard format to be integrated, analyzed and finally visualized.

Visualization:

For data visualization we used Circos, a powerful standalone tool. The circular layout representations greatly enhance the visualization of scientific results.

Data Integration:

We aim to integrate our data in online epigenome and mobilome database for the rice and arabidopsis data.

• GBrowse graphical interface: 1 : pop up menu 2 : part of detailed report, both available for each feature.

  Fig1Fig1B

• Quick search: search through all databases and display search results as tables. Conclusion: Playing on different filter parameters and tools, this versatile workflow will be useful in the detection of active transposable elements. It has already been run on different data sets, from Arabidopsis thaliana and Oryza sativa.


Greener Aviation

M Forte

E Piot

J Perraud
email: jean.perraud@onera.fr

D Hue

J Duchemin
email: jeremy.duchemin@onera.fr

O Vermeersch
email: olivier.vermeersch@onera.fr

Research activities of ONERA on laminar airfoils in the framework of JTI Clean Sky SFWA-ITD: transition prediction

Keywords: laminarity, transition prediction, NLF, HLFC, surface defaults, surface tolerance

The present paper deals with recent research activities carried out in ONERA regarding laminar/turbulent transition prediction within JTI Clean Sky SFWA-ITD (Smart Fixed Wing Aircraft Integrated Technology Demonstrator). Two strategies were implemented to predict the location of the transition on laminar wings, airfoils or nacelles. Either the transition location is directly predicted with the elsA Navier-Stokes solver in which transition criteria are activated, or a mean flow is computed using a fully turbulent simulation and this data is used as an input for boundary-layer and stability computations (database approach or exact computations). This paper addresses also the issue of the impact of surface defaults on transition prediction, which is a prerequisite to increase the maturity of laminarity technology.

INTRODUCTION

Considering the continuously increasing demand in air transport, worldwide actions have been taken about environmental protection in the aviation industry. In this context, reduction of pollutant emissions remains one of the biggest issues raised by the International Civil Aviation Organization (ICAO). As less drag would imply a lower fuel consumption, numerous investigations aim at reducing the skin friction drag which usually reaches from 30% to 60% of the airplane total drag. Because of the rather large skin friction gap between laminar and turbulent boundary layers, such figures justify the scientific interest in Natural Laminar Flow (NLF) and Hybrid Laminar Flow Control (HLFC). These concepts aim at extending the laminar boundary layer region as far as possible from the leading edge by optimizing airfoil shapes and, in the latter case, by applying suction onto the laminar region. For more than three decades, the different scenarios to maintain a laminar boundary layer on wings are investigated by ONERA scientists, which developed dedicated tools for modeling these phenomena and predicting the transition location. Experimental activities have been conducted in the framework of national EUfunded research projects and numerical transitional models have been consolidated and implemented in our fluid solvers (elsA Navier-Stokes solver [START_REF] Cambier | The ONERA elsA CFD Software: Input from Research and Feedback from Industry[END_REF], 3C3D boundary layer code [START_REF] Perraud | Descriptif et mode d'emploi du code 3C3D[END_REF]). The present paper deals with recent research activities carried out regarding transition prediction within JTI Clean Sky and more specifically in the framework of the SFWA-ITD (Smart Fixed Wing Aircraft Integrated Technology Demonstrator). A summary of the main results obtained with ONERA tools are given on three industrial applications, and references are given to papers in which all the details can be found: a business-jet model equipped with laminar wings (collaboration with Dassault Aviation (DAv)), a two-dimensional airfoil model designed to study laminar boundary layer effect on the buffeting (with DAv) and a laminar nacelle (with SAFRAN Nacelles). Then, activities on the effect of twodimensional imperfections (such as steps, gaps or waviness) on the transition location are presented.

TRANSITION PREDICTION ON SMOOTH BODIES

Different strategies can be followed to predict the location of the transition in numerical simulations: the use of transition criteria in RANS simulations can give good results, provided that the near-wall mesh is fine enough to compute the boundary layer correctly, since those criteria are mainly based on the estimation of the integral parameters [START_REF] Cliquet | Application of Laminar-Turbulent Transition Criteria in Navier-Stokes Computations[END_REF]. Another strategy entails mean flow computation using a fully turbulent RANS simulation; this data is then used as an input for boundary-layer and stability computations, which give the amplification of boundary layer unstable waves until the laminar/turbulent transition threshold is reached [START_REF] Arnal | Laminarturbulent transition prediction in 3D flows[END_REF].

Business-jet model equipped with laminar wings

Within SFWA 213 workpackage, a Low Speed Business Jet equipped with laminar wings was designed by Dassault Aviation. A high speed wind tunnel test campaign was performed by Dassault Aviation at the European Transonic Windtunnel (ETW) in Cologne (Germany) in October 2015 (see Dassault Aviation's paper [START_REF]Design, Manufacturing and Testing at ETW of a Laminar Wing Business Jet Model at High Reynolds Number[END_REF]).

The main objectives of these tests were to assess the laminarity of the wings, the high speed performance of the U-shape horizontal tail plan, and innovative mock-up implementation.

ONERA and Dassault Aviation teams analyzed results in close collaboration. Especially, ONERA performed RANS computations with the elsA solver of the LSBJ configuration, the laminarturbulent transition being taken into account through the AHD-C1 transition criterion.

Transition criteria represent a compromise between accuracy and simplicity. In order to calculate the transition location induced by TS instabilities, the Arnal-Habiballah-Delcourt (AHD) criterion has been developed, while transition due to crossflow waves is captured by C1 criterion [6][7]. Figure 1 illustrates for instance one comparison between CFD (black lines) and WTT results (TSP images) of the transition location on the upper and lower surfaces of the EE wing in cruise conditions performed by Dassault Aviation and ONERA. CFD computations were performed before the tests. A very good match is obtained, and a journal paper is in preparation on this subject.

Figure 1: Comparisons between CFD and WT data for transition location

Laminar buffet and 2D laminar aerofoils

Within SFWA 112 workpackage, ONERA, in cooperation with Dassault-Aviation, analysed buffet onset for 2D aerofoils representative of the LSBJ configuration. Tests were performed in ONERA S2MA wind tunnel WT tests for natural and tripped transitions (see Fig. 2). Three different aerofoils designed by Dassault-Aviation for extended laminar flows were tested. The objectives of these tests were to evaluate the laminar behaviour of these aerofoils, to study buffet in laminar conditions and to acquire a detailed database for CFD software validation. RANS computations of the most interesting experimental configurations were performed with the elsA solver, the transition location being captured by the AHD criterion.

Effect of surface imperfections on transition location

In real applications, the surface of the aircraft is not perfectly smooth. Splices, rivets, manufacturing defaults can create surface imperfections which may deeply affect the transition process, and even ruin the positive effect of a laminar design. ONERA scientists have a long-time expertise in this topic and have compiled a lot of experimental data concerning the effect of two-dimensional imperfections such as steps, gaps or waviness on the transition location, either in a laboratory wind tunnel or during flight tests performed with industrial partners. From these experiments, ONERA proposed several transition triggering criteria, which are complementary to the famous Nenni & Gluyas criteria [START_REF] Nenni | Aerodynamic Design and Analysis of an LFC Surface[END_REF] and may thus allow more accurate surface tolerance requirements. ONERA also developed a ∆N approach [START_REF] Perraud | Laminar-Turbulent transition prediction in the presence of surface imperfection[END_REF], which models how the transition location may be progressively moved toward the leading edge depending on several surface imperfections parameters. SFWA activities were focused on gaps and holes.

Modeling of gap and holes effect by criterion and ∆ ∆ ∆ ∆N approach

From 2009 to 2015 experiments have been conducted in a research wind tunnel at ONERA in order to investigate the effect of gaps and holes on laminar-turbulent transition at low speed (U ∞ < 80 ms -1

). A two-dimensional model based on an ONERA-D profile has been successively equipped with several inserts including different surface defaults at 10% chord. Their effect on transition was investigated through hot wire measurements for sweep angles ϕ = 0° (TSinduced transition) and 60° (crossflow-induced transition). Some typical results are summarized below for the unswept configuration. The results are plotted using the critical Reynolds number R b defined by:

ν b U R c b =
, where b is the width of the gap and U c is the critical velocity for which transition is triggered. Figure 4 shows the evolution of R b versus the aspect ratio of the gap h/b. The dashed line corresponds to R b =15000, which is the critical value suggested by the Nenni-Gluyas criterion. Except for the low values of h/b, this criterion largely overestimates the Reynolds number corresponding to transition tripping. The results also show that the aspect ratio h/b = 0.2 is a value of interest: Re b seems to have a constant value (~ 6000) as long as the aspect ratio is higher than 0.2 thus, the assumption that depth is not a critical parameter (as assumed in the Nenni & Gluyas criterion) seems to be verified only for gaps sufficiently deep. The tripping limits are also plotted on the second plot of Figure 4: This plot shows the values of the ratio h/δ 1c versus the ratio b/δ 1c for each tested gap under the critical conditions, which means that δ 1c is the value of the displacement thickness at the location of the gap (x/c = 10%) and for the critical free-stream velocity. All the tripping limits gathered on this figure depict a specific L-shaped curve which separates the domain where the boundary layer is laminar from the domain where transition is tripped by the gap. This critical curve seems to have two asymptotes which are h/δ 1c ~ 2 and b/δ 1c ~ 18 meaning that a gap has to verify these two geometrical criteria to have a chance to trip transition. .

Figure 4 : Critical Reynolds Number R b versus aspect ratio and tripping limits for all tested gaps (φ=0°).

These experiments were also used to derive a ∆N approach dedicated to gaps. Steady laminar Navier-Stokes computations were undertaken using the ONERA code elsA, to compute the mean flow on a flat plate with a rectangular gap downstream of the leading edge. The boundary layer profiles were then used as input data for local, linear stability calculations to produce Nfactor curves. From the results of such parametric studies, it was shown that the N-factor modification induced by a gap results in the superposition of two components: on one hand a Npeak value which represents a drastic amplification just past the roughness and may induce transition onset at the gap position and, on the other hand, a Nfar value which represents a shift from the natural N-factor curve and which may induce a progressive upwind move of the transition position. All these results have been detailed in a conference paper [START_REF] Forte | Experimental and numerical study of the effect of gaps on laminar-turbulent transition of incompressible boundary-layers[END_REF] and in a journal paper [START_REF] Beguet | Modeling of transverse gaps effects on boundary-layer transition[END_REF].

In collaboration with Dassault-Aviation, the same kind of experiment was performed using isolated holes in order to address the lack of bibliographical data on this issue. A dedicated experimental campaign has thus been prepared, in order to study transition triggering by holes of various size (diameter and height) and tool prints in subsonic conditions.

Application to flight tests of ONERA methodology for surface tolerance prediction

From 2010, ONERA participated to the preparation and analysis of the HTP (Horizontal Tail plane) flight experiments on the F7X aircraft with Dassault Aviation, in the scope of SFWA 316 workpackage. In some of the flight tests, the HTP was perfectly smooth, in other flight tests calibrated imperfections were installed on the surface. These surface imperfections were placed parallel to the leading edge and could be considered as two-dimensional. Transition was detected with an infra-red camera. ONERA was especially involved into the qualification of surface quality. ONERA tools and methodology, similar to the one presented in the section before, were applied to the flight test configuration and results were compared to the experimental ones, with a specific focus on backward and forward facing steps. In order to take into account the effect of surface defaults on the amplification of instabilities inside the boundary layer, the ∆N method has been applied. Since two-dimensional defaults destabilized longitudinal waves, an envelope curve restricted to the amplification of only Tollmien-Schlichting instabilities has first been computed. In a second time, the modified N-factor (N+∆N) taking into account surface influence has been calculated, and a very good agreement was obtained with infra-red measurements.

CONCLUSION

ONERA was deeply involved into JTI Clean Sky SFWA activities dedicated to laminarity, in close collaboration with industrial partners. ONERA improved its numerical methodology and tools for laminar/turbulent transition prediction, which are based on the 3C3D boundary layer code and on the elsA Navier-Stokes solver. On smooth wings or airfoils, fast transition prediction location can be obtained with RANS simulations coupled with a transition criterion. If a more accurate transition location is needed, the database method or linear stability calculations have to be performed based on boundary layer simulations. Challenges that will be tackled by ONERA in the near future will be the implementation of the database method into the elsA solver, transition prediction on deformable wings, and optimal design of laminar wings. ONERA also has a specific methodology for predicting surface tolerance, based on a ∆N approach.

These tools will be continuously improved, for all kinds of surface imperfections that may exist in industrial configurations, taking into account the compressibility effects as accurately as possible in transonic flows.

Figure 2 :Figure 3 :

 23 Figure 2: Picture of the 2D laminar aerofoil within ONERA S2MA wind-tunnel

ACKNOWLEDGMENTS

Work presented in this paper has been undertaken within the Joint Technology Initiative "JTI Clean Sky", Smart Fixed Wing Aircraft Integrated Technology Demonstrator "SFWA-ITD" project (contract N° CSJU-GAM-SFWA-2008-001) financed by the 7th Framework Program of European Commission.


Pierre-Antoine Cuniasse
email: pierre-antoine.cuniasse@sncf.fr

Christine Buisson
email: christine.buisson@entpe.fr

Joaquin Rodriguez
email: joaquin.rodriguez@ifsttar.fr

Emmanuel Teboul
email: emmanuel.teboul@sncf.fr

David De Almeida
email: david.dealmeida@sncf.fr

Analyzing railroad congestion in a dense urban network through the use of road traffic network fundamental diagram concept

Transilien, the branch of SNCF in charge of operating the main urban railroad network in the area of Paris, faces a regular increase of passengers flows. The planning of railway operations is made carefully: simulation runs permit to assess the timetable stability. However, many disturbance appears and cause trains delays. Due to the nature of the railroad network those delays are cumulative and an on-line update of the timetable is not always successful in maintaining the trains schedule. In this tensed context, operators are searching solutions to better use the infrastructure capacity and better service quality. A needed step towards this objective is a better understanding of the phenomena of disruptions. In particular because the expansion of congestion is not clearly understood until now. This paper explores the possibility to transpose a traffic flow theory tool, the network fundamental diagram, in the field of dense railroad traffic. Railroad traffic is different of road traffic by many ways: railways are a planned system, traffic volume does not satisfy the continuum hypothesis, stations force stops and the signalization system brings a discrete behavior. Despite those big differences we show how to build a network fundamental diagram for a railroad system and how to interpret some obtained shapes for those diagrams. These diagrams gives us some means to compare planned and reality. We also identify the limits that need to be overcome to take benefits of the road traffic tools in railroad traffic analysis.

INTRODUCTION

In dense urban regions, the railway system is frequently operating not as scheduled. Indeed, it is a common fact that theoretical schedules of trains are not respected and that in reality, the train circulation system is operating differently than what is planned. Deviations from schedules are observed in two cases: 1/ exceptional situations with the occurrence of an external event (due to meteorological conditions for example) 2/ everyday situations where the schedule is not strictly respected (for example a small delay in one station induces a bigger in the next one). Whatever the cause of the delay, railroad users are impacted and this is decreasing their confidence into railway mode, destroying therefore the efforts of modal shift to a more energy-efficient mode.

To face this problem, railway companies add margins to the theoretical schedule to increase the probability of schedule adherence; however, the added margins reduce the availability of track capacity for additional rail transport services. [START_REF] Uic | Links between railway infrastructure capacity and the quality of operations[END_REF][START_REF]Capacity. International Union of Railways[END_REF][START_REF] Goverde | Punctuality of Railway Operations and Timetable Stability Analysis[END_REF]. Afterwards, the robustness of the obtained schedule is frequently evaluated with microscopic simulation tools [START_REF] Nash | Railroad simulation Using OpenTrack[END_REF][START_REF] Radtke | Handling of railway operation problems with RailSys[END_REF]. Nevertheless, those increased efforts towards more realistic theoretical schedules are not effective. This is even truer if, like in dense urban regions, both political orientations toward sustainable transportation and the economic grow cumulate in an increasing demand.

Therefore, the lack of positive results for those approaches and the absence of willpower for building an explanatory model leads us to search for another intellectual scheme to reduce the differences between observed and planned schedules in urban dense railway operations. Analysis tools initially developed for other congested transportation modes can be considered and their transfer to the railway context has to be examined. For the present research we choose the road transportation mode and proceed according to a three steps procedure: 1/ use the analysis methods applied to road traffic data to quantify the impact of congestion on the railway system, 2/ develop models inspired by road traffic models, 3/ make a cross comparison between modelling results and observations to test the validity of those models.

As a start, we consider a recently re-discovered tool for road traffic data analysis (the network fundamental diagram) and study its transfer to the case of the railway network of the Paris region, one of the most congested railway systems in Europe. The main aim of this paper is therefore to build the network fundamental diagram of railway lines (NFD-R) and to study the link between congestion and the shape of the NFD-R.

The paper is organized as follows. From a brief literature overview of road traffic network fundamental diagrams, we identify the main reported causes of dispersion of their points. A description of the railroad network and of the data collection technique is thus given. The results and their analysis is presented before a discussion. The paper ends with conclusions and recommendations for future research.

NETWORK FUNDAMENTAL DIAGRAM : STATE OF THE ART FOR ROAD TRAFFIC

NETWORKS

For more than half a century, road traffic is measured through local detectors (double electromagnetic loop detectors) that are able to measure at a given point the time between the passing of two successive vehicles and the speed of each of them. Macroscopic variables characterizing the traffic flow seen as a whole can be extracted from those observations. Those variables are the spatial mean speed and the flow (number of vehicles passing the detector during a given period). The density (number of vehicles for a given length of link) can also be extracted through a rough approximation of the fraction of time during which the loop is occupied by a vehicle.

Expanding those measurements to compute mean accumulation and flow of vehicles for a given period over a whole network is not a new idea. After the initial publication of this idea in [START_REF] Godfrey | The Mechanism of a Road Network[END_REF], Herman and Prigogine propose in [START_REF] Herman | A two-fluid approach to town traffic[END_REF] a formulation of the relationship between mean speed and mean concentration. This idea was further developed in many papers of those authors and collaborators, where the idea was to examine the relationship between the mean flow on a road network and its total accumulation. Those papers are devoted to the free flow part of the diagram: the mean flow increases with the accumulation. In particular [START_REF] Mahmassani | Investigations of Network-Level Traffic Flow Relationships: Some Simulations Results[END_REF] evokes the possibility of random events of lane blockages in a rectangular grid urban network simulated with NETSIM. The bigger the ratio between time with lane blockage and total simulation time is, the lower the slope of the relationship. Even with such lane blockages, it is worth noticing that the congested part of the relationship between global flow and global accumulation was not observed.

In (9) Daganzo hypothesized that above a given threshold in total accumulation (critical value of accumulation), an increase of this accumulation leads to a decrease of the total flow. This was experimentally observed for the first time by Geroliminis and Daganzo one year after [START_REF] Geroliminis | Existance of Urban-Scale Macroscopic Fundamental Diagrams: Some Experimental Findings[END_REF][START_REF] Geroliminis | Macroscopic Modeling of Traffic in Cities[END_REF] for the congested center of Yokohama city. This relationship between total flow and total accumulation averaged during periods of typically a few minutes on a (sub-) network was initially named by Daganzo and Geroliminis Macroscopic Fundamental Diagram. Following the authors of [START_REF] Saberi | Exploring the Properties of Network-wide Flow-Density Relations in a Freeway Network[END_REF] in their introduction, we do prefer the term "Network fundamental diagram" which is to our opinion more precise, "macroscopic" referring to the variables (flow and concentration) by opposition with the microscopic ones like inter-vehicular time or individual speed.

Since 2008, network fundamental diagram exploration, (driven either experimentally or analytically or through simulation data), has known many developments as a particular sub domain of road traffic flow studies. A recent and promising direction of research is the control of the most congested parts of a city network. It consists in splitting the city network into more than one reservoir and controlling the entry and exit flows of those reservoirs to decrease congestion (see for example [START_REF] Haddad | On the stability of traffic control in two-region urban cities[END_REF][START_REF] Keyvan-Ekbatani | Exploiting the fundamental diagram of urban networks for feedback-based gating[END_REF]). Recently, one of these groups of authors proposed a method to partition the network into homogeneous zones to reduce the scatter of the diagram and increase the network controllability [START_REF] Ji | On the spatial partitioning of urban transportation networks[END_REF]. Some authors study the impact of detectors types and location either on experimental data (see for example [START_REF] Buisson | Exploring the Impact of the Homogeneity of Traffic Measurements on the Existence of Macroscopic Fundamental Diagrams[END_REF]) or on analytical expressions (see [START_REF] Courbon | Cross-comparison of macroscopic fundamental diagram estimation methods[END_REF]). More generally, we have to stress that there is common agreement among various authors that the measurement methodology impact significantly the shape of the NFD. Some papers relate the scatter and /or the hysteresis near capacity point with the variability of density among links at a given instant. This can be derived from data [START_REF] Geroliminis | Properties of a well-defined macroscopic fundamental diagram for urban traffic[END_REF] or from simulation, in particular when looking at the traffic assignment impacts [START_REF] Knoop | Routing Strategies based on the Macroscopic Fundamental Diagram[END_REF].

Globally, as was first put into evidence in ( 16) and further on mentioned in [START_REF] Cassidy | Macroscopic Fundamental Diagrams for freeway networks : Theory and observation[END_REF], the NFD for a highway network is not well defined in the sense of [START_REF] Ji | On the spatial partitioning of urban transportation networks[END_REF]. Indeed, the congestion level of the highway network is not homogeneous. In recently published papers, the impact of combining various traffic states among various parts of the considered network for computing the NFD has been explored in parallel by two teams: [START_REF] Saberi | Exploring the Properties of Network-wide Flow-Density Relations in a Freeway Network[END_REF] and also by [START_REF] Daganzo | Macroscopic relations of urban traffic variables: bifurcations, multivaluedness and instability[END_REF].

To the best of our knowledge, the mains concepts of road traffic flow, in particular the network fundamental diagram, were seldom transposed in other congested modes of transportation. A paper present the adaptation of Personal Car Equivalent concept to take into account the difference in mechanical abilities among trains in computing the capacity of a railroad [START_REF] Lai | Development of Base Train Equivalents for Headway-Based Analytical Railway Capacity Analysis[END_REF]. We will now present the transposition of the concept of Network Fundamental Diagram in the case of the dense and highly congested railroad network surrounding Paris.

DATA AND METHODS Network

With a population of 11.5 millions inhabitants, Paris urban area is the most dense conglomeration in Europe. This density implies high levels of transportation demand and mass transit issues. This area is covered with two railroad networks: the subway and the RER network which transports about 7 millions of passengers a day (RER stand for Local Express Network). For this paper we choose to focus on the SNCF (french national railroad) operated part of the RER. Today, the growing demand makes the RER over-saturated. Figure 1 represents the RER network which is mostly unnested. In this paper, we focus on the north part of B line and the central part of C line.

The interest of those lines lies on their very high train traffic.

Data

Beacons are located over the railway network and detect the passage of trains and record their identifiers and events times. The beacon density over the network is variable and its order of magnitude is more than one beacon per station. Table 1 is a description of the two analysed areas: the northbound of B line between Paris gare du nord, CDG Aeroport terminal 2 and Mitry Claye, and the central section of C line between Champ de Mars and Choisy le Roi. The results of this paper are based on a SNCF database which gathers all the informations from the beacons and the corresponding theoretical times of passage over the beacons. Since each lines also got two directions, this gives four sets of data per day. Days with reduced demand (like week-end and holidays) or a special traffic event (planned public work...) are rejected. In the end, thirty-three days of data are available for B line and thirty-four for C line. The list of days used for each line is presented in table 2.

RESULTS AND ANALYSIS

First observations of the network fundamental diagram for railroad

For one day and one set of data, the number of beacon events over time and the number of trains were counted and filtered with a running average (window time: 10 minutes, step time: 1 minute). We consider that at any given instant, the total number of trains present in a line is directly proportional to the average train concentration. We further assume that the number of trains passing over all the beacons of a line during a given period is proportional to an average flow on the same section. Note that this is not completely true in our case, since the distance between two consecu-

FIGURE 1

The Paris-area RER network, where the SNCF parts are in color. The shaded area is the Paris inner city. The scale is not respected. tive beacons is not homogeneous over the line. Then, when plotting the number of events reported by all the beacons of a line against the number of trains of this line during the same period, we construct the NFD-R (Railroad) diagram.

Figure 3 presents the NFD-R for 2 lines: the B line -North direction (left column) and the C line -West direction (right column). To facilitate incident identification, markers were colored according to hours. The top diagrams of figure 3 represent the planned situation (T stands for Theoretical) on a typical day. The diagrams below present the results obtained in operation for various days. Depending on the day, a dispersion is observed or not. In any cases, most points are aligned along a (slopped) line. We now examine the slope of the diagrams and, if present, the possible causes of spreading.

We can see from figure 3 that various shapes were observed for real diagrams. We visually selected days and lines where the realized NFD-R was almost linear (for example, diagrams (B-R2), (B-R3), (C-R2) and (C-R3) from figure 3 were rejected). This set of NFD-R is further referred as the "filtered NFD-R". The numbers of filtered NFD-R are given in table 2. For those days and those lines, we made two polynomial regressions called y 1 and y 2 (see figure 4). Some of

TABLE 1

Details of analysed sections: the North-B line got a "Y" shape which is here described in three sections. This table provide the station name, the distance between consecutive stations, and the number of beacons in each station sector. This number can change according to the direction. List of days selected for analysis. The red-bold days are for B and C lines and the green-italic ones for C line only.

these regression are represented in figure 5. Since the second regression model y 2 is not linear, the R-value is not a good choice to compare the two relationships and the χ 2 residuals are a better choice. Table 2 presents the mean value of χ 2 for each cases; the two relationship gives very similar results. The linear regression was rejected because there is a strong reason for the curve to pass through the (0,0) point (zero

β 1 concentration flow y1 = β1.x β 2 γ 2 < 0 concentration y2 = β2.x + γ2.x 2 FIGURE 4
The two types of fit; y 1 is a proportional fit where the slope β 1 represents an average speed and is expressed in [measure.min -1 .train -1 ]; y 2 is a quadratic fit where γ 2 gives a curvature depending to its sign. A negative γ implies a concave curve which means that the average speed decrease when the train concentration increase.

trains implies zero beacon's measurements).

Comparison between theoretical and real NFD-R

It is well known that the slope of the left, uncongested part of a fundamental diagram is homogeneous to the speed in free flow conditions. To estimate a proxy for the railroad network, we use the β 1 value which is homogeneous to a speed. β 1 (from y 1 ) is then used to compare theoretical and realised NFD-R over all NFD-R (all diagrams must be included to evaluate the global performance of a day; the β 1 value represents in this case a performance indicator and not striclty the slope of a graph). Figure 6 represents the frequency of the difference between the β 1 values obtained for the actual and planned timetables (β realisedβ theoretical ) for the two studied lines. This figure shows that the trains of B line (resp. C line) go often more slowly than planned. C line appears also more constant than B line. The lower size of the B line filtered sample compared to the C line might be an explanation for the higher spreading. The lower values of the slopes for the realized NFD-R compared to the theoretical NFD-R reflects a lower operating speed of the system, compared to what is planned. This has to be linked with the results obtained in simulations and presented in [START_REF] Mahmassani | Investigations of Network-Level Traffic Flow Relationships: Some Simulations Results[END_REF], where the slope of the NFD of a urban road traffic network is lower when then duration of red phases of traffic lights is higher.

Coming back to filtered NFD-R, we can see in figure 7 that theoretical and real diagrams do not fit with the same shape when considering the quadratic regression model y 2 : γ 2 value are nearly always bigger for the theoretical timetables than for the timetables resulting from operations. A first exploration has shown that in some cases, the occurrence of those abnormal points is linked in time with unplanned events. We have to better establish this point and to look at the impact of the stations and the accelerations they generate on the NFD-R.

.train -1 ], γ = [min -1 .train -2 ] -6 -4 -2 0 2 •10 -2 0 10 20 60% β R 1 -β T 1 frequ B-line -3 -2 -1 0 1 •10 -2 0 10 20 89% β R 1 -β T 1 C-line

DISCUSSION AND CONCLUSION

We have shown in this paper that a network fundamental diagram (NFD-R) can be built for dense railroad Systems. This diagram is usually constructed for road traffic network both urban or freeways, and the literature is numerous where analysis of its shape is presented. Here we proposed a way to build it from identification by beacons of the passing of trains and we applied it to two lines where trains frequency is particularly high. With this diagram, We highlight the existence of congestion phenomena. For most of the days, the slope of the NFD-R, which is an equivalent of the average speed, is lower than planned. This congestion phenomenon can appear in several ways: the first is discussed in this article and consist in a simple decrease of the speed while a NFD-R still exhibits a linear shape. On some of the NFD-R, on the contrary, we observe another shape, with dots below the sloped line. Those points will be analyzed in a further publication. For the first category of days, we proposed two polynomials fits (one linear and one parabolic, both with a null value at origin). Regarding the residuals, they gives similar results. We observed that the quadratic model realized a small improvement over the linear one. Nevertheless, figure 7 shows a difference of curvatures between theoretical and real diagrams. We get positive values of γ (especially for theoretical NFD-R) when we expected it to be null or negative. An hypothesis is that the planning process leads to an optimistic timetabling where trains are asked to go faster when their concentration increases. If validated, a light on the timetabling process should be helpful to understand these counter-intuitive results.

The research presented in this paper was only a first step of a larger project. The research directions to explore now are numerous. Among them one can cite a better exploration of the data: a more precise study must be undertaken prior to eliminate the influence of variable beacons density and reinforce our results. An exploration of the data of other lines will be needed to confirm the validity of our first results. Also the database of unplanned events has to be examined jointly with the points of the NFD-R located below the slope line. Simulating at a large scale a railroad network with stations and realistic decelerations and acceleration will permit us to better understand the real NFD-R. We might also in a farer future explore through the NFD-R tool the feasibility of using cordon regulations to better cope with undesired congestion. The transposition of the recent findings of road traffic theory might help us in finding better ways of identifying and reducing congestion if the Paris suburban railroad network, provided that we keep in mind the

  FIGURE 3 examples of NFD-R. (B-T) & (C-T) represent theoretical diagrams, (B-R1), (B-R2), (B-R3), (C-R1), (C-R2) and (C-R3) represent real diagrams.

20 hoursyy 1 y 2 (B-T) β 1 =

 2021 FIGURE 5 examples of NFD-R, fitted with two mathematical relationship. (B-T) and (B-R): theoretical and real NFD-R for the B line in south direction. (C-T) and (C-R): theoretical and real NFD-R for the C line in east direction. β = [min -1 .train -1 ], γ = [min -1 .train -2 ]

FIGURE 6

 6 FIGURE 6Comparison of β value in y 1 fit for all NFD-R. For B line and for 60% of days β R < β T . For C line and for 89% of days β R < β T .

FIGURE 7

 7 FIGURE 7 Comparison of γ value in y 2 fit for filtered NFD-R

FIGURE 8

 8 FIGURE 8 local (a) and network (b) theoretical fundamental diagrams of railroad operations in absence of stations. The bold lines of figure a presents the various equilibrium states associated with the train operations blocking rules. Three states are presented; uncongested when trains operates with a free flow speed, congested states, where the speed is strictly limited, and the fully congested points where the trains are stopped. In figure b, three points illustrates the network points associated with two states. Point A is a combination of uncongested states, whereas points B and C are the combination of one uncongested and one congested state.

TABLE 2

 2 Average value for the residuals in different cases.

	Line direction	B south T B south R B north T B north R C east T C east R C west T C west R
	Days	16	16	14	14	24	24	24	24
	y˙1	0.0116	0.0129	0.0101	0.0102	0.0117	0.0162	0.0090	0.0136
	y˙2	0.0108	0.0127	0.0089	0.0101	0.0113	0.0160	0.0087	0.0126

Acknowledgements

This research is conducted with the benefit of a PhD grant from Agence Nationale de la Recherche Technologique, France. Authors want to thank Ludovic Leclercq and Winnie Daamen for fruitful discussions during the preparation of this paper.

significant differences among the two modes.