Christophe Kassiotis
email: christophe.kassiotis@enpc.fr

Adnan Ibrahimbegovic

Rainer Niekamp

Hermann G Matthies

partitioned algorithm, nonlinear stability proof and validation examples

Keywords: fluid-structure interaction, partitioned iterations, nonlinear stability proof

    Nonlinear fluid-structure interaction problem. Part I: implicit partitioned algorithm, nonlinear stability proof and validation examples

Introduction

Among multi-physics problems that are currently entering the mainstream of scientific research in computational mechanics (e.g. see [START_REF] Ibrahimbegovic | Engineering Structures Under Extreme Conditions: Multiphysics and Multi-scale Computer Models in Non-linear Analysis and Optimal Design[END_REF][START_REF] Oden | Research directions in computational mechanics[END_REF]), perhaps the most frequently studied are the problems of fluid-structure interaction (e.g.. see [START_REF] Barcelos | A Schur-Newton-Krylov solver for steadystate aeroelastic analysis and design sensitivity analysis[END_REF][START_REF] Bathe | A mesh adaptivity procedure for CFD and fluid-structure interactions[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Degroote | Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction[END_REF][START_REF] Deparis | Fluid-structure algorithms based on Steklov-Poincaré operators[END_REF][START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF][START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems[END_REF][START_REF] Charbel Farhat | Mixed explicit/implicit time integration of coupled aeroelastic problems: three-field formulation, geometric conservation and distributed solution[END_REF][START_REF] Felippa | Synthesis tools for structural dynamics and partitioned analysis of coupled systems[END_REF][START_REF] Felippa | Stabilization of staggered solution procedures for fluid-structure interaction analysis[END_REF][START_REF] Fernández | A Newton method using exact Jacobians for solving fluid-structure coupling[END_REF][START_REF] Förster | On the geometric conservation law in transient ow calculations on deforming domains[END_REF][START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF][START_REF] Heil | An efficient solver for the fully coupled solution of large-displacement fluid-structure interaction problems[END_REF][START_REF] Hübner | A monolithic approach to fluid-structure interaction using space-time finite elements[END_REF][START_REF] Joosten | Analysis of the block gauss-seidel solution procedure for a strongly coupled model problem with reference to fluid-structure interaction[END_REF][START_REF] Legay | An eulerian-lagrangian method for fluid-structure interaction based on level sets[END_REF][START_REF] Küttler | A solution for the incompressibility dilemma in partitioned fluid-structure interaction with pure Dirichlet fluid domains[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF][START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF][START_REF] Hermann | Algorithms for strong coupling procedures[END_REF][START_REF] Mehl | An Eulerian approach for partitioned fluid-structure simulations on Cartesian grids[END_REF][START_REF] Perić | Modelling fluid-induced structural vibrations: reducing the structural risk for stormywinds[END_REF][START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems-Part II: energy transfer analysis and three-dimensional applications[END_REF][START_REF] Michael | Treatment of acoustic fluidstructure interaction by localized Lagrange multipliers and comparison to alternative interface-coupling methods[END_REF][START_REF] Tezduyar | Space-time finite element techniques for computation of fluid-structure interactions[END_REF][START_REF] Tezduyar | Interface projection techniques for fluid-structure interaction modeling with moving-mesh methods[END_REF][START_REF] Tezduyar | Multiscale sequentially-coupled arterial FSI technique[END_REF][START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF][START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Wang | Fluid-structure interaction by the discontinuous-Galerkin method for large deformations[END_REF] among others). The fluid-structure interaction is already an interesting problem in its own right with a vast number of important applications. However, in this work we use it as the model problem for testing the novel paradigm of solution procedure based upon the direct coupling of different codes developed for a particular sub-problem (i.e. either solid or fluid mechanics) into a single code. In particular, we seek to provide the guarantees for the robustness of such a computation approach in fully nonlinear setting, where implicit schemes are used for each sub-problem, and we derive (by a formal proof) the convergence criterion for partitioned scheme iterations.

We note in passing that the current tendencies for solving fluid-structure interaction problems are mostly oriented towards the monolithic schemes, where both sub-problems are discretized in space and time in exactly the same manner resulting with a large set of (monolithic) algebraic equations to be solved simultaneously with no need to distinguish between the "fluid" and the "structure" part. To provide the unified discretization basis for monolithic approach, the most frequent choice is to use the stabilized finite elements for fluids (first proposed by Hughes and co-authors [START_REF] Franca | Stabilized finite element methods[END_REF], and followed by [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF][START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF][START_REF] Legay | An eulerian-lagrangian method for fluid-structure interaction based on level sets[END_REF][START_REF] Hübner | A monolithic approach to fluid-structure interaction using space-time finite elements[END_REF][START_REF] Perić | Modelling fluid-induced structural vibrations: reducing the structural risk for stormywinds[END_REF][START_REF] Tezduyar | Space-time finite element techniques for computation of fluid-structure interactions[END_REF][START_REF] Tezduyar | Interface projection techniques for fluid-structure interaction modeling with moving-mesh methods[END_REF][START_REF] Takizawa | Wall shear stress calculations in spacetime finite element computation of arterial fluidstructure interactions[END_REF][START_REF] Wang | Fluid-structure interaction by the discontinuous-Galerkin method for large deformations[END_REF]), which can easily be combined with standard finite elements for nonlinear structure mechanics (e.g. see [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Olgierd | The Finite Element Method[END_REF][START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF]); less frequently used possibility is to provide the finite volume scheme discretization of structure sub-problem [START_REF] Mehl | An Eulerian approach for partitioned fluid-structure simulations on Cartesian grids[END_REF] and thus provide the monolithic basis constructed by finite volume method.

We depart from these trends of using monolithic approach, and we set to develop the partitioned approach for nonlinear fluid-structure interaction. We gain in this manner the first important advantage of using standard discretization scheme that is the most suitable for a particular sub-problem. The chosen model problem will consider finite volume method for fluid and finite element method for structure. We hope in this manner to gain in our approach generality with mixing different discretization methods, as well as combining the existing codes already developed for each sub-problems.

We also note that the partitioned approach was the favorite choice of early developments (e.g. see [START_REF] Belytschko | Mixed methods for time integration[END_REF][START_REF] Felippa | Synthesis tools for structural dynamics and partitioned analysis of coupled systems[END_REF][START_REF] Felippa | Stabilization of staggered solution procedures for fluid-structure interaction analysis[END_REF][START_REF] Hughes | Implicit-explicit finite elements in nonlinear transient analysis[END_REF]). A number of those developments were motivated by the efficiency, which was gained by treating a (large) part of the fluid-structure set of equations by explicit schemes (see [START_REF] Belytschko | Mixed methods for time integration[END_REF][START_REF] Hughes | Implicit-explicit finite elements in nonlinear transient analysis[END_REF]). The stability of such an approach is proved mostly for linearized stability criteria (e.g. see [START_REF] Felippa | Synthesis tools for structural dynamics and partitioned analysis of coupled systems[END_REF][START_REF] Felippa | Stabilization of staggered solution procedures for fluid-structure interaction analysis[END_REF]), which still leaves the potential danger of appearance of so-called "arrested" instability (e.g. see [START_REF] Belytschko | An overview of semidiscretization and time integration procedures[END_REF]), where the disproportionate increase of response in the unstable regime is stopped by inelastic response of the structure that softens, but not before a non-negligible error is introduced.

We will eliminate herein any of these problems related to arrested instability by using the implicit schemes non only for each sub-problem, but also for the partitioned coupling. The resulting algorithm will thus provide the computational robustness, which will be a great interest not only for fluid-structure interaction but also for any other multi-physic problem of current interest where one would like to re-use the available codes for any particular sub-problem in a more general framework.

The main advantage of code-coupling approach for fluid-structure interaction concerns the fact that the coupling is limited only to the fluid-structure interface. Therefore, the main difficulty is reduced to enforcing the interface matching with respect to two different discretization schemes, finite element versus finite volume, as well as two different time integration schemes and different time steps. We thus split the presentation of our work in two parts, pertaining respectively, to time and to space discretization for fluid and for structure and their matching at the interface. We will deal with the interface matching for different space discretization in Part II, along with other related issues pertaining to the computational efficiency enhancements by nested parallelization. In present paper (Part I), we discuss how to accommodate any particular (implicit) scheme that ensures the unconditional stability for either fluid or structure motion computation, and how to ensure that the unconditional stability extends to partitioned solution of the fluid-structure interaction problem. The outline of the paper is as follows. In Section 2 we give the theoretical formulation of the chosen interaction problem of the incompressible fluid and the finite deformation elastic structure. The partitioned solution algorithm for this interaction problem is presented in Section 3, and the formal proof of the stability in fully nonlinear setting in Section 4. In Section 5, we present the results of the benchmark examples and in Section 6, some closing remarks. The complementary presentation of the space discretization matching issues and nested parallel computing are discussed in Part II of this paper.

Theoretical formulation of fluid-structure interaction problem

We consider the interaction problem of motion of fluid (denoted further with f ) and structure (denoted with s). For the sake of generality, it is supposed that theses problems are nonlinear and time dependent. This interaction is handled by the classical Direct Force-Motion Transfer (DFMT, see [START_REF] Felippa | Synthesis tools for structural dynamics and partitioned analysis of coupled systems[END_REF][START_REF] Michael | Treatment of acoustic fluidstructure interaction by localized Lagrange multipliers and comparison to alternative interface-coupling methods[END_REF]) that can formally be expressed in terms of the corresponding Steklov-Poincaré operator (see [START_REF] Deparis | Fluid-structure algorithms based on Steklov-Poincaré operators[END_REF]):

S i : H 1 2 (Γ ) → H -1 2 (Γ ) u i → λ i with i ∈ { f , s} (1)
The Steklov-Poincaré operator gives the evolution of dual field λ i for an imposed primal field u i on the interface space and time domain Γ × [0, T ]. In fact, this operator requires the computation of the fluid and structure problem on the complete space-time domain, governed by Navier-Stokes equations for incompressible flow and nonlinear dynamics equations for structural motion. For the problems of this kind, the Steklov-Poincaré operator is not available analytically, but rather as the results of numerical approximations. These results not only depend on the chosen model, material properties and boundary conditions, but also on discretization techniques, time integration algorithms, discrete equations solvers, etc.

The interface matching considered in the following is based on two classical mechanics principles:

(i) Continuity of primal quantities or the perfect matching condition over the interface:

u s = u f = u on Γ × [0, T ] (2)
where u denotes the value of the primal variable at the interface. In the continuum setting, the time derivatives of this condition give the equivalent equations for velocity v = u and acceleration a = v:

v s = v f = v; a s = a f = a on Γ × [0, T ] (3)
However, as the result of time discretization, these conditions are no longer equivalent. (ii) Equilibrium of dual quantities or action-reaction principle, which implies:

λ f + λ s = 0 on Γ × [0, T ] (4)
The action-reaction principle in (4) can be reformulated using the Steklov-Poincaré operator defined in (1), resulting with the so-called Steklov-Poincaré formulation (see [START_REF] Deparis | Fluid-structure algorithms based on Steklov-Poincaré operators[END_REF]):

Find: u on Γ × [0, T ] , so that: S f (u) + S s (u) = 0 (5)
Using the inverse of the first Steklov-Poincaré operators allows to rewrite the equilibrium of dual quantities as the following fixed-point equation that concerns only the unknown at the interface:

Find: u on Γ × [0, T ] , so that: u = S s -1 -S f (u) (6)
We note that the inverse of Steklov-Poincaré operator S -1 is no more expensive to compute than S . The fixed-point equation can be reformulated in order to get a root equation:

Find: u on Γ × [0, T ] , so that: S s -1 -S f (u) -u = 0 (7)
where all the requested quantities are defined at the interface. We recall that these Steklov-Poincaré operators are associated with the semi-discrete form of the continuum equations. Namely, the fluid problem is defined by an ALE formulation (e.g. see [START_REF] Demirdžić | Space conservation law in finite volume calculations of fluid flow[END_REF][START_REF] Hughes | Lagrangian-Eulerian finite element formulation from incompressible viscous flows[END_REF]) of the Navier-Stokes equations. The latter considers the fluid mesh motion that u m is imposed by the motion of the interface u and can be written as:

R m (u m ; u) := K m u m -D m u = 0 (8)
where D m is a projection/restriction operator and K m governs the extension of the boundary displacement either by a diffusion process or a pseudo-solid equation (e.g. see [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]). The fluid flow in this moving domain is described by the semi-discrete Navier-Stokes equations:

R f (v f , p f ; u m ) := M f v f + N f (v f -um )v + K f v f + B f p f -f f B T f v f = 0 (9)
where v f , v f and p f are fluid velocity, its derivative and fluid pressure M f is a (positive definite) mass matrix, N f is a (non-symmetric) advection matrix, K f is the matrix with diffusion terms, B f stands for the gradient matrix, and f f is the driving force on the flow. The discretization process leading to fluid equations of motion and incompressibility constraint in ( 9) is carried out by Finite Volume Method [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. The Finite Element Method [START_REF] Olgierd | The Finite Element Method, The Basis[END_REF] is used for the structure sub-problem, resulting with structural semi-discrete equations of motion:

R s (u s ; λ ) := M s üs + f int s (u s ) -f ext s (λ ) = 0 ( 10
)
where M s is the mass matrix, f int s is the internal force which is highly non-linear if large deformation or complex material behavior is used (e.g.. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF]), and f ext s is the external force vector.

Algorithm 1 Steklov-Poincaré operator for fluid: S f Require: Fluid state variable at time T N 1: Impose displacement of mesh at fluid-structure interface:

u m = u on Γ × [T N , T N+1 ] 2: Solve mesh intern nodes displacement: R m (u m ; u) = 0 on Ω f × [T N , T N+1 ] 3: Solve fluid problem in ALE formulation: R f (v f , p f ; u f ) = 0 on Ω f × [T N , T N+1 ] 4: Get boundary traction force at the interface: λ = -σ f n on Γ × [T N , T N+1 ]
With this notation on hand we can provide the explicit form of the Steklov-Poincaré operator for fluid and for structure, which can be expressed in terms of Algorithm 1 and 2. 3 Discrete force motion transfer (DFMT) -Explicit versus implicit scheme for interface matching

We assume in the subsequently that each sub-problem (either fluid flow or structural motion) has been integrated by an implicit time-stepping scheme (wide variety is available in [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF][START_REF] Ibrahimbegovic | Nonlinear solid mechanics: Theoretical formulations and finite element solution methods[END_REF]). We are thus only left with integrating the interface matching unknowns. Each of Steklov-Poincaré operators involves the complete domain for either fluid and solid, and furthermore, each results in a set of nonlinear equations to be solved. Therefore, an iterative scheme has to be defined to compute the solution of the coupled problem. To that end, one can use any of the formulation in ( 5), ( 6) or [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF] and corresponding iterative schemes, such as: i) the one in [START_REF] Belytschko | An overview of semidiscretization and time integration procedures[END_REF] with Picard like iterative schemes [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF]; ii) the one in [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF] with Newton [START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF][START_REF] Fernández | A Newton method using exact Jacobians for solving fluid-structure coupling[END_REF] or quasi-Newton [START_REF] Gerbeau | A quasi-Newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows[END_REF][START_REF] Tezduyar | Space-time finite element techniques for computation of fluid-structure interactions[END_REF][START_REF] Degroote | Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction[END_REF] iterations; iii) the one in (5) with non-linear Richardson strategy [START_REF] Deparis | Fluid-structure algorithms based on Steklov-Poincaré operators[END_REF][START_REF] Fernández | Domain decomposition based Newton methods for fluid-structure interaction problems[END_REF].

Generalized Conventional Serial Staggered (GCSS) algorithm

This explicit approach, yet called weak coupling, tries to provide the best possible approximation of the interface matching solving only once Steklov-Poincaré operator for each sub-problem and in each time step. More precisely, we consider the evolution of the fluidstructure interaction problem in a time-interval (or window) [T N , T N+1 ] of size ∆t. The idea is to solve at each time step a single Picard-iteration of the fixed-point equation in [START_REF] Belytschko | An overview of semidiscretization and time integration procedures[END_REF]:

u N+1 = S s -1 -S f (u N ) (11)
However, this leads to bad conservation properties at the interface, which can be improved by the addition of a better predictor P. This results with the Algorithm 3.

Algorithm 3 Generalized Conventional Serial Staggered 1: Given: initial time T = T 0 , final time T max , window size ∆t, initial interface displacement u 0 . 2: while T < T max do 3: Predict displacement:

u P N+1 = P(u N , uN , u N-1 , . . . ) 4:
Solve problem f:

λ N+1 = S f (u P N+1 ) 5:
Solve problem s: u N+1 = S s (λ N+1 ) 6:

N ← N + 1 and T ← T + ∆t 7: end while Note that we do not impose the way the time integration is performed, but only that the fluid and the structure part be collocated at the end of the same time windows. When considering equal time step size for fluid and structure, this Direct Force-Motion Transfer algorithm is named Conventional Serial Staggered (DFMT-CSS) (see [START_REF] Charbel Farhat | Mixed explicit/implicit time integration of coupled aeroelastic problems: three-field formulation, geometric conservation and distributed solution[END_REF][START_REF] Felippa | Synthesis tools for structural dynamics and partitioned analysis of coupled systems[END_REF]). We can also consider the so-called Sub-cycled Conventional Staggered Scheme (DFMT-SCSS) where time steps selected for integration of fluid flow and structure motion are not the same size. Usually, the characteristic time scale of the fluid is smaller, and thus we consider integration for the fluid part with many small time steps on the window [T N , T N+1 ].

The interface displacements for the structure at time T N+1 is therefore u N+1 while for the fluid, the final value is set by the predictor as a function of u N . For that reason, the interface matching condition at T N+1 , cannot be fulfilled by this explicit scheme and we obtain:

u s,N+1 = u N+1 = u f ,N+1 = u P N+1 ( 12
)
Moreover there is no reason that the interface matching condition be fulfilled for the velocity. As a result, an energy error is introduced by this kind of exchange at the interface which can be estimated by computing the energy transfer to the interface for both fluid and structure [START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems-Part II: energy transfer analysis and three-dimensional applications[END_REF][START_REF] Farhat | Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems[END_REF].

The first analysis of this explicit DFMT partitioned strategy [START_REF] Felippa | Stabilization of staggered solution procedures for fluid-structure interaction analysis[END_REF] printed out the upper limit on the time step size beyond which each numerical simulation of structure motion explicitly coupled to fluid flow diverges. This upper limit depends on the fluid/structure density ratio and on the speed of sound in the fluid medium. This criterion directly applied to incompressible flows with infinite wave speed, predicts immediate instability for any chosen time step. The instability phenomenon are noticed in [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF][START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF][START_REF] Daniel | Accelerated iterative substructuring schemes for instationnary fluid-structure interaction[END_REF]. and a simplified model is proposed [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF] for their prediction. The instability of explicit interface matching is clearly shown by the numerical examples in Sec. 5.

Direct Force-Motion Transfer Block Gauss-Seidel (DFMT-BGS) algorithm for implicit matching

By enforcing the continuity of primal variables at the interface we can eliminate the energy errors that characterize the explicit interface matching. This ought to be done by iterating on the following residual to reduce its value below the chosen tolerance:

r N+1 := u s,N+1 -u f ,N+1 ≃ 0 ≤ TOL ( 13
) t t ∆t s ∆t f P f P s T N T N+1 ×(k N ) λ | [T N ,T N+1 ] u| [T N ,T N+1 ]
Fig. 1 Block Gauß-Seidel coupling algorithm for fluid (P f ) structure (P s ) interaction problems; this iterative scheme is applied (k N ) times until convergence on a window [T N , T N+1 ].

In this way we obtain an implicit algorithm requiring more than one iteration to enforce the interface matching condition. The chosen order of iterations, as presented in Fig. 1, corresponds to the Block-Gauß-Seidel algorithm for fluid-structure interaction problem. We also show in Fig. 1 that not only the value at synchronization points T n or T n+1 , but also the interpolated evolution of variables have to be exchanged in the entire time

-interval t ∈ [T n , T n+1 ].
Contrary to explicit algorithms which generate spurious energy at the interface, the present implicit interface matching algorithm enforce the same evolution of the primal variables at the fluid-structure interface. In other words, an iterative solution for equation ( 6) is performed by using the Picard iteration:

u (k+1) N+1 = G u (k) N+1 ; G = S s -1 • -S f ( 14
)
where S f and S s are Steklov-Poincaré operators for fluid and structure defined previously.

The Picard iterations will continue until convergence of interface residual is achieved:

r (k) N+1 = u (k) s,N+1 -u (k) f ,N+1 = G u (k) N+1 -u (k) N+1 ( 15
)
Such a Block-Gauß-Seidel algorithm with implicit interface matching, further denoted as DFTM-BGS, can be presented as a natural generalization of the explicit algorithms. We can thus write:

Algorithm 4 Direct Force-Motion Transfer Block-Gauß-Seidel 1: Given: initial time T = T 0 , final time T max , window size ∆t, initial interface displacement u 0 . 2: while T < T max do 3: (k) = 0 4: Predict displacement: u (0) N+1 = P(u (kmax) N , u(kmax) N , u (kmax) N-1 , . . . ) 5: repeat 6: Perform Picard iteration: G u (k) N+1 7: Compute residual: r (k) N+1 = G u (k) N+1 -u (k) N+1
8: Update interface primal variable: u

(k+1) N+1 = u (k) N+1 + r (k) N+1 9: do (k) ← (k) + 1 10: until r (k-1) N ≥ TOL 11:
N ← N + 1 and T ← T + ∆t 12: end while As illustrated in Figure 2, this fixed-point algorithm based on Picard iterations for the time has the main drawback that the search directions for u and λ variables do not exploit any information from the fixed-point function G nor the Steklov-Poincaré operators S f and S s . Therefore, quite a few iterations may be needed to reach the convergence.

In order to improve the convergence of the DFMT-BGS method, we can use a more accurate update: u

(k+1) N+1 = u (k) N+1 + H r (k) N+1 ( 16
)
The choice of matrix H above should be made to improve the method convergence. Our favorite choice for constructing H is using a secant methods which can keep the cost of each iteration as low as possible. In particular, an approximation of H by a scalar ω employed herein with optimal value at each iteration ω (k) obtained by Aitken's relaxation. For such a choice the convergence rate is able to reach α ≃ 1.6.

e (k) u u (k-2) u (k-1) u (k) u (k+1) λ λ (k-1) λ (k) λ (k+1) -S f S s
Fig. 2 Implicit for interface matching by non-linear Block Gauß-Seidel.

We note in passing that various other improvements for constructing H matrix have been attempted (e.g. see [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Hermann | Algorithms for strong coupling procedures[END_REF][START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Degroote | Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction[END_REF]), but we did not find them indispensable for the studies presented herein where dynamic Aitken relaxation proved to be sufficient. Moreover, for such a choice we are able to provide the formal proof of stability of DFMT-BGS in fully nonlinear framework.

Proof of stability and convergence of DFMT-BGS

Reformulation of the fluid-structure interaction problem in a Differential Algebraic Equation framework

The main goal of this section is to confirm the stability of the partitioned approach for coupled fluid-structure interaction problems with implicit interface matching based upon the DFMT-BGS iterative scheme. At the outset, we assume that the approximate choice of stable and accurate integration schemes is made to solve for fluid flow and structure motion 1 . We can thus focus only upon the interface matching computations. To that end, we will first recast the DFMT-BGS algorithm formulated in the general framework of differentialalgebraic equations (DAE). This allows to follow in the footsteps of the proof of nonlinear stability of the partitioned algorithm provided in [START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF].

The first sub-problem pertaining to the fluid flow on a moving domain, discretized by Finite Volume Method, is defined jointly by equations ( 8) and ( 9), which can be written as:

0 = r f (x f (t), x s (t), y f (t), λ (t)) =   K m u m -D m u s M f v f + N f (v f -um )v + K f v f + B f p f -f f -D T f λ B T f v f   (17)
In ( 17) above v f is fluid velocities, and u m is the mesh displacements, which are gathered together in x f = (v f , u m ), while their time derivatives as well as the pressure field are placed in y f = (v f , p, um ). The forces at the interface are denoted with λ and matrix D f is the result of their interpolation across the interface from the fluid side. The residual r f gathers the mesh motion for ALE, the coupled discretized momentum equation and incompressibility condition. Rewriting the incompressibility condition using the acceleration in order to reduce the order of the DAE associated with incompressible fluid problem leads to:

0 = r f (x f (t), x s (t), y f (t), λ (t)) =    K m um -D m us M f v f + N f (v f -um )v + K f v f + B f p f -f f -D T f λ -B T f M -1 f N f (v f -um )v + K f v f + B f p f -f f -D T f λ    (18)
Considering now the structure motion, we will gather the primal variables, the solid displacements u s and velocities us , in the same vector x s = (u s , us ), whereas the acceleration is denoted as y s = üs . The residual form of the equation of motion can then be written:

0 = r s (x s (t), x f (t), y s (t), λ (t)) (19)
The last equation can be written explicitely as:

r s (x s (t), x f (t), y s (t), λ (t)) := M s üs + f int s (u s ) -f ext s -D T s λ (20)
where the D s indicates the force distribution at the interface on the structure side. The interface matching equation, which accounts for interaction, can be stated in terms of acceleration:

0 = r λ (x f (t), x s (t), y f (t), y s (t)) := -D s üs + D f v f (21)
With this notation in hand, the proposed DFMT-BGS algorithm with implicit interface matching can be stated as follows: first solve the fluid sub-problem together with the continuity equation at the interface:

∂ t x f (k) = f r f (x f (k) , x s (k-1) , y f (k) ) 0 = r f (x f (k) , x s (k-1) , y f (k) , λ (k) ) 0 = r λ (x f (k) , x s (k-1) , y f (k) , y s (k-1) ) (22)
and then solve the structure sub-problem with imposed forces at the interface:

∂ t x s (k) = f r s (x f (k) , x s (k) , y s (k) ) 0 = r s (x f (k) , x s (k) , y s (k) , λ (k) ) (23)

Error propagation, stability and convergence of DFMT-BGS algorithm

In the examples to follow we solved each nonlinear sub-problem by a dedicated solver (e.g. conserving algorithm for the structure (see for instance [START_REF] Ibrahimbegovic | Energy conserving & decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations[END_REF]), segregated approach like PISO for the fluid part [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]). The final stability results is valid in a more general context as long as sub-problem computation remains stable and convergent. Thus we only need to confirm the convergence of the iterative DFMT-BGS procedure for the interface matching. We introduce the following notations: the subscript N for the restriction of the function of time t to the interval [T N-1 , T N ], h = T N -T N-1 , the superscripts (k) for the iteration counter of the DFMT-BGS for interface matching procedure and the ∆ symbol denoting the distance of numerical approximation to the exact solution, ∆ x

(k) N = x (k) N -x ⋆ .
Moreover, the same time steps for each sub-problem of the fluid-structure interaction are not needed; for instance the fluid problem can be solved with many smaller time steps as opposed to a single time step for structure. The time interpolation of the evolution on the window has to be considered a priori.

The stability of DFMT-BGS operator split procedures is mainly governed by error propagation from one window, say N, to the following N + 1. In [START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF], a stability study of the partitioned approach for differential-algebraic systems of equations (DAE) that characterize the multi-body system dynamics is given. The stability proof given in [START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF] considers the case Lagrange multipliers or dual, rather than primal variable values are exchanged at the interface from one iteration to the next. In [START_REF] Arnold | Constraint partitioning in dynamic iteration methods[END_REF], the stability criterion corresponding to the block-Gauß-Seidel algorithm such as the one used herein is given [START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF], but without any proof. In the following a detailed proof for DFMT-BGS algorithm that cannot be found in the literature for the fluid-structure interaction context is given.

We note in passing that the same kind of results can be used for other coupled problems where the time integration schemes of the sub-problems are different, such as thermomechanics [START_REF] Kassiotis | Multiscale in time and stability analysis of operator split solution procedure applied to thermomechanical problems[END_REF] or with different time scales for mechanics and thermal component as well as a generalized non-linear operator split for problems with internal variables [START_REF] Kassiotis | Stable splitting scheme for general form of associated plasticity including different scales of space and time[END_REF]. However, the proof given herein is entirely original, as the dual quantities are exchanged in each iteration.

The main result can be stated as follows:

Theorem 1 There exists C ∈ | + such that for all (k) > 1 max n ∆ x f (k) n + ∆ x s (k) n + ∆ y f (k) n + ∆ y s (k) n + ∆ λ (k) n < C • max n µ k-2 ∆ x f (0) n + ∆ x s (0) n + µ k-1 ∆ y s (0) n
where: µ = α + O(H) and α the contraction constant characterizing the operator split procedure employed to solve this nonlinear interaction problems, which can be written as:

α = max n ∂ y s r λ ∂ y s r s -1 ∂ λ r s -1 ∂ y f r λ ∂ y f r f -1 ∂ λ r f
The main condition on contraction property with α < 1 will guarantee the convergence of the operator split procedure for the given time window when (k) -→ ∞. We set now to apply this stability criterion to the DFMT-BGS for fluid-structure interaction presented in the previous section. We thus obtain:

(i) Fluid sub-problem on a moving domain discretized with FV method:

∂ y f r f =   K m 0 0 ∂ um (N f (v f -um )v) M f B f -B T f M -1 f ∂ um (N f (v f -um )v) 0 -B T f M -1 f B f   (24)
and

∂ λ r f =   0 -D T f B T f M -1 f D T f   (25)
(ii) Matching condition that corresponds to the continuity at the interface:

∂ y s r λ = -D s ; ∂ y f r λ = 0 D f 0 (26)
(iii) Structural problem solved with a FE method:

∂ y s r s = M s ; ∂ λ r s = -D T s ( 27
)
Combining the given corresponding values of Jacobian computed for all these equations leads to the corresponding value of the contraction constant α for fluid-structure interaction, that can be written after simplification:

α = max n ∂ y s r s -1 ∂ λ r s ∂ y f r λ ∂ y f r f -1 ∂ λ r f -1 ∂ y s r λ = max n M -1 s D T s D f M -1 f 1 + B f B T f M -1 f B f -1 B T f M -1 f D T f -1 D s
This criterion shows that the numerical computation is highly linked to the material properties and more generally to the model chosen for the fluid and the structure subproblems. As the first approximation an estimate for α can be provided by the mass ratio between the fluid M f and the solid part M s weighted by some geometrical conditions for the field transfer defined in D f and D s . When the mass ratio increase, the scheme can become unstable. The terms in B f show the influence of the incompressibility condition, resulting in the stability domain reductionand through the added-mass effect (see [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF][START_REF] Förster | On the geometric conservation law in transient ow calculations on deforming domains[END_REF][START_REF] Tallec | Fluid structure interaction with large structural displacements[END_REF]). In order to the domain of stability for iterative DFMT-BGS procedure, Aitken's relaxation can be used as presented in the previous section. The stability of such modified algorithm can then be proved following the recepies from [START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF] for pre-conditioned algorithm. This proof will not be given herein, but only the corresponding improvement shown in the numerical examples.

Proofs for stable error propagation

The proof of the stability error propagation theorem is given as follows:

(i) The error bound for one iteration of the operator split procedure for an approximation in the neighborhood of the solution is given in Lemma 1. (ii) In Lemma 2, a recursive application of Lemma 1 gives a bound for (k) iteration of the operator split procedure. (iii) The proof is concluded by the application of the two Lemmas with suitable arguments.

Lemma 1 Consider U γ 0 a neighborhood of the solution x f ⋆ , x s ⋆ , y s ⋆ :

U γ 0 = x f , x s , y s x f -x f ⋆ + x s -x s ⋆ + y s -y s ⋆ γ 0
There exists (C, H 0 , γ 0 ) ∈ | 3

+ such that:

∀ (x f (0) , x s (0) , y s (0) ), ( x f (0) , ỹs (0) , λ (0) ) ∈ U γ 0 2 , ∀H < H 0 ,   δ x f (1) δ x s (1) δ y s (1)     CH CH CH CH CH CH C C α +CH     δ x f (0) δ x s (0) δ y s (0)   +   δ x f (0) (T n ) δ x s (0) (T n ) 0  
with δ denoting the distance between two approximations, i.e. δ x (k) = x (k)x(k) and

α = α + O(1) ∆ x f (0) + ∆ x s (0) + ∆ x s (0) + ∆ x f (0) + ∆ xs (0) + ∆ xs (0)
Proof (of Lemma 1) By inserting (x f (0) , x s (0) , y s (0) ) ∈ U γ 0 in the proposed algorithm, we obtain the evolution equation of the fluid sub-problem [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF] as the first step of the staggered scheme: 1) , x s (0) , y f (1) ) ) , x s (0) , y f (1) , λ (1) ) ) , x s (0) , y f (1) , y s (0) )

     ∂ t x f (1) = f r f (x f (
0 = r f (x f ( 1
0 = r λ (x f ( 1
with x f (1) (T N ) = x f (0) (T N ) (28)
For the second sub-system in [START_REF] Förster | On the geometric conservation law in transient ow calculations on deforming domains[END_REF], the evolution of the structure sub-problem under the loading λ (k) is recovered: ) , x s (1) , y s (1) ) ) , x s (1) , y s (1) , λ (1) )

∂ t x s (1) = f r s (x f ( 1
0 = r s (x f ( 1
with x s (1) (T N ) = x s (0) (T N ) (29)
The same kind of evolution is found for any other initial value in the neighborhood of the solution:

( x f (0) , xs (0) , ỹs (0) ) ∈ U γ 0 .
With a regularity assumption which guarantees that f r f and f r s satisfy Lipschitz conditions w.r.t. to their arguments (x f , x s , y f , y s ), 2 , we can obtain the following bound after one time integration over the window N + 1 for δ x f (1) = x f (1)x f (1) and δ x s (1) = x s (1) xs (1) :

δ x f (1) δ x f (0) (T N ) + O(H)( δ x f (0) + δ y f (1) ) δ x s (1) δ x s (0) (T N ) + O(H)( δ x s (0) + δ y s (1) ) (30)
Our goal is now to provide a bound to the difference in the primal unknown that will be imposed at the next iteration δ y s (1) . For a fixed time t, the algebraic equations are summarized to F(0) = F(1) = 0 with:

F(θ ) =    r f (x f
(1),θ , x s (0),θ , y f (1),θ , λ (1),θ ) r λ (x f (1),θ , x s (0),θ , y f (1),θ , y s (0),θ ) r s (x f (1),θ , x s (1),θ , y s (1),θ , λ (1),θ )

   , θ ∈ {0, 1} (31)
where x θ can be written as

x θ = (1 -θ )x + θ x.
The identity F(1) -F(0) = 1 0 F ′ (θ ) dθ = 0 further gives:

1 0 independent of θ     ∂ y f r f ∂ λ r f 0 ∂ y f r λ 0 0 0 ∂ λ r s ∂ y s r s     δ y f (1) δ λ (1) δ y s (1)   +   0 ∂ y s r λ 0   δ y s (0)   dθ + O(1) δ x s (0) + δ x f (1) + δ x s (1)
= 0

In the equation above the arguments x f (1) and x s (1) of the Jacobian that are placed in the neighborhood U γ 0 of size O(γ 0 ) of the solution (x f ⋆ , x s ⋆ ) are neglected. If γ 0 is sufficiently small, then we can easily solve the previous equation, provided the algebraic equations give a non-singular matrix 3 .

-

  ∂ y f r f ∂ λ r f 0 ∂ y f r λ 0 0 0 ∂ λ r s ∂ y s r s   -1   0 ∂ y s r λ 0   =     ∂ y s r s -1 ∂ λ r s ∂ y f r λ ∂ y f r f -1 ∂ λ r f -1 ∂ y s r λ     (32)
thus the contractivity coefficient α can be defined as:

α = max n ∂ y s r s -1 ∂ λ r s ∂ y f r λ ∂ y f r f -1 ∂ λ r f -1 ∂ y s r λ
Using this inverse calculation in [START_REF] Hughes | Implicit-explicit finite elements in nonlinear transient analysis[END_REF] gives the following bounds for the difference between the two coupled problems after one iteration:

δ y f (1) O(1) δ x f (0) + δ x s (0) + δ y s (0) + δ x f (1) + δ x s (1) δ y s (1) α δ y s (0) + O(1) δ x f (0) + δ x s (0) + δ x f (1) + δ x s (1)
where:

α = α + O(1) ∆ x f (0) + ∆ x s (0) + ∆ λ (0) + ∆ x f (0) + ∆ xs (0) + ∆ λ (0)
Inserting the inequality above in the Lipschitz condition in [START_REF] Hübner | A monolithic approach to fluid-structure interaction using space-time finite elements[END_REF] gives:

δ x f (1) O(H) δ x f (0) + δ x s (0) + δ y s (0) + δ x f (0) (T N ) δ x s (1)
O(H) δ x f (0) + δ x s (0) + δ y s

(0) + δ x s (0) (T N ) δ y s (1) ( α + O(H)) δ y f (0) + O(1) δ x f (0) + δ x s (0)
Rewriting the equation above in matrix form completes the proof of Lemma 1.

Lemma 2 Provided that the assumptions of the previous Lemma are satisfied and assume that α < 1 and C > α, one can write:

∃ Ĉ ∈ | + such that ∀k > 1, ∀H H 0   δ x f (k) δ x s (k) δ y s (k)     1 + ĈH δ x f (0) (T N ) 1 + ĈH δ x s (0) (T N ) Ĉ   +   Ĉ(4 Ĉ + 1)H μk-2 Ĉ(4 Ĉ + 1)H μk-2 4 ĈH μk-1 Ĉ(4 Ĉ + 1)H μk-2 Ĉ(4 Ĉ + 1)H μk-2 4 ĈH μk-1 4 Ĉ μk-1 4 Ĉ μk-1 μk + ( μ -α) k     δ x f (0) δ x s (0) δ y s (0)   with μ = α + 2CH α 2C + √ H Proof (of Lemma 2)
Although the successive values of (x f (k) , x s (k) , y s (k) ) remain in the neighborhood U γ 0 of the solution, Lemma 1 shows that iteration error is mainly governed by the matrix:

J =   CH CH CH CH CH CH C C α +CH   (33)
Recursive application of Lemma 1 leads to:

  δ x f (k) δ x s (k) δ y s (k)   J k   δ x f (0) δ x s (0) δ y s (0)   + k-1 ∑ i=0 J i   δ ∆ x f (0) (T N ) δ ∆ x s (0) (T N ) 0   (34)
Hence, the goal is to express in a relatively simple way the bounds for the elements of J k and for the first column of ∑ J i . The classical method used to obtain such bounds seeks to first transform the matrix into corresponding diagonal form J and then to bound the power of the eigenvalues. This leads to tedious but straightforward calculations; for more details the reader is invited to consider the Lemma 3.2. in [START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF].

Proof (of Theorem 1)

The errors on the interest variables like ε x f ,N on x f for the fluid part, ε x s ,N on x s for the structure part, and ε y s ,N on y s used for the coupling, are split in two terms: e .,N+1 representing error propagation from one window

N (t ∈ [T N-1 , T N ]) to the next N + 1 (t ∈ [T N , T N+1 ]); d .,N+1 corresponds to local error contribution on the window of interest N + 1.
The proof is organized as follows: applying Lemma 2 with suitable arguments, yields estimates for error propagation e (see proof (i)) and local error d (see proof (ii)). Then these estimates are combined to bounds set for the global errors ε (see proof(iii)). We can then show by induction that the global error bound is always verified (see proof (iv)). The constants µ and α are explicitly stated at the end of the proof of Theorem 1 (v).

(i) Estimate of propagation error contribution is the first part of the proof, where Lemma 2 with the following suitable arguments is applied:

  x f (0) xs (0) ỹs (0)   ←-    x f (k max ) N x s (k max ) N y s (k max ) N    and   x f (0) x s (0) y s (0)   ←-   x f ⋆ N x s ⋆ N y s ⋆ N   (35)
The symbol, '←-' indicates the initial guesses (obtained with zero-th, first or second order predictors) from one window N to the next one N + 1. For instance, in the simplest predictor is the constant function that leads to an error of size O(H):

(x f (0) N+1 , x s (0) N+1 , y s (0) N+1 ) = (x f (k max ) N , x s (k max ) N , 0)
The values (x f

(k max ) N , x s (k max ) N , y s (k max ) N
), in equation ( 35) above, obtained by numerical integration on the previous window N, can be considered as the best representation of the exact solution of the problem

(x f ⋆ N , x s ⋆ N , y s ⋆ N
) on this window. The choice of initial values [START_REF] Ibrahimbegovic | Energy conserving & decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations[END_REF] gives by definition the propagation error in the N + 1-th windows for the (k)-th iteration: δ x f (k) = e x f ,N+1 , δ x s (k) = e x s ,N+1 and δ y s (k) = e y s ,N+1 .

It is assumed that the chosen initial guess operator ←satisfies Lipschitz condition:

∃L * ∈ | such that    ∆ x f (0) L * ε x f ,N ∆ x s (0) L * ε x s ,N ∆ y s (0) L * ε y s ,N (36)
It will extrapolate x f and x s continuously from one window to another: δ x f (0) (T N ) ε x f ,N and δ x s (0) (T N ) ε x s ,N . Therefore, the application of Lemma 2 to the (k)th iteration on window N with α < 1 and a window size H small enough so that µ < 1 leads to:

  e x f ,N+1 e x s ,N+1 e y s ,N+1     1 +C * 1 H 1 +C * 1 H C * 1 H 1 +C * 1 H 1 +C * 1 H C * 1 H C * 1 C * 1 α *     ε x f ,N ε x s ,N ε y s ,N   ( 37
)
with

C * 1 ∈ | + and α * := L * ( μk + ( μ -α) k ).
(ii) Estimate of local error contribution consists again in applying Lemma 2 with suitable arguments:

  x f (0) xs (0) ỹs (0)   ←-   x f ⋆ N x s ⋆ N y s ⋆ N   and   x f (0) x s (0) y s (0)   =   x f ⋆ N+1 x s ⋆ N+1 y s ⋆ N+1   ( 38
)
Since the solution of the problem

(x f ⋆ n , x s ⋆ n , λ ⋆ n ) is a fixed-point of the iteration sequence, the local error contributions are measured by (δ x f (k) , δ x s (k) , δ y s (k) ) = (d x f ,N+1 , d x s ,N+1 , d λ ,N+1
). Furthermore, the use of the solution of the problem as the initialization of the iterative sequence yields (δ

x f (0) , δ x s (0) , δ y s (0) ) = (∆ x f (0) , ∆ x s (0) , ∆ y s (0) ), x f (0) (T N ) -x f (0) (T N ) = 0 and xs (0) (T N ) -x s (0) (T N ) = 0.
Thus the application of Lemma 2 then gives:

  d x f ,N+1 d x s ,N+1 d y s ,N+1   μk-2 ∆ x f (0) + ∆ x s (0) + μk-1 ∆ λ (0)   C * 2 H C * 2 H C * 2   (39)
with a positive constant C * 2 . (iii) these estimates are combined with the bounds set for the global errors ε N+1 e N+1 + d N+1 :

  ε x f ,N+1 ε x s ,N+1 ε y s ,N+1     1+C * 1 H 1+C * 1 H C * 1 H 1+C * 1 H 1+C 1 H * C * 1 H C * 1 C * 1 α *     ε x f ,N ε x s ,N ε y s ,N   + μk-2 ∆ x f (0) + ∆ x s (0) + μk-1 ∆ y s (0)   C * 2 H C * 2 H C * 2   (40)
If the contraction condition α * < 1 is fulfilled, the behavior of such coupled error recursions is known (see [START_REF] Deuflhard | One-step and extrapolation methods for differential-algebraic systems[END_REF], Lemma 2) and can be written as follows:

max N ε x f ,N + ε x s ,N + ε y s ,N C • max N μk-2 ∆ x f (0) + ∆ x s (0) + μk-1 ∆ y s (0) (41)
with a positive constant C.

As the total error is the expressed by the sum ε x f ,N + ε x s ,N + ε y s ,N , the last expression completes the proof of Theorem 1. (iv) As the initial guess operator ←leads to initial errors ∆ x f (0) , ∆ x s (0) and ∆ y s

(0)
of size O(H), the right-hand side of the last equation remains bounded for all H < H 0 , if H 0 is sufficiently small:

max N ε x f ,N + ε x s ,N + ε y s ,N C • max N μk-2 ∆ x f (0) + ∆ x s (0) + μk-1 ∆ λ (0)
γ 0 [START_REF] Legay | An eulerian-lagrangian method for fluid-structure interaction based on level sets[END_REF] This further shows that errors is bounded by γ 0 , and that approximate numerical solution remains in the neighborhood U γ 0 of the solution. (v) The constants α and μ which appear in proofs of Lemma 1 and Lemma 2 are given as follows:

α = α + O(1) max n ( ε x f ,N + ε x s ,N + ε λ ,N ) + O(H) = α + O(H) and μ = α + O(H)
These results give the following order of magnitude for α * :

α * = L * ( μk + ( μ -α) k ) = L * ((α + O(H) k + O(H k )) (43)
The expression above confirms the criterion of stable error propagation for α * < 1. For the iterative DFMT-BGS with an appropriate guess operator, the stability is guaranteed when α < 1 and when the window size is small enough.

Validation and numerical examples

Lid-driven cavity flow with flexible bottom

Having recognized the value of this problem for establishing the standard benchmark in fluid-structure interaction, we will provide a sufficiently detailed presentation of the results obtained in order to promote the comparison. The lid-driven cavity flow considers 2D fluid flow problem in a square domain. The imposed boundary conditions at three out of four sides are the zero value of velocity. Only at the top of the cavity we impose the velocity with a given nonzero value. The geometry and boundary conditions are depicted in further details in Figure 3.

For the case of fluid-structure interaction problem [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Gerbeau | A quasi-Newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows[END_REF] the flexible membrane is placed at the bottom boundary (see Figure 3(b)). The flow is governed with Navier-Stokes equations and we consider a Saint-Venant Kirchoff material able to undergo finite deformation for the structure part. The material properties are chosen as follows: the fluid density is ρ f = 1kg • m -3 , the kinematic velocity ν f = 0.01m • s -2 , the structure density is ρ s = 500kg • m -3 , the Young modulus E = 250Pa and the Poisson ratio ν = 0.

The flow has a dominant convective term, which is usually the most difficult to solve accurately. Furthermore, when the imposed velocity is sufficiently small, the flow is laminar and incompressible. In fact, this lid-driven cavity flow (with a rigid bottom) problem is traditionally used as a validating example for CFD codes. Therefore extensive literature is v(y = 1m,t) (a) Lid-driven fluid flow in cavity v(y = 1m,t) (b) Modified lid-driven cavity for FSI Fig. 3 The lid-driven cavity example dedicated to its study, ever since the early works in computational fluid dynamics [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF] to more recent reviews [START_REF] Bruneau | The 2D lid-driven cavity problem revisited[END_REF] using different fluid solvers [START_REF] Hortmann | Finite volume multigrid prediction of laminar natural convection: Bench-mark solutions[END_REF][START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. It is also proposed as a test case for the fluid solver of OpenFOAM [START_REF] Opencfd | Openfoam home page[END_REF] we employ in this work. In the present test, the cavity is discretized by finite volume mesh of 32 × 32 cells.

For the lid-driven cavity flow in Figure 3(a), all boundary nodes are constrained at the Dirichlet boundary with imposed velocity, and pressure field remains undetermined up to a constant. If this does not lead to any problem when only fluid flow is considered with no interaction with the structure, a special care has to be taken for the fluid-structure interaction case where the imposed velocity condition should satisfy incompressibility condition and where the exact value of pressure is needed to define structure boundary conditions. This difficulty was studied in [START_REF] Küttler | A solution for the incompressibility dilemma in partitioned fluid-structure interaction with pure Dirichlet fluid domains[END_REF]. One way to overcome this difficulty for fluid-structure interaction lid-driven cavity problem, as proposed in [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Küttler | A solution for the incompressibility dilemma in partitioned fluid-structure interaction with pure Dirichlet fluid domains[END_REF] is to unconstrain two nodes on each size of the cavity (see Figure 3(b)). One subtle point concerns the chosen discretization techniques: for stabilized FE approximation used in [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Küttler | A solution for the incompressibility dilemma in partitioned fluid-structure interaction with pure Dirichlet fluid domains[END_REF] constraints are removed at the nodes, whereas for the FV method used herein the corresponding constraints at the boundary pertain to face of cells. Another modification of lid-driven cavity problem, besides taking into account the fluid-structure interaction, concerns the time-dependent velocity boundary condition defined as:

v • e x = 1 -cos 2π t T char ( 44
)
where T char = 5s. For such harmonic function, the solution of the fluid flow within the fluidstructure interaction problem exhibits an oscillating behavior that is reached after a short transition period. The maximum value of the Reynolds number in the cavity reaches Re = 200, and thus the flow can be considered as laminar. This example represents a very good benchmark test for the case when the flow is mainly driven by incompressibility. Moreover, it allows to obtain pretty good results with a rather coarse mesh; namely with the Reynolds number not bigger than Re = 200, a mesh with 32x32 cells and a second order FV solver is sufficient to get an accurate solution [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF][START_REF] Hortmann | Finite volume multigrid prediction of laminar natural convection: Bench-mark solutions[END_REF]. Due to such a small mesh size, computations are fairly inexpensive, especially for a fluidstructure interaction problem (for each time step of ∆t = 0.1s, the flow computation takes 1.08 × 10 -1 s and the structure motion takes 2.95 × 10 -3 s when run on a single processor). The problem quickly reaches a harmonic steady-state solution, which provides a perfect platform to test the energy creation or dissipation for partitioned solution to fluid-structure interaction problem.

For the fluid flow only, we can stick with the problem where boundary conditions are imposed on the velocity field, with no condition on the pressure field p. These boundary conditions are also called zero gradient (or Neumann zero). In this case any constant can be added to the computed pressure field, with no change for the fluid part since the gradient of the pressure field in (9) will filter out any constant pressure. However, this can lead to a non-unique pressure, and most of the iterative solver will fail for such a problem if no remedy is proposed. This is even a bigger problem for fluid-structure interaction, where the exact pressure is needed to obtain the corresponding structural motion.

If the fluid flow case only is considered there are two possible modifications (see [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]) of this problem in order to specify a unique pressure field for Neumann zero boundary. Namely, we can impose either the fixed value (e.g. zero) of the pressure in a arbitrary cell, or we can consider the pressure field with the mean value equal to an imposed constant.

However, for the present version of the lid-driven cavity that accounts for fluid-structure interaction problem, the issue of pressure computation becomes more demanding since not only the gradient of the pressure field, but also its actual values must be known at the fluid-structure interface. Unfortunately, a number of previous works (see [START_REF] Wall | Fluid-structure interaction based upon a stabilized (ALE) finite element method[END_REF][START_REF] Küttler | A solution for the incompressibility dilemma in partitioned fluid-structure interaction with pure Dirichlet fluid domains[END_REF][START_REF] Gerbeau | A quasi-Newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows[END_REF]) did not clearly specify the imposed condition on the pressure field. For example in [START_REF] Bathe | A mesh adaptivity procedure for CFD and fluid-structure interactions[END_REF], where the pressure field is imposed at the flow inlet and outlet, the obtained result is totally different from the one given by the previous works predicting a high pressure at the bottom of the cavity leading to an average negative displacement of the structure in the e y direction, and requiring the use of a stiffer material in order to remain in the acceptable displacement range.

In order to illustrate the importance of the chosen pressure condition in fluid-structure interaction analysis we take the case of an almost rigid structure at the bottom and compute the total force (from pressure and viscosity) applied to its center (nodal force at (0.5, 0.0)). These results are presented in Figure 4, showing very different time evolutions of pressure force exerted on the structure for different choices of pressure boundary conditions. For further computations, the pressure was set to zero at the flow outlet. Other boundaries are considered to be Neumann zero for the pressure. The fluid domain is discretized with fourth order FV in space, and the time integration is carried out by implicit Euler scheme. The momentum equation is solved by PBiCG, with a DILU preconditioning whereas the pressure correction problem is worked out by iterative PCG, with DIC preconditioning. Two iterations of the PISO correction algorithm are performed at each time step, and the precision required for the iterative algorithms for pressure and velocity is 10 -8 . As the mesh is initially orthogonal, there is no reason to require initially non-orthogonal correctors. At a later stage the deformation of the bottom will eventually produce non-orthogonal meshes and we specify that two non-orthogonal corrections should be performed at each time step.

The fluid domain is subject to structure motion as the bottom of the domain. We assume that no cavitation takes place and that fluid follows the deformation of the structure. The fluid mesh deformation is handled by a smoothing process, based upon the solution to Laplacian equation with the coefficient that depends on the distance to the bottom. On vertical walls, the points are allowed only to move vertically, whereas elsewhere the boundary nodes are fixed.

For the structure FE model a spatial discretization with 16×Q8 (quadratic) elements is chosen. The time integration is carried out by implicit generalized HHT α-scheme in Note that the inertia of the flexible bottom leads to a positive mean value of structure displacement. Moreover, this further induces a large decrease in pressure and accordingly in the expected force amplitude. A tolerance TOL = 10 -7 is set for testing the convergence of the DFMT-BGS algorithm. With a fixed under-relaxation of ω = 0.25, a constant decrease of the residual with a low order is observed (see Figure 10). Aitken's relaxation allows to reduce the mean number of iterations required from 30 to 17 (see Figure 8). Aitken's relaxation allows also to improve the order of convergence. The convergence exhibits a less smooth behavior, with a faster decrease of the residual with the increase of relaxation parameters (for instance the 4 th iteration in Figure 12(a) and 10(a)). The relaxation parameter value is given in Figure 12 for two chosen times 39s and 41s used in Figure 10 to represent the convergence of the residual; The characteristic oscillations depicted in [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF] are also observed. The use of predictors does not change the order of convergence, but reduces the initial residual (see Figure 11). For example the characteristic number of iterations to reach required precision decreases from 17 iterations for a zero order to 7 for a second order predictor (Figure 9).

At the end of this section, we also present the results obtained for explicit DFMT coupling algorithms applied to the lid-driven cavity with flexible structure at the bottom. The added mass effects characteristic of explicit algorithms for the case of an incompressible flow is observed in Figure 13, causing each computation to diverge sooner or later.

Our result divergence was not in agreement with the results convergence presented in [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF] for the least accurate predictors. The main reason, we believe, is due to the way the enforcement of the true incompressibility condition for the fluid problem is really enforced. Namely, with FV and PISO algorithm used herein, we can strongly enforce incompressibility condition, which is less the case when using stabilized FE discretization as in [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF].

The choices recognized to be important for triggering earlier instability of explicit coupling are: increasing the order of predictor (see Figure 13(a)), decreasing the window size for synchronization, (see Figure 13 Fig. 12 The lid-driven: evolution of relaxation parameter for two characteristic time steps ratio between fluid and solid with larger ρ f /ρ s , the instability occurs earlier; larger fluid velocities and more flexible structure trigger instability more quickly. This is in fully agreement with the theoretical and numerical observations made in [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF][START_REF] Förster | On the geometric conservation law in transient ow calculations on deforming domains[END_REF]. The general conclusion is that accurate time integration of the sub-problem and coupling algorithms are more sensitive to the artificial "Added-Mass Effect" and diverge more easily.

Flexible appendix in a flow

This example first proposed in [START_REF] Wall | Fluid-structure interaction based upon a stabilized (ALE) finite element method[END_REF], and subsequently also studied in [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF][START_REF] Hübner | A monolithic approach to fluid-structure interaction using space-time finite elements[END_REF][START_REF] Hermann | Algorithms for strong coupling procedures[END_REF] has also been used as a benchmark for fluid-structure interaction. It is close to the benchmark proposed in [START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF] itself constructed from the traditional CFD benchmark proposed in [START_REF] Schäfer | Benchmark computations of laminar flow around a cylinder[END_REF]. We consider a fixed square bluff body, with a flexible appendix attached to it, immersed within an incompressible flow (see Figure 14) filling the whole domain. At a sufficiently long distance from this body, the flow is uniform with an imposed velocity v; the corresponding value of the Reynolds number with respect to the characteristic size of the obstacle is Re = 330. For this Reynolds number, the flow exhibits a transient behavior with vortices separating from the corner of the square. These vortices induce alternative drop and increase in the pressure field behind the rigid bluff body at a frequency that depends on the Reynolds number and the shape of the bluff body [START_REF] Roshko | Of the Development of Turbulent Wakes from Vortex Streets[END_REF]. The vortex shedding induces oscillations of the flexible appendix.

The thickness of the appendix as well as the material properties are chosen so that its first eigen-frequency is close to the frequency of the vortex shedding. The material parameters are respectively ρ f = 1.18×10 -3 kg.m -3 for the density and ν f = 1.54×10 -1 m 2 .s -1 for the for the fluid (air at 20 • C). The imposed velocity at the left hand side is v = (51.3m.s -1 , 0).

For the solid part, consider a density ρ s = 0.1kg.m -3 , Young modulus is of E s = 2.5 × 10 6 N.m -2 and the Poisson ratio coefficient of ν s = 0.35. The first eigen-frequency of the

: v • n = 0 outflow p = 0 ρ s , E s , ν s v = v f ρ f , ν f slip: v • n = 0
Fig. 14 The benchmark used for FSI problems

The flexible appendix is discretized with 20 nine-node elements, with quadratic polynomes description of large structural displacement. Neo-Hookean and Saint-Venant-Kirschoff materials are used. The time discretization is carried out by a generalized HHT-α scheme with the following parameters: ρ ∞ = 1 2 ; β = 4 9 ; γ = 5 3 ; and α = 2 3 At each iteration, the linear system is solved by the direct solver for asymmetric matrices.

The interface matching computation is carried out by the proposed DFMT-BGS solver, that can easily converge since the channel is open and the flow is not mainly driven by incompressibility. The Aitken relaxation technique is used with an initial value of 0.5, that rapidly increases to 1. No more than 4 iterations are required to reach the required tolerance on the interface displacement residual that is set to 1 × 10 -7 . The results for pressure and velocity field at several instants are given in Figure 15. The deformed shapes of the appendix also given in Figure 15, reveal the oscillations dominated by the first mode.

The displacement at the free-end of the appendix is plotted in Figure 16 for both Saint-Venant-Kirschoff and Neo-Hookean solid materials. The two results are very close, since the appendix deformations remain small, despite its large displacement and rotations. The long term response (see Figure 17) indicates an almost harmonic response dominated by first eigen-frequency of the structure. In Figure 17 we also give a comparison with the results from the literature in term of the maximum amplitude of motion. Despite a wellknown sensitivity of the computed result with respect to the initial condition [START_REF] Hübner | A monolithic approach to fluid-structure interaction using space-time finite elements[END_REF], we get the answers obtained very close to the previous results from the literature based upon a FE discretization for both fluid and solid parts and obtained either by a monolithic [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF] or a partitioned approach [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF].

Contrary to the lid-driven cavity with a flexible bottom, the small number of iterations required to solve the fluid-structure interaction problem of the oscillating appendix suggests that an explicit coupling can also be used for solving this problem. The results from explicit DFMT algorithm presented in Figure 18 for the free-end displacement compare to a reference solution obtained with an implicit computation.

Using better predictors is supposed to reduce the errors made in term of residual and energy. In Figure 19 coupling: long term response with implicit coupling. Maximum amplitude comparison with [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF] (dotted line), [START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF] (solid line), [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF] (dashed line). second order predictor. All the results confirm the trends we expected, with a decrease of errors when the predictor order increases. The final study is then considered with respect to the size of time steps. In Fig 20, the maximum residual error on the time interval t ∈ [0, 15s] is presented as a function of the time step size. The error is observed to decrease with a decreasing time step size. However, when the time steps become too small, the added mass effect triggers the divergence of the computation. Thus, only the less sensitive schemes with a zero order predictor are able to solve the coupled problem with the smallest time step.

Remark:

The extension to three-dimensional cases does not require the introduction of theory, but requires tools to solve efficiently the problem with a large number of d-o-f, with the use of paralleling for instance. A generalization of the two dimensional case presented

Conclusion

In this work we examine partitioned solution approach for nonlinear fluid-structure interaction problems. The partitioned approach is preferred for its modularity and the possibility of re-using existing software developed for each subproblem (see Part II). The partitioned approach used here is based on the Direct Force-Motion Transfer. Both explicit and implicit coupling algorithms for multi-physics problems are detailed. An explicit strategy leads to the so-called "added mass effect", and for that justifies the use of more costly implicit solvers for the case of incompressible fluid flows.

In this work, the problem of enforcing the fluid-structure interface matching is handled by the fixed-point strategy (DFMT-BGS) with an adaptive relaxation parameter. This strategy shows a sufficiently robust performance, especially for the example where the flow is not highly constrained by incompressibility. In fact, we showed by direct proof the stability of the implicit DFMT-BGS algorithm which is valid for the fully nonlinear fluid-structure interaction problem.

Algorithm 2 1 s

 21 Poincaré-Steklov operator for structure: S -Require: Solid state variable at time T N 1: Impose boundary traction force at fluid-structure interface: σ s n = λ on Γ × [T N , T N+1 ] 2: Solve structure problem: R s (u s ; λ ) = 0 on Ω s × [T N , T N+1 ] 3: Get boundary displacement at the interface: u = u s on Γ × [T N , T N+1 ]

0 Fig. 4

 04 Fig.4Force acting at the center of the flexible bottom for some fixing pressure strategies

Fig. 5

 5 Fig. 5 Lid driven cavity with a flexible bottom: snapshots with pressure field and streamlines for different time steps

Fig. 8

 8 Fig. 8 Number of iterations per time step for the BGS solver with fixed and Aitken's relaxation (order 0 predictor)

Fig. 9

 9 Fig.9Number of iterations per time step for the BGS solver with predictor of order 0, 1 and 2 (Aitken's relaxation)

Fig. 10 2 )Fig. 11

 10211 Fig. 10 Residual convergence for fixed and Aitken's relaxation with order 0 predictor

25 (

 25 (b)), increasing the order of fluid integrator, (see Figure 13(c)), use of non-collocated algorithm like the DFMT-ISS (see Figure 13(d)and predictor O(∆ 1) Aitken and predictor O(∆t) Aitken and predictor O(∆t 2 ) Fixed relaxation ω = 0.and predictor O(∆ 1) Aitken and predictor O(∆t) Aitken and predictor O(∆t 2 ) Fixed relaxation ω = 0.25 (b) Time t = 41.0s

Fig. 15 16Fig. 17

 1517 Fig. 15 Oscillating appendix in flow: velocity and pressure field snapshots

  Fig.17Implicit coupling: long term response with implicit coupling. Maximum amplitude comparison with[START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF] (dotted line),[START_REF] Matthies | Partitioned strong coupling algorithms for fluid-structure interaction[END_REF] (solid line),[START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF] (dashed line).

Fig. 18

 18 Fig.18Comparison between explicit and implicit coupling algorithm based computations for the displacement of the appendix.

  Fig. 19 Energy error at the interface for different explicit coupling schemes.

Fig. 20 Fig. 21

 2021 Fig. 20 Maximum residual error for explicit coupling schemes with different time step sizes and predictors.

It is also assumed that fluid flow and structure motion are computed by the corresponding software products, to be coupled for solving a fluid-structure interaction problem; the details of software coupling are discussed in Part II of this paper.

Regarding this regularity assumption, see the remark at the end of this proof.

This requires non-singular Jacobians, see[START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF][START_REF] Brenan | Numerical solution of initial-value problems in Differential-Algebraic Equations[END_REF] for more details

u y (0.50,0) u y (0.25,0) Fig. 6 Displacements at the fluid-structure interaction interface at the middle and the quarter lengths: u(0.5, 0) and u(0.25, 0). order to maximize numerical damping [START_REF] Förster | Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows[END_REF]. The non-linear algebraic equations is solved at each time step by the Newton iterative algorithm with a prescribed tolerance of 10 -8 . A non-symmetrical direct solver is used at each iteration.

The typical results for the pressure field, streamlines and structure deformation are plotted in Figure 5. In Fig 6 the time history of the displacement is displayed when no interaction is considered. Namely, the forces are computed by the CFD-based component, and applied to the structure, but the structure motion is not imposed on the fluid domain, and thus there is no need to compute any mesh motion. The same kind of results are presented in Figure 7 for implicit coupling computations with a window size of ∆t = 0.1s. The results largely differ from the one obtained with a weak interaction model. The displacement remains exactly the same for all predictors and relaxation techniques, once convergence to a residual norm less than 10 -7 is obtained. The results are close to those obtained in [START_REF] Wall | Partitioned analysis approach of the transient coupled response of viscous fluids and flexible structures[END_REF][START_REF] Gerbeau | A quasi-Newton algorithm based on a reduced model for fluid-structure interaction problems in blood flows[END_REF][START_REF] Küttler | A solution for the incompressibility dilemma in partitioned fluid-structure interaction with pure Dirichlet fluid domains[END_REF] (see Figure 7 f), which is quite sufficient to get an accurate representation of the flow, and also of the fluid loading on the structure. The model accuracy is comparable to the one used in [START_REF] Dettmer | A fully implicit computational strategy for strongly coupled fluid-solid interaction[END_REF], built with 4300 finite elements. The PISO algorithm and the Euler implicit scheme with a time step ∆t = 0.004s are used. The PCG solver is used for pressure correction step and mesh motion equations and PBiCG for the momentum predictor.


Nancy Honicker nancy.honicker@univ-paris8.fr Taking the Veil (an essay based on interviews with veiled Muslim students at the author teaches) At the heart of Christianity are the words spoken by Jesus to his disciples in an upper room in Jerusalem, a few days before his death: I am the way, the truth, and the life: no man cometh unto the Father, but by me. For the Muslim, there is but one God, Allah, and Mohammed is his prophet.

Any mingling of the divine and the human -the Word made flesh -is sin. Can these two mutually exclusive positions be reconciled in a shared love of God and humankind? I admit to not having the answer in regards to doctrine, but living and working in a community with a large Muslim population, I have had the opportunity to observe Islam in action and to get to know those who practice it in their daily lives. Rather than enter into a debate about what may be irreconcilable differences, I prefer to look at how Muslims live their faith. The context, in this case, is French. I am an associate professor in a state-run French university, the University of Paris 8, located a few miles north of the center of Paris, in the town of Saint-Denis, at the very heart of those dangerous suburbs where violence erupted in November 2005. Until 1789, Saint-Denis, with its magnificent basilica, was the seat of the spiritual authority of the French monarchy. Today its churches are almost empty and the number of conversions to Islam is on the rise.

In a narrow street a few steps from the basilica, the Tawhid Cultural Center opens its doors to all those Muslims who seek to deepen their faith through greater separation from secular society. At all hours of the day, clusters of men gather on the sidewalk in front of the center. When it is the hour of one of the five daily Muslim prayers, the street is crowded, and often traffic is blocked as men and women flow in and out. Some of the men, though not a majority, wear the kamis, a tunic that descends to mid-calf, in imitation of the dress of the Prophet. The women are strictly veiled, with a white band covering that part of their forehead left uncovered by the veil. Though I have never visited the center, when I pass by, I get a sense of proselytising verve. Yet, movements such as "Tawhid," which means the oneness of God, remain on the fringes when it comes to bringing new converts to Islam.

To understand the appeal of Islam, however, one must first understand that today in France, religious practice in general -and that of Catholicism in particular -is considered "ringard" -a word which is difficult to translate, for it means both old-fashioned and stale. Though France has an excellent Catholic daily, La Croix, and a very active and innovative Catholic press, religious practice is for the most part desultory. A majority of the French may be nominally Catholic, few go to church, and fewer still, especially among the younger generations, understand the fundamentals of their faith. I recently took a survey among my students to ask how many knew the significance of Ascension Day and Pentecost, both public holidays in France. Not a one had a clue.

When I began teaching at the University of Paris 8 in 1991, there was not a single veiled students on campus. Today I have at least one veiled student in every class; in some, a majority of students are Muslim, which does not mean most women wear the veil. Many of these students are foreign, from Algeria or Morocco, but many are also French, and some are converts to Islam. To understand the appeal of this religion, I decided to turn to the students themselves, beginning with Laetitia, a recent Catholic convert to Islam. And I will begin by stating that Laetitia's conversion does not surprise me. She was born in the town of Saint-Denis, with its large Muslim population. She was also born into one of those nominally Catholic families, who never set foot in church. Despite the absence of any religious education, early on Laetitia was stirred by yearnings to learn more about God. As her family had placed her in a Catholic primary school -more for reasons of discipline and security than of religion -she was able to turn to the nuns to learn more about the Catholic faith. All her decisionsto attend catechism, to be baptized at age eleven, and to take first communion -were made on her own. Once she entered junior high school, many of her classmates were Muslim, and among them she met a girl who became her best friend. Though they did not practice the same religion, they shared the same love of God, which became the basis of a friendship that continues strong today.

Within her family Laetitia had been alone with her faith in God. Through Islam she came into contact with a shared discipline, a very important word in the context of French, and much of Western, society today: there were the fixed hours for prayer, the same prayers recited by millions all over the world, the same positions, standing or kneeling, in imitation of the Prophet, the same month-long fast, followed each evening by festive meals with family and friends. This is not to say that such a sense of shared community is absent from Christianity, but certainly it is absent from the practice of most French Catholics today. In many ways, Islam, the second religion of France, presents itself as a welcoming faith with the power to bring the believer closer to God through the structure religion brings to daily life. Laetitia insists on this dimension in her practice of the five daily prayers, a discipline which has become a gentle prodding, a reminder to take time out to think of God.

After much soul searching, Laetitia converted to Islam in March 2004. Like all Muslims, she accepts Jesus as a prophet, but today she sees Christianity as incomplete and rejects the divine nature of Christ. She would like to wear the veil as an outward sign of her religious convictions but refrains, in part for fear of hurting her mother, who has not easily accepted her daughter's choice.

While Laetitia was studying abroad in spring of 2005, her father died. Today her mother fears she may have lost a daughter as well. To reassure her, Laetitia has chosen not to take the veil.

Living daily among women wearing the veil, I observe that it comes in many shapes and sizes, many colors and textures, and many different styles. In the street where I live, I often pass a woman covered in black from head to toe, her face covered with a veil of black gauze. At the university, however, the style that dominates is one of "spiritual elegance," which has taught me to recognize the beauty of the veil and the inner beauty of the students who wear it. While interviewing Muslim students to prepare this article, I met Fatima, a thirty-year-old Moroccan student from Agadir, born into a practicing Muslim family, who has been studying in France since 2002. When she was sixteen, she was required to cover her head for medical reasons. After a few months, when it was safe to remove the scarf, she had grown so used it that she decided to keep it on: the veil had become the reminder of God's presence in her life. Sitting across from Fatima, looking into her thoughtful eyes, I feel that presence as well.

The day we meet her head is covered by a silk scarf the color of fine French champagne.

The entire oval of her face is visible, her forehead, her dark eyes, a long thin nose and full, finely sculpted lips. The veil covers her ears and neck and is expertly fastened beneath the chin. She is also wearing a long-sleeved tunic of light beige cotton and matching slacks. She is, to my eyes, the epitome of "spiritual elegance", a young woman at peace with herself, expressing inner and outer harmony, a force for good in the world. In her presence one forgets the divisions of Christian, Muslim or Jew.

Yet no religion is free of fear, and like the fire-and-brimstone brand of Christianity, Islam too knows how to instil fear in its believers. One of my students told me it was a sin not to wear the veil and the unveiled woman puts herself in danger of hellfire. In certain contexts the veil can also participate in the ongoing "culture wars." Miniskirts and low-cut tops become "green lights," and women who choose to wear them are declaring themselves "for sale." Looking around at her scantily clad classmates, the woman wearing the veil may find it all to easy to be tempted by a sense of moral superiority. Even without going so far, many Muslim women are made uncomfortable by certain aspects of French society which they deem superficial and demeaning to women.

Walking in the streets of Paris, where naked women meet and greet me from billboards all over the city, buying my newspaper at a kiosk where pornography takes precedence over the daily news, I can understand their point of view. The French may consider me puritanical and prudish, but at times I am literally hurt -I feel a stab of pain to the heart -by the use made of women's bodies in advertising in France, hurt as well by pornography on display to adults and children alike. This is freedom of expression, one of the pillars of democracy, yet I ask myself: would these images be so blatantly displayed in a truly democratic society? Thus, for some young Muslim women, the veil becomes protection and a sign of sexual purity. Laetitia told me she hopes to know and love one man alone. Though she does not wear the veil, her modest attire and her embrace of Islam have sustained her in her choice. Naima, another student, adds that Islam considers men the weaker sex and by choosing the veil, women help protect them from themselves. For Fatima, the veil is a "positive obstacle": the respect of a Muslim man for a veiled woman includes recognition of his weakness.

On one point, all agree: the presence of the veil in the life of a Muslim woman brings her closer to God.

In her 1995 memoir Dreams of Trespass, Fatima Mernissi describes the life of a young girl born in 1940, in Fez, Morocco. Her father, a powerful patriarch, had the means to keep a harem, where his wives and daughters, though privileged and wealthy, lived among themselves, deprived of instruction and contact with the outside world. This was the traditional world of Islam, where the public sphere was reserved for men, and women spent their lives behind high walls and locked doors. Fatima Mernissi herself could have grown up illiterate had not the madrassahs of Morocco been opened to women in 1943, during the period of the country's struggle for independence from France. Today she is a world-renown writer and sociologist, having published numerous books on the relation of women to Islam. * To her, no development in the past century has been more important than the entry of Muslim women into the public sphere. And we are far from seeing the end of the changes this movement has brought about. In this context, the veil can become a symbol of liberation, and women may very well be the future of Islam.

Finally

  , to return to doctrine, if I were to sum up the significance of the Islamic veil to the students with whom I have spoken, I would say it is the Muslim woman's answer to Saint Paul's exhortation to "pray without ceasing." My students and I may have chosen different paths, yet we are all seeking -as best we know how -the same God of love. * such as The Forgotten Queens of Islam, Beyond the Veil, and Islam and Democracy: Fear of the Modern World.