An introduction to the simulation of fluid and structure dynamics - PowerPoint PPT Presentation

Turbulence The Navier-Stokes equations become unstable when the transport term ( u · ∇ u ) is dominant with respect to the viscous term (2 µ D ( u )), i.e. when Re ≫ 1. In such a situation, energy is transferred from large eddies (large vortex structures) to smaller ones up to a characteristic scale, so called Kolmogorov scale where they are dissipated by viscous forces. l K ∼ L ∗ Re − 3 / 4 Kolmogorov scale A direct numerical simulation (DNS) aims at simulating all relevant scales up to the Kolmogorov scale. Therefore, the mesh size has to be h ≈ l K . � L � 3 ∼ Re 9 / 4 3D simulations # Dofs ∼ h Unfeasible even for nowadays computers for Re ∼ 10 4 ÷ 10 6 .

Kinetic energy decay 2 ̺ f | u | 2 and It is instructive to look at the fluid kinetic energy E = 1 perform a Fourier analysis with respect to space wavenumbers: � � R 3 E ( x ) e i x · k d x dA , ˆ E ( k ) = with S ( k ) = {| k | = k } S ( k ) log E(k) inertial region E ( k ) = c 2 k ε / 3 k − 5 / 3 production dissipation tranfer of energy of energy of energy L −1 −1 log(k) l k (integral scale) (Kolmogorov scale) E ( k ) ∼ ǫ 2 / 3 k − 5 / 3 is one of the main results of Kolmogorov The decay ˆ theory of turbulence, valid for homogeneous isotropic turbulence

Turbulence and boundary layer The golf ball The roughness of a golf ball energizes the boundary layer retarding flow separation. The Pressure induced drag is then reduced. Even if viscous drag is likely to be greater than that on a smooth ball the overall resistance is diminished.

Turbulence and boundary layer Laminar turbulent transition As the boundary layer grows even an initially laminar flow (i.e. with no turbulence) will become turbulent. The point where this change of behavior occurs is called laminar-turbulent transition point. The location of the Laminar-turbulent transition is very important to determine the drag of submerged bodies like a yacht bulb, since it affects the resistance. Mathematical models of the transition are still a partially open problem.

The level of simulation of turbulent flows ◮ Direct NS simulation Used only for special studies. In general is unfeasible ◮ Large Eddy Simulation We model only the smallest scale turbulence. Rather costly for high Reynolds. ◮ Differential Models. Based on the Raynolds averaging paradigm where turbulence is accounted for by modifying the viscosity µ → µ + µ t . We add advection-diffusion-reaction equations for quantities like the turbulence kinetic energy ( κ ) and dissipation ( ǫ ) and µ t = µ t ( κ, ǫ ). They need several empirical closure relation. Most used in engineering applications. Often indicated as RANS (Reynolds Averaged Navier Stokes). ◮ Algebraic models Also based on RANS. We use algebraic relations to link µ t to the mean flow velocity. Applicable only to attached flow and special geometries.

Numerical approximation of the Navier-Stokes equations

Finite Volumes They are based on the equations written in conservative form ∂ u ∂ t + div F ( u , p ) = S The domain is subdivided into a mesh of polygons (finite volumes, also called control volume) C i , i = 1 , ... N . By applying the divergence theorem � � � d u ( t ) + F ( u ( t ) , p ( t )) · n d γ = S ( t ) dt C i ∂ C i C i

Finite Volumes (cont) We now approximate u with its mean value on each control volume and the fluxes at the boundary of the volume are approximated with suitable averages between values at adjacent volumes (or using the boundary conditions) � d dt | C i | u i + F h ( u i , u i , f , p i , p i , f ) · n f = | C i | S i ( t ) f ∈ ∂ C i C i C il l l C i C il

Finite element approximation of the Navier-Stokes equations We will consider only flows in laminar regimes (no turbulence models) Navier-Stokes equations  ∂ u  ∂ t + ̺ f u · ∇ u − div σ f ( u , p ) = f f ,  ̺ f in Ω , t > 0 ,       div u = 0 , in Ω , t > 0 , u = g , on Γ D , t > 0 ,     σ f · n = d , on Γ N , t > 0 ,     u = u 0 , in Ω , t = 0 with σ f ( u , p ) = 2 µ D ( u ) − p I

Weak formulation Let us define the following functional spaces: Γ D (Ω)] d ≡ { v ∈ V , v = 0 on Γ D } V = [ H 1 (Ω)] d V 0 = [ H 1 Q = L 2 (Ω) Remark: if Γ D = ∂ Ω, the pressure is defined only up to a constant. The pressure space should be Q = L 2 (Ω) \ R . We multiply the continuity equation by a test function q ∈ Q and the momentum equation by a test function v ∈ V 0 and integrate over the domain. Integration by parts of the stress term: � � � − div σ f · v = σ f : ∇ v − ( σ f · n ) · v Ω Ω ∂ Ω ✯ 0 � � � � ✟ ✟✟✟✟✟✟ 2 µ D ( u ) : ∇ v − p I : ∇ v − ( σ f · n ) · v − ( σ f · n ) · v = � �� � � �� � Ω Ω Γ D Γ N = p div v = d

Weak formulation Navier Stokes equations – weak formulation Find ( u ( t ) , p ( t )) ∈ V × Q , u (0) = u 0 , u ( t ) = g ( t ) on Γ D such that � � ∂ u � � � � � ∂ t + u · ∇ u · v + 2 µ D ( u ) : ∇ v − f · v + d · v ̺ f p div v = Ω Ω Ω Ω Γ N � q div u = 0 Ω ∀ ( v , q ) ∈ V 0 × Q Weak solutions exist for all time ([Leray ’34], [Hopf ’51]). Uniqueness is an open issue in 3D.

Finite element approximation Introduce finite dimensional spaces of finite element type V h ⊂ V ; V h 0 = V h ∩ V 0 ; Q h ⊂ Q Observe that pressure functions need not be continuous. Moreover let us denote by u h 0 ∈ V h and g h ∈ V h (Γ D ) suitable approximations of the initial datum and Dirichlet boundary datum Finite Element approximation (continuous in time) Find ( u h ( t ) , p h ( t )) ∈ V h × Q h , u h (0) = u h 0 , u h ( t ) = g h ( t ) on Γ D such that � � ∂ u h � � � ̺ f ∂ t + u h · ∇ u h · v + 2 µ D ( u h ) : ∇ v − p h div v Ω Ω Ω � � = f · v + d · v Ω Γ N � q div u h = 0 Ω ∀ ( v , q ) ∈ V h 0 × Q h

Algebraic formulation (b) φ j ◮ { φ i } N u i =1 : basis of V h 0 φ i ψ l } N b ◮ { φ ( b ) j =1 : basis of V h \ V h 0 (shape functions u j corresponding to boundary nodes) ◮ { ψ l } N p l =1 : basis of Q h Expand the solution ( u h ( t ) , p h ( t )) on the finite element basis N b � N u � u g i ( t ) φ ( b ) u h ( x , t ) = u i ( t ) φ i ( x ) + ( x ) j i =1 j =1 � �� � = g h known term N p � p h ( x , t ) = p i ( t ) ψ l ( x ) l =1 Vectors of unknown dofs (nodal values if Lagrange basis functions are used) U ( t ) = [ u 1 ( t ) , . . . , u N u ( t )] T , P ( t ) = [ p 1 ( t ) , . . . , p N p ( t )] T

Algebraic system Replacing in the Galerkin formulation the expansions of u h and p h on the finite element basis and testing with shape functions φ i and ψ l we are led to the following system of ODEs  M d U  dt + A U + N ( U ) U + B T P = F u ( U ) , t > 0  B U = F p , with � � M ij = ̺ f φ j φ i mass matrix A ij = 2 µ D ( φ j ) : ∇ φ i , stiffness matrix Ω Ω � N ( U ) ij = !! depends on U ; non linear term u h · ∇ φ j · φ i , Ω � B li = − ψ l div φ i , divergence matrix Ω � � � � ∂ g h � � ( F u ( U )) i = f · φ i + d · φ i − ̺ f ∂ t + u h · ∇ g h · φ i − 2 µ D ( g h ) : ∇ φ i Ω Γ N Ω Ω � ( F p ) l = div g h ψ l Ω

A simpler problem – Stokes equation Let us consider for the moment the steady state Stokes problem  − 2 µ div D ( u ) + ∇ p = f ,  in Ω ,    div u = 0 , in Ω ,  u = g , on Γ D ,    (2 µ D ( u ) − p I ) · n = d , on Γ N , Weak formulation: find ( u , p ) ∈ [ H 1 ] d × L 2 , u = g on Γ D , such that � ∀ v ∈ [ H 1 Γ D ] d a ( u , v ) + b ( v , p ) = F ( v ) , ∀ q ∈ L 2 b ( u , q ) = 0 , with � � a ( u , v ) = 2 µ D ( u ) : ∇ v , b ( v , p ) = − p div v , Ω Ω � � F ( v ) = f · v + d · v Ω Γ N

Stokes problem The Stokes problem is well posed (admits a unique solution ( u , p ) which depends continuously on the data). In particular, the bilinear form b ( · , · ) satisfies the important property (LBB condition) b ( v , q ) q ∈ L 2 sup inf � v � H 1 � q � L 2 ≥ β > 0 v ∈ H 1 Γ D which implies that ∀ q ∈ L 2 , ∃ v q ∈ H 1 Γ D : b ( v q , q ) > 0. Finite elements approximation: find ( u h , p h ) ∈ V h × Q h , u h = g h on Γ D , such that � a ( u h , v h ) + b ( v h , p h ) = F ( v h ) , ∀ v h ∈ V h , 0 (*) b ( u h , q h ) = 0 , ∀ q ∈ Q h Corresponding algebraic formulation � � � � � � � A U + B T P = F u B T A U F u = ⇒ = B 0 P F p B U = F p The algebraic system has the typical structure of a saddle point problem

Spurious pressure modes If these exists a pressure p ∗ h such that b ( v h , p ∗ h ) = 0 , ∀ v h ∈ V h , this pressure will not be seen by the first equation of (*). It follows that, if ( u h , p h ) is a solution of (*), then also ( u h , p h + p ∗ h ) is a solution to the system and we loose uniqueness of the pressure. Such a pressure p ∗ h is called a spurious pressure mode A necessary and sufficient condition to avoid the presence of spurious pressure modes is that the finite elements spaces V h × Q h satisfy the b ( v h , q h ) inf-sup condition: q h ∈ Q h sup inf � v h � H 1 � q h � L 2 ≥ β h > 0 v h ∈ V h 0 Observe that the inf-sup condition is satified by the continuous spaces V 0 × Q but not necessarily by the finite element spacese V h , 0 ⊂ V 0 and Q h ⊂ Q . At the algebraic level, the inf-sup condition is equivalent to the request that Ker( B T ) = ∅

Spurious pressure modes Spaces that satisfy the inf-sup condition are called compatible. Assume the spaces ( V h , 0 , Q h ) are compatible with a constant β h independent of h . Then, the optimality of the Galerkin projection holds: � � � u h − u ex � H 1 + � p h − p ex � L 2 ≤ C v h ∈ V h � u ex − v h � H 1 + inf inf q h ∈ Q h � p ex − q h � L 2 Examples of spaces that satisfy the inf-sup condition P bubble P 2 / P 1 / P 1 Q 2 / Q 1 1 Examples of spaces that do not satisfy the inf-sup condition P 1 / P 1 P 1 / P 0 Q 1 / Q 1

Examples of spurious modes u=1 Couette flow: .n=0 σ .n=0 exact solution: u 1 = y ; u 2 = p = 0 σ u=0 Pressure Pressure IsoValue IsoValue -9.94024 -3.7491 -8.52021 -3.21351 -7.57352 -2.85645 -6.62683 -2.4994 -5.68015 -2.14234 -4.73346 -1.78528 -3.78677 -1.42822 -2.84008 -1.07116 -1.8934 -0.714103 -0.94671 -0.357045 -2.32106e-05 1.34456e-05 0.946664 0.357072 1.89335 0.71413 2.84004 1.07119 3.78673 1.42825 4.73341 1.78531 5.6801 2.14236 6.62679 2.49942 7.57347 2.85648 9.94019 3.74913 Pressure computed with P 1 / P 0 Pressure computed with P 1 / P 1 The issue of spurious pressure modes appears in the full Navier-Stokes equations as well. It is related to the incompressibility contraint.

Back to the Navier-Stokes equations: Temporal discretization A simple time integration scheme: Implicit Euler with semi-implicit treatment of the convective term:  u n − u n − 1  + u n − 1 · ∇ u n − div(2 µ D ( u n )) + ∇ p n = f n ̺ f ∆ t in Ω , n = 1 , 2 , . . .  div u n = 0 u 0 = u 0 , u n = g n on Γ D − p n n + 2 µ D ( u n ) · n = d n on Γ N In algebraic form this leads to a linear system to solve at every time step, of the form � 1 � � � � � 1 ∆ t M + A + N ( U n − 1 ) B T U n F n ∆ t M U n − 1 u + = P n F n B 0 p Observe the structure of the matrix, equal to the one for the Stokes problem.

Projection methods The idea of projection methods is to avoid the costly solution of the saddle point problem and to split the computation of velocity and pressure at each time step. Fractional step approach (assume homogeneous Dirichlet boundary cond.) ∂ u NS eqs. ̺ f ∂ t + L 1 u + L 2 u = f ���� ���� viscous+transport terms incompress. constraint Split the computation into  u n − u n − 1 ˜  u n = f n ̺ f + L 1 ˜ Step I: (transport-diffusion eq.) ∆ t  u n | ∂ Ω = 0  u n − ˜ u n  + ∇ p n = 0  ̺ f  ( L 2 -proj. on divergence ∆ t Step II: div u n = 0 free functions)    u n | ∂ Ω · n = 0

Some numerical issues ◮ High Reynolds flow may need stabilization of the convective term to avoid spurious oscillations in the velocity ◮ Non-conforming discretization may be computationally attractive ⇒ stabilized schemes to circumvent the LBB condition ◮ The numerical problem is usually quite large ( > 10 6 unknown). An iterative solution of the linear system is mandatory → good preconditioners are required ◮ To afford large computation in acceptable times we need parallel computing (and parallel, scalable preconditioners as well) ◮ Turbulence modeling is still an active research problem

Some references on the Navier-Stokes problem A.J. Chorin, J.E. Marsden. A mathematical introduction to fluid mechanics. Third edition. Springer-Verlag, 1993 H. Elman, D. Silvester, A. Wathen. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Oxford University Press, 2005 J.H. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Springer, 2002 J.L. Geurmond, P. Minev, J. Shen An overview of projection methods for incompressible flows , Comput. Methods Appl. Mech. Engrg. 195 (2006) 60116045 P.M. Gresho, R.L. Sani. Incompressible Flow and the Finite Element Method. (Two volumes), Wiley, 2000. W. Layton. Introduction to the numerical analysis of incompressible viscous flows. Society for Industrial and Applied Mathematics (SIAM), 2008. A. Prohl Projection and quasi-compressibility methods for solving the incompressible Navier-Stokes equations. B.G. Teubner, 1997 A. Quarteroni. Numerical models for differential problems . Springer-Verlag, 2009.

Free surface flow

Free Surface Flow Γ s �������� �������� �������� �������� �������� �������� �������� �������� �������� �������� �������� �������� �������� �������� Γ �������� �������� c �������� �������� �������� �������� �������� �������� �������� �������� Γ �������� �������� �������� �������� open Γ b γ fo We consider a fluid domain composed by ◮ the bottom surface Γ b ( t ) = { ( x , y , − h ( t )) | ( x , y ) ∈ Ω } ◮ the free surface Γ s ( t ) = { ( x , y , η ( t )) | ( x , y ) ∈ Ω } ◮ the lateral boundary Γ l ( t ) = { ( x , y , z ) | ( x , y ) ∈ γ f , z ∈ ( − h , η ( t )) } .

The Navier-Stokes equations (revisited) � � Du Dt = − 1 ∂ p ∂ x + ∇ · ( ν ∇ u ) + ∂ ν ∂ u + fv , ρ ∂ z ∂ z � � Dv Dt = − 1 ∂ p ∂ y + ∇ · ( ν ∇ v ) + ∂ ν ∂ v − fu , ρ ∂ z ∂ z � � Dw Dt = − 1 ∂ p ∂ z + ∇ · ( ν ∇ w ) + ∂ ν ∂ w − g , ρ ∂ z ∂ z ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z = 0 , ◮ ∇ is the 2D Laplace operator (while � ∇ is the 3D operator); ◮ ν is the kinematic viscosity coefficient, ν = µ/ρ ◮ f the Coriolis coefficient ◮ D / Dt = ∂/∂ t + u · � ∇

Mathematical models Traditionally methods used for free surface dynamics fall into one of these categories Integral equations boundary elements or panel methods. Interface capturing methods Volume of fluid. MAC methods. We solve the equations for both water and air. We can account for overturning and breaking waves. Interface tracking methods. The surface is described by a function η = η ( x , y , t ). We need to solve the equation only for the water, but we cannot describe overturning waves

The VOF method The domain Ω comprised both water and air. The interface is tracked by means of an additional scalar variable c which is transported with the fluid ∂ c ∂ t + ( u · ∇ ) c = 0 in Ω with initial datum � 1 in water c = 0 in air Then ρ = c ρ w + (1 − c ) ρ air µ = c µ w + (1 − c ) µ air Remark: Care must be taken in the discretization of the transport equation to avoid excessive “smearing” of the interface.

Interface tracking: conditions on Γ s ( t ) ∂η ∂ t + u ∂η ∂ x + v ∂η ∂ y = w , (kinematic condition) ν ∂ u ∂ z = ( ρ a C w | W | W ) (dynamic condition) where ◮ W is the wind velocity ◮ ρ a is the air density ◮ (Dean, Dalrymple, 1991)  1 . 2 × 10 − 6 , if | W | ≤ W c (= 5 . 6 m / s ) ,  � � 2 C w = 1 − W c 1 . 2 × 10 − 6 + 2 . 25 × 10 − 6  , if | W | > W c , | W |

Conditions on Γ b (fixed bottom) and Γ l u ∂ h ∂ x + v ∂ h ∂ y + w = 0 (kinematic condition) ν ∂ u ∂ z = f 2 ( u ) , (dynamic condition) The function f 2 accounts for the friction on the bottom and may take a similar form to that used on the free surface. It approximates the action of the boundary layer at the bottom that in this type of formulation is usually not resolved. On Γ l we may have any admissible condition for a Navier-Stokes problem. However, computationally it may be convenient to devise non reflective conditions for the elevation.

Integral form of the free surface equation We integrate the continuity equation along the vertical coordinate � η � η ∂ u ∂ v ∂ y d z + w | Γ s = 0 ∂ x d z + − h − h Using the kinematic condition at free surface ⇒ � η � η ∂η ∂ t + ∂ u d z + ∂ v d z = 0 ∂ x ∂ y − h − h Remark: the free surface is represented by a single valued function of x and y and t : no overturning or breaking waves are allowed in this formulation.

The non-hydrostatic model We split the pressure into a constant term (atmospheric pressure, may be set to zero) a hydrostatic component and a non-hydrostatic correction p = p a + ρ g ( η − z ) + gq . � � D u Dt −∇ · ( ν ∇ u ) − ∂ ν ∂ u + ∇ q + g ∇ η = f ∂ z ∂ z � � Dt −∇ · ( ν ∇ w ) − ∂ ν ∂ w + ∂ q Dw ∂ z = 0 , ∂ z ∂ z ∇ · u + ∂ w ∂ z = 0 � η ∂η ∂ t + ∇ · u d z = 0 − h When the ratio between basin height and wave length H / L << 1 horizontal viscous term are small compared to the vertical ones (Grener et al 2000), so may be neglected.

The 2D equivalent: Boussinesq equations With certain manipulations a 2D approximation may be devised: ∂η ∂ t + ∇ · [( η + h ) u ] = 0 , ∂ t ) − h 2 ∂ u ∂ t + u · ∇ u + g ∇ η = h 2 ∇∇ · ( h ∂ u 6 ∇∇ · ∂ u ∂ t , Eliminating the terms in blue we obtain the shallow water equations. Pressure is recovered via Bernoulli’s relation p = g ρ ( η − z ) + ρ 2(2 hz + z 2 ) ∇ · ∂ u ∂ t . Remark: However, for boat dynamics the 2D equations are inadequate.

Waves produced by a ship in a channel Hydrostatic case Non hydrostatic case

Lagrangian treatment of time derivative Dt ( t n +1 ) ≃ u n +1 − u n ( X ) We use the approximation D u where the foot of ∆ t the characteristic line X ( s ; t , x ) is obtained by solving an ODE d X ( s ; t , x ) = V ( s , X ( s ; t , x )) for s ∈ (0 , t ) , ds X ( t ; t , x ) = x . Discretization: backward (composite) Euler (other possibility: Runge-Kutta scheme) . . . x . . . . X Care must be taken if X falls outside the domain.

How to account for the presence of the boat? A possibility is to describe the boat external shape by means of a suitable function Ψ = Ψ( x , y , t ) and impose that η ≤ Ψ at all times. Ψ η γ Ψ ω

Some References on free-surface and ship hydrodynamics R. Azcueta. Computation of turbulent free-surface flows around ships and floating bodies. Ship Technology Research , 49(2):4669, 2002 U. Bulgarelli, C. Lugni, M. Landrini. Numerical modelling of free-surface flows in ship hydrodynamics Int. J. Numer. Meth. Fluids , 43(5):465-481, 2003. C.C. Mei. The applied dynamics of ocean surface waves. World Scientific Publishing, Singapore, 1989 E. Miglio and A. Quarteroni. Finite element approximation of Quasi-3D shallow water equations. Comput. Methods Appl. Mech. Engrg. 174(3-4):355-369 , 1999 N. Parolini and A. Quarteroni. Mathematical models and numerical simulations for the Americas Cup Comput. Methods Appl. Mech. Engrg. 194(9-11):1001-1026, 2005 G. B. Witham. Linear and non linear waves, Wiley, 1974

Structural models and their finite element approximation

Structural models – kinematics Dynamics of structures is more conveniently described in Lagrangian coordinates, i.e. equations are written in the reference configuration Ω s 0 (e.g. the initial configuration). Lagrangian map : x = x ( ξ , t ) = L t ( ξ ) ◮ ξ : coordinate of material point ξ − in reference configuration Ω s = x Ω s η t 0 ◮ x : coordinate of material point x in deformed configuration ξ ◮ η ( ξ , t ) = x − ξ displacement of material point Kinematics u ( ξ , t ) = ∂ x velocity : ∂ t ( ξ , t ) = ˙ x ( ξ , t ) = ˙ η ( ξ , t ) a ( ξ , t ) = ∂ 2 x ∂ t 2 ( ξ , t ) = ¨ acceleration : x ( ξ , t ) = ¨ η ( ξ , t )

Eulerian versus Lagrangian Lagrangian frame (mostly used for solids) Kinematics is described in the reference configuration position x = x ( ξ , t ) velocity u = u ( ξ , t ) = ˙ x ( ξ , t ) acceleration a = a ( ξ , t ) = ¨ x ( ξ , t ) Eulerian frame (mostly used for fluids) Kinematics is described in the current (deformed) configuration position (trajectories) x ( ξ , t ) : x = u ( x , t ) , x ( ξ , 0) = ξ ˙ velocity u = u ( x , t ) a = a ( x , t ) = D u Dt ( x , t ) = ∂ u acceleration ∂ t + u · ∇ x u

Measures of strain Let us introduce the deformation gradient tensor F ij = ∂ x i F = ∇ ξ x = I + ∇ ξ η ∂ξ j Given an infinitesimal material line segment d ξ in Ω s 0 , it will be transformed through the Lagrangian map into d x = F d ξ . ◮ Right Cauchy-Green tensor C = F T F Measures the change in � d ξ � due to the motion. � d x � 2 � d ξ � 2 = d ξ T F T F d ξ = d ξ T � d ξ � 2 = d x T · d x � d ξ � C d ξ � d ξ � = v T Cv = λ 2 ( v ) � d ξ � 2 with v = d ξ / � d ξ � unitary vector. x=x( ,t) ξ The change in length of d ξ depends on its orientation. v||d || ξ λ ξ (v)||d || s s Ω 0 Ω t

Equations of motion Equations in the reference configuration ( V t is a materal volume) � � � � ◮ d ̺ s u d x = d ̺ 0 ̺ 0 J ̺ s u d ξ = s a d ξ = s ¨ η d ξ ���� dt dt V t V 0 V 0 V 0 ̺ 0 s � � f d x = J ( ξ ) f ( x ( ξ , t )) d ξ ◮ � �� � V t V 0 f 0 ( ξ ) ◮ � � � � Nanson’s f. J σ s F − T div ξ σ 0 t dA = σ s n dA = n 0 dA 0 = s d ξ � �� � ∂ V t ∂ V t ∂ V 0 V 0 σ 0 s ̺ 0 η − div ξ σ 0 s = f 0 momentum equation s ¨ where we have introduced the first Piola-Kirchhoff stress tensor (also called nominal stress tensor) σ 0 s = J σ s F − T . Observe that the divergence div ξ is taken with respect to the reference coordinate ξ .

Constitutive relations Differently than for fluids, in an elastic body the internal stress depends on of deformation gradients (instead of gradients of deformation rates). Cauchy elastic material ◮ σ 0 σ 0 s = σ 0 s is a continuous function of F = ⇒ s ( F ) s ( F ) Q T for all Q orthogonal. ◮ Frame indifference σ 0 s ( QF ) = σ 0 = ⇒ Implies that σ 0 s can be expressed as a funtion of C = F T F only: σ 0 s = σ 0 s ( C ) ◮ Incompressible materials: one adds the constraint det F = 1 and a corresponding pressure p (Lagrange multiplier) s = � σ 0 s − pJ F − T σ 0 σ s = ˜ σ s ( η ) − p I , = ⇒

Examples of constitutive relations Incompressible materials ( J = 1) ◮ Neo-Hookean: W ( C ) = µ 2 (tr C − 3) s = 2 F ∂ W ∂ C − pJ F − T = µ F − p Cof F σ 0 ◮ Mooney-Rivlin: W ( C ) = µ 1 2 (tr C − 3) − µ 2 2 ( I 2 ( C ) − 3) � � s = 2 F ∂ W ∂ C − pJ F − T = σ 0 µ 1 + µ 2 tr( F T F ) F − µ 2 FF T F − p Cof F Compressible materials 2 tr( E ) 2 + µ tr( E 2 ) ◮ Saint-Venant Kirchhoff: W ( C ) = λ with E = 1 2 ( C − I ) and ( λ, µ ) Lam´ e constants s = 2 F ∂ W ∂ C = F ∂ W σ 0 ∂ E = F ( λ (tr E ) I + 2 µ E )

Equations of elasto-dynamics Compressible materials s = 2 F ∂ W ̺ 0 η − div ξ σ 0 s ( η ) = f 0 , σ 0 s ¨ ∂ C Incompressible materials � s = 2 F ∂ W ̺ 0 η − div ξ σ 0 s ( η , p ) = f 0 , σ 0 s ¨ ∂ C − p Cof F J = 1 Possible boundary conditions ◮ imposed displacement: η = g on Γ D ◮ imposed nominal stress: σ 0 s n 0 = d on Γ N Example: incompressible Neo-Hookean material  ̺ 0 η − div ξ [ µ F ( η ) − p Cof F ( η )] = f 0 ( η ) in Ω 0 s ¨   s   in Ω 0 J ( η ) = 1 s  η = g on Γ D    σ 0 s n 0 = d on Γ N

Finite element approximation – compressible case Consider a finite element space V h ⊂ V (e.g. piecewise continuous polynomials) defined on a suitable triangulation T h of Ω s 0 . Finite element formulation: find η h ( t ) ∈ V h , η h = g h on Γ D , such that � � � � ∂ 2 η h ̺ 0 σ 0 f s ∂ t 2 · φ h + s ( η h ) : ∇ ξ φ h = 0 · φ h + d · φ h dA 0 s Ω s Ω s Ω s Γ N 0 0 0 for all φ h ∈ V h , vanishing on Γ D Introduce a basis { φ i } N s i =1 of V h and expand the displacement η h on the η h ( t , ξ ) = � N s basis: i =1 η i ( t ) φ i ( ξ ). Then, we obtain a system of non-linear ODEs M s ¨ d s + K s ( d s ) = F s , t > 0 (1) with ◮ Solution vector: d s ( t ) = [ η 1 ( t ) , . . . , η N s ( t )] T � ◮ Mass matrix: ( M s ) ij = 0 ̺ 0 s φ i φ j Ω s � 0 σ 0 ◮ Stiffness (nonlinear) term: ( K s ( d s )) i = s ( η h ) : ∇ ξ φ i Ω s

Time discretization - Newmark scheme A popular numerical scheme for system (1) the second-order Newmark s and ˙ scheme, were truncated Taylor expansions for d n d n s are used to express the acceleration ¨ d n s in terms of the displacement d n s . 1 1 ) + 1 − 2 β with ζ n = ¨ s − ζ n , β ∆ t 2 ( d n − 1 + ∆ t ˙ d n − 1 d n − 1 ¨ d n β ∆ t 2 d n s = s s s 2 β We obtain 1 s + M s ζ n β ∆ t 2 M s d n s + K s ( d n s ) = F n At each time step we have to solve a nonlinear system in the displacement vector d n s . This can be done, for instance, with a Newton-Rapson iteration. The scheme is unconditionally stable for β > 1 / 4.

Some references on elasticity and its numerical approximation J. Bonet, R.D. Wood. Nonlinear continuum mechanics for finite element analysis. Second edition. Cambridge University Press, 2008. T. Belytschko, W.K. Liu, B. Moran. Nonlinear finite elements for continua and structures. John Wiley & Sons, Ltd., Chichester, 2000. P.G. Ciarlet Mathematical Elasticity, Vol. I , Elsevier, 2004 O.C. Zienkiewicz, R.L. Taylor. The finite element method. Vol. 2. Solid mechanics. Fifth edition. Butterworth-Heinemann, Oxford, 2000 R.W. Ogden. Nonlinear elastic deformations. Halsted Press [John Wiley & Sons, Inc.], New York, 1984

Coupled fluid-structure problem

Eulerian versus Lagrangian description We consider now an incompressible fluid interacting with an elastic structure featuring possibly large deformations. As we have seen, ◮ Fluid equations are typically written in Eulerian form ◮ Structure equations are typically written in Lagrangian form on the reference domain Ω s 0 . WARNING : The fluid equations are defined in a moving domain Ω f t Current Configuration Reference configuration (at time t) Ω s Ω s 0 η t Γ t Γ 0 A) Ω f Ω f 0 t ( u ,p) x 2 ζ 2 x 1 ζ 1

Fluid structure coupling conditions At the common interface Γ t (evolving in time), we impose ◮ Continuity of velocity (kinematic condition) ◮ Continuity of the normal stress (dynamic condition) Warning : the continuity of stresses on Γ t involves the physical stresses, i.e. the Cauchy stress tensors. u ( t , x ) = ∂ η (cont. velocity) ∂ t ( t , ξ ) , with x = L t ( ξ ) σ f ( u , p ) n f = − σ s ( η ) n s (cont. normal stress) J ( η ) σ 0 1 s ( η ) F T ( η ) between the Cauchy Remember the relation σ s ( η ) = stress tensor and the nominal stress tensor. n n f Ω s t η Normal convention: n f = − n s Γ t f n s Ω t

The coupled fluid-structure problem  ∂ u  ∂ t + ̺ f div( u ⊗ u ) − div σ f ( u , p ) = f f  ̺ f in Ω f t ( η ) ,   div u = 0 ∂ 2 η ̺ 0 ∂ t 2 − div ξ σ 0 s ( η ) = f s in Ω s 0 , 0 , s u ( t , x ) = ∂ η with x = L t ( ξ ) ∂ t ( t , ξ ) , on Γ t σ f ( u , p ) n f = − σ s ( η ) n s . on Γ t

Sources of non-linearity The coupled fluid-structure problem is highly non-linear. Sources on non-linearity are: ◮ Convective term div( u ⊗ u ) of the Navier-Stokes equations ◮ Non-linearity in the structure model when using non-linear elasticity ◮ The fluid domain is itself an unknown

Energy inequality • Taking as test functions ( v , q , φ ) = ( u , p , ˙ η ) in the weak formulation of the coupled problem we can derive an Energy inequality Fluid Structure Energy (kinetic + elastic) t ) + ̺ 0 E FS ( t ) ≡ ̺ f 2 � ∂ η 2 � u ( t ) � 2 ∂ t ( t ) � 2 s 0 ) + W ( η ( t )) L 2 (Ω f L 2 (Ω s Then, for an isolated system (no external forces) � T � E FS ( T ) + 2 µ D ( u ) : D ( u ) d Ω dt ≤ E FS (0) Ω f 0 t � T � The term 2 µ t D ( u ) : D ( u ) d Ω dt represents the energy dissipated Ω f 0 by the viscosity of the fluid

Energy Inequality Key points in deriving an energy inequality: ◮ Perfect balance of work at the interface � � 0 ) · ∂ η ( σ f · n f ) · u = − ( σ 0 s ( η ) · n s ∂ t Γ t Γ 0 ◮ No kinetic flux through the interface � � � ∂ u ∂ t · u = ̺ f d | u | 2 − ̺ f | u 2 | w · n (time der.) ̺ f 2 dt 2 Ω f Ω f Γ t � t t � ̺ f div( u ⊗ u ) · u = ̺ f | u 2 | u · n (convective term) 2 Ω f Γ t t where w is the velocity at which the interface moves. � η , the kinetic flux ̺ f Γ t | u 2 | ( u − w ) · n vanishes. Since w = u = ˙ 2

Some References [DGeal04] Q. Du, M.D. Gunzburger, L.S. Hou, J. Lee. Semidiscrete finite element approximations of a linear fluid-structure interaction problem. SIAM J. Numer. Anal. 42(1):1–29, 2004. [DGeal03] Q. Du, M.D. Gunzburger, L.S. Hou, J. Lee. Analysis of a linear fluid-structure interaction problem. Discrete Contin. Dyn. Syst. 9(3)633–650, 2003. [LeTM01] P. Le Tallec and J. Mouro. Fluid structure interaction with large structural displacements. Comput. Methods Appl. Mech. Engrg. , 190:3039–3067, 2001. [N01] F. Nobile. Numerical approximation of fluid-structure interaction problems with application to haemodynamics , PhD Thesis n.2458, ´ Ecole Polytechnique F´ ed´ erale de Lausanne, 2001.

ALE formulation for fluids in moving domains

How to deal with moving domains The fluid equations are defined on a moving domain. How can we treat them numerically? Idea Use a moving mesh, that follows the deformation of the domain. T s w Γ t,h h f f f T T h,0 t,h A t The mesh follows the moving boundary and is fixed on the artificial (lateral) sections.

Eulerian / Lagrangian / ALE frame of reference Γ Γ w w t t L f f Ω Ω L Ω Ω A out out in Γ Γ in Γ t Γ t t t t L t t L t Γ w Γ w t L t (ALE) (Lagrangian) t t A L Γ w 0 Γ in f Γ out Ω u u 0 0 0 Γ w 0 Neither an Eulerian nor a Lagrangian frame are suited to our case. We want to follow the interface in a Lagrangian fashion, while keeping an “Eulerian-type” description of the remaining part of the boundary. A possible answer: Arbitrary Lagrangian Eulerian (ALE) frame: Introduce a (fixed) frame of reference which is mapped at every time to the desired physical domain. The equations of motion are then recast in such reference frame.

✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ☎ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ✆ ☎ ☎ ✆ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ✆ ✆ ✆ ✆ ✆ ✆ � �✁ ✆ ☎ ✆ ✂✄ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ✆ ALE formulation (assuming known the domain deformation) The moving domain is recast at each time t to a fixed configuration Ω f 0 through the ALE mapping A t : ✁ x=x (t, ξ ) A t Ω 0 A t : Ω f → Ω f 0 − t , ξ Ω t x ( ξ , t ) = A t ( ξ ) w ( x , t ) = ∂ A t ◮ domain velocity ∂ t ◦ A − 1 t ( x ) � � � � ∂ u = ∂ u � � ◮ ALE derivative + w · ∇ x u � � ∂ t ∂ t x ξ � � ∂ J A t � ◮ Euler expansion = J A t div w , J A t = det( ∇A t ) � ∂ t ξ ◮ ALE Transport formula � � � � � � d ∂ u � , ∀ V t ⊂ Ω f u ( x , t ) = + u div w � t dt ∂ t V t V t ξ

Navier Stokes equations in ALE form  � � ∂ u   � ̺ f + ̺ f (( u − w ) · ∇ ) u − div σ f ( u , p ) = f � ∂ t in Ω t ξ   div u = 0 Hybrid formulation: the spatial derivative terms are computed on the deformed configuration Ω f t , whereas the time derivative is computed on the reference configuration Ω f 0 . Using the Euler formula, the momentum equation can be equivalently written in conservative form � � ̺ f ∂ J A t u � + ̺ f div(( u − w ) ⊗ u ) − div σ f ( u , p ) = f � J A t ∂ t ξ

NS-ALE weak formulation 0 )] d defined on the reference v ∈ [ H 1 0 (Ω f Consider a test function ˆ configuration, and let us “map” it onto the deformed configuration: v ◦ ( A t ) − 1 . v = ˆ We multiply the momentum eq. by v and integrate by parts Non conservative formulation � � � � � � ∂ u � ̺ f + ( u − w ) · ∇ u · v + σ f : ∇ v + div u q = f · v � ∂ t Ω f Ω f ξ t t Again, we can write an alternative conservative formulation by using the ALE transport formula Conservative formulation � � � d ̺ f u · v + ̺ f div(( u − w ) ⊗ u ) · v + σ f : ∇ v + div u q = f · v dt Ω f Ω f Ω f t t t

Finite Element ALE formulation We introduce a triangulation T h , 0 of the reference domain Ω f 0 and choose finite element spaces ( V h , 0 , Q h , 0 ) on T h , 0 for velocity and pressure resp. T s w Γ t,h h f f f T h,0 T t,h A t At time t , let T h , t be the image of T h , 0 through the mapping A t . Similarly, we define finite element spaces ( V h , t , Q h , t ) on T h , t as the map of ( V h , 0 , Q h , 0 ) through A t . Remark. Assume A t is a piecewise linear function on T h , 0 . Then, T h , t is a triangulation of Ω f t . Moreover, if V h , 0 (resp Q h , 0 ) is the space of piecewise polynomials of degree p on T h , 0 , then V h , t (resp Q h , t ) is the space of piecewise polynomials of degree p on T h , t

Finite Element ALE formulation – conservative find ( u h , p h ) ∈ ( V h , t , Q h , t ), for each t > 0, such that � � � d ̺ f u h · v h + ̺ f div(( u h − w h ) ⊗ u h ) · v h + σ f : ∇ v h +div u h q h = f · v h dt Ω f Ω f Ω f t t t for all ( v h , q h ) ∈ ( V h , t , Q h , t ) Algebraic formulation of the momentum equation d + B T ( t ) P ( t ) = F ( t ) dt ( M ( t ) U ( t )) + A ( t ) U ( t ) + N c ( t , U − W ) U ( t ) � �� � � �� � � �� � stiffness transport ALE time der. � � � with N c ( U ) ij = Ω div ( u h − w h ) ⊗ φ j · φ i ◮ The two formulations (conservative and non conservative) are equivalent at the time-continuous level. But they lead to different time discretization schemes. ◮ Observe that all matrices depend on time!

How to construct the ALE mapping in practice T s w Γ t,h h f f f T T h,0 t,h A t Assume the displacement η n of the boundary is known at time t n . We can solve for instance a Laplace equation for the displacement of the internal nodes � ∆ x n = 0 , in Ω f 0 x n = ξ + η n , on ∂ Ω f 0 In practice, we solve it by the finite element method: x n = � N f i =1 x n i φ i , and X n = [ x n N f ] T satisfying 1 , . . . , x n K m X n = G n Remark. the computed mapping x n ( ξ ) = A t n ( ξ ) will be a piecewise polynomial (for instance piecewise linear)

Some References [D83] J. Donea, An arbitrary Lagrangian-Eulerian finite element method for transient fluid-structure interactions Comput. Methods Appl. Mech. Engrg. , 33:689–723:1982. [FGG01] C. Farhat, P. Geuzainne, C. Grandmont, The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids, J. Comput. Physics 174:669–694, 2001. [FN99] L. Formaggia, F. Nobile, A stability analysis for the Arbitrary Lagrangian Eulerian formulation with finite elements, East-West J. Num. Math. , 7:105–132, 1999. [FN04] L. Formaggia, F. Nobile, Stability analysis of second-order time accurate schemes for ALE-FEM, Comput. Methods Appl. Mech. Engrg. , 193:4097–4116, 2004. [EGP09] S. ´ Etienne, A. Garon, D. Pelletier, Perspective on the geometric conservation law and finite element methods for ALE simulations of incompressible flow, J. Comput. Physics , 228:2313–2333, 2009. [JLBL06] J-F. Gerbeau, C. Le Bris, T. Leli` evre, Mathematical Methods for the Magnetohydrodynamics of Liquid Metals (Ch. 5), Oxford Univ. Press, 2006

Some References [HLZ81] T.J.R. Hughes, W.K. Liu, T.K. Zimmermann, Lagrangian-Eulerian finite element formulation for incompressible viscous flows, Comput. Methods Appl. Mech. Eng. 29(3):329–349, 1981. [LF96] M. Lesoinne, C. Farhat, Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact in aeroelastic computations, Comput. Methods Appl. Mech. Engrg. 134:71–90, 1996. [N00] B. Nkonga, On the conservative and accurate CFD approximations for moving meshes and moving boundaries, Comput. Methods Appl. Mech. Eng. , 190:1801–1825, 2000. [TL79] P.D. Thomas, C.K. Lombard, Geometric conservation law and its application to flow computations on moving grids, AIAA 17:1030–1037, 1979.

Algorithms for fluid-structure interaction with large displacements

FSI problem - a three field formulation F l ALE ( w n ; u n , p n ) = 0 in Ω f Fluid problem n S t ( η n , ˙ η n ) = 0; in Ω s Structure problem 0 M esh ( A n , ˙ A n ) = 0 in Ω f ALE mapping 0 Coupling conditions u n ◦ A n = ˙ η n kinematic cond. on Γ 0 σ f ( u n , p n ) n f = − σ s ( η n ) n s dynamic cond. on Γ t A n = ξ + η n , A n = ˙ ˙ η n geometric cond. on Γ 0 w n = ˙ A n ◦ ( A n ) − 1 , Ω t n = A n (Ω f 0 ) Fully coupled non-linear problem in the unknowns ( u n , p n , A n , ˙ η n ) A n , η n , ˙

Partitioned versus Monolithic solvers ◮ We refer to Monolithic schemes when we try to solve the non-linear problem all at once. We assemble and solve a “huge” linear system in the unknowns u n , p n , A n , ˙ η n . A n , η n , ˙ ◮ We refer to Partitioned strategies when we solve iteratively the fluid and structure subproblems (never form the full matrix). This allows one to couple two existing codes, one solving the fluid equations and the other the structure problem. ◮ There are other options in between: for instance one can “formally” write the monolithic problem and then use a block preconditioner to solve it. In this way one recovers partitioned procedures (see e.g. [Heil, CMAME ’04])

Partitioned fluid-structure algorithms

A simple loosely coupled partitioned FSI algorithm Given the solution ( u n − 1 , p n − 1 , A n − 1 , ˙ η n − 1 ) at time A n − 1 , η n − 1 , ˙ step t n − 1 : 1. Extrapolate the displacement and velocity of the interface: e.g. η n = ˙ η n = η n − 1 + ∆ t ˜ ˜ η n − 1 , η n , ˙ ˜ ˙ on Γ 0 � M esh ( A n , ˙ A n ) = 0 in Ω f 0 2. Compute the ALE mapping: A n = ξ + ˜ A n = ˜ η n , ˙ η n ˙ on Γ 0 3. Solve the fluid problem with Dirichlet boundary conditions � F l ALE ( w n ; u n , p n ) = 0 , in Ω f n = A n (Ω f 0 ) u n = ˜ η n ◦ ( A n ) − 1 , on Γ n = A n (Γ 0 ) ˙ 4. Solve the structure equation with Neumann boundary conditions � η n ) = 0 , S t ( η n , ˙ in Ω s 0 σ s ( η n ) n s = − σ f ( u n , p n ) n f , on Γ n 5. Go to next time step

The staggered algorithm by Farhat and Lesoinne (’98) A popular algorithm proposed by C. Farhat and M. Lesoinne, CMAME 1998, uses staggered grids: ◮ The fluid equation is discretized with a second order backward difference scheme (BDF2) and colocated at t n +1 / 2 ◮ The ALE mapping is piecewise linear in [ t n − 1 / 2 , t n +1 / 2 ] ◮ The structure equation is discretized by the mid-point scheme and collocated at t n +1 . u n−1/2 p n−1/2 u n+1/2 p n+1/2 σ n+1/2 n−1/2 η f . . . n−1 n−1 n n n+1 n+1 η η η η η η

η n ): Given the solution ( u n − 1 / 2 , p n − 1 / 2 , A n − 1 / 2 , η n , ˙ 1. Extrapolate the interface displacement and velocity at t n +1 / 2 : η n +1 / 2 = η n + ∆ t η n +1 / 2 = ˙ η n , η n , 2 ˙ ˙ on Γ 0 � M esh ( A n +1 / 2 ) = 0 2. Compute the ALE mapping at t n +1 / 2 : A n +1 / 2 = ξ + η n +1 / 2 + A n − 1 / 2 t n +1 / 2 − t and A ( t ) = A n +1 / 2 t − t n − 1 / 2 t ∈ [ t n − 1 / 2 , t n +1 / 2 ] , ∆ t ∆ t 3. Solve the fluid problem with BDF2 at t n +1 / 2 � F l GCL ALE ( w ; u n +1 / 2 , p n +1 / 2 ) = 0 , in Ω f n +1 / 2 = A n +1 / 2 (Ω f 0 ) η n +1 / 2 ◦ ( A n +1 / 2 ) − 1 , u n +1 / 2 = ˙ on Γ n = A n (Γ 0 ) 4. Solve the structure equation with Mid-Point at t n +1 � η n +1 ) = 0 , S t ( η n +1 , ˙ in Ω s 0 σ s ( η n + η n +1 ) n s = − σ f ( u n +1 / 2 , p n +1 / 2 ) n f , on Γ n +1 / 2 2 5. Go to next time step

On partitioned procedures Partitioned procedures may have serious stability problems in certain applications. Dangerous situations appear for incompressible fluids when the mass of the structure is small compared to the mass of the fluid The instability is due to a non perfect balance of energy transfer at the interface. For instance, in the algorithm by Farhat & Lesoinne, we have: Work done in [ t n , t n +1 ] by � n + 1 2 ) n f ) · ∆ t u n + 1 ( σ f ( u n + 1 2 , p n + 1 fluid = W 2 2 f Γ n + 1 2 � η n +1 + ˙ ( σ s ( η n +1 + η n η n ) n s ) · ∆ t ˙ n + 1 structure = W 2 s 2 2 Γ n + 1 2 n + 1 n + 1 Ideally, one would like to have W 2 + W 2 = 0. s f n + 1 n + 1 Energy unbalancing: ∆ W = W + W 2 2 s f � η n +1 − ˙ η n 2 ) n f ) · ˙ ( σ f ( u n + 1 2 , p n + 1 ∆ W = − ∆ t 2 = O (∆ t 2 ) 2∆ t Γ n +1 / 2

A fully implicit partitioned FSI algorithm To have a perfect energy balance, we need to satisfy exactly at each time step the kinematic and dynamic coupling conditions. A possible remedy: fixed point approach: subiterate at each time step, eventually adding a relaxation step ◮ if � η n − ˜ η n � > tol η n ← ω η n + (1 − ω )˜ η n set ˜ repeat points 2-4 until convergence ◮ The fixed point algorithm is stable, in general. ◮ However, whenever the partitioned procedure is unstable, the fixed point may need a very small relaxation parameter ω and the convergence may be very slow.

A fully implicit partitioned FSI algorithm - Dir. Neumann Given the solution ( u n − 1 , p n − 1 , A n − 1 , η n − 1 , ˙ η n − 1 ) at time step t n − 1 , set η n − 1 and η n η n and for k = 1 , 2 , . . . 0 = η n − 1 + ∆ t ˙ η n ˙ 0 = ˙ 1. Compute the ALE mapping: � k , ˙ M esh ( A n A n k ) = 0 ˙ A n k = ξ + η n A n η n k − 1 , k = ˙ k − 1 , on Γ 0 2. Solve the fluid problem with Dirichlet boundary conditions � F l ALE ( w n k ; u n k , p n in Ω f n , k = A n k (Ω f k ) = 0 , 0 ) η n u n k − 1 ◦ ( A n k ) − 1 , on Γ n , k = A n k = ˙ k (Γ 0 ) 3. Solve the structure equation with Neumann boundary conditions � k , ˜ η n η n in Ω s S t (˜ ˙ k ) = 0 , 0 k ) n s = − σ f ( u n η n k , p n k ) n n , σ s (˜ on Γ n , k η n η n k − η n k − 1 � > tol , set η n k + (1 − ω ) η n 4. Relaxation step: If � ˜ k = ω ˜ k − 1 η n (and similarly for ˙ k ) and go to step 1.

Load transfer between fluid and structure At the discrete level, there is a proper way of exchanging information. i } N b Conforming finite element spaces at the interface: let { φ b i =1 be the basis functions (fluid and structure) corresponding to the interface nodes; N b N b � � u i φ b η i φ b u h | Γ = i | Γ , η h | Γ = i | Γ i =1 i =1 and U b = [ u 1 , . . . , u N b ], d b = [ η 1 , . . . , η N b ] are the vectors of dofs . ◮ Fluid Dirichlet b.c. u h = ˙ η h U b = ˙ d b = ⇒ nodal values u i = ˙ η i = ⇒ ◮ Structure Neumann b.c. σ s · n s = − σ f · n f � − σ f ( u h , p h ) n f · φ b = < Res( u h , p h ) , φ b F b s = − F b F s , i = i > = ⇒ i f � �� � Γ n − F f , i LOAD i ��� ��� ��� ��� ��� ��� Γ ��� ��� ��� ��� ��� ��� NS

Non-conforming finite element spaces at the interface N bf ◮ { φ b i } i =1 : basis functions (fluid) corresponding to the interface nodes. u h | Γ = � N bf i =1 u i φ b i | Γ j } N bs ◮ { ψ b j =1 basis functions (structure) corresponding to the interface nodes. η h | Γ = � N bs j =1 η j ψ b j | Γ η n Fluid Dirichlet b.cs: u n h = Π h ˙ h where Π h is a suitable operator (can be j = � N bf interpolation, L 2 projection, etc.). Let C ij such that Π h ψ b i =1 C ij φ b i | Γ N bs � U b = C ˙ d b = ⇒ u i = C ij ˙ η j = ⇒ j =1 Structure Neumann b.cs.: σ s n s = − Π ∗ h ( σ f n f ) (with Π ∗ h adjoint of Π h ) � N bs � ( σ f · n f ) · Π h ψ b C ij < Res( u h , p h ) , φ b F b s = − C T F b F s , j = − ⇒ j = i > , = f � �� � Γ i =1 − F f , i