A monolithic semi-implicit algorithm for fluid-structure interaction - PowerPoint PPT Presentation

A monolithic semi-implicit algorithm for fluid-structure interaction problem at small structural displacements Cornel Marius MUREA, joint work with Soyibou SY Laboratoire de Math´ ematiques, Informatique et Applications, Universit´ e de Haute-Alsace e-mail: cornel.murea@uha.fr http://www.edp.lmia.uha.fr/murea/ FreeFem++ Workshop, Paris, September 15, 2009

Initial (left) and intermediate (right) geometrical configuration Ω S A D A D Γ t B C B C Γ 0 Ω F Ω F Σ Σ Σ Σ 3 3 1 1 0 t Σ Σ 2 2 ∂ Ω F ∂ Ω F 0 = Σ 1 ∪ Σ 2 ∪ Σ 3 ∪ Γ 0 , t = Σ 1 ∪ Σ 2 ∪ Σ 3 ∪ Γ t .

Linear elasticity equations � � T : Ω S We denote by u S = 0 × [0 , T ] → R 2 the structure u S 1 , u S 2 displacement. ρ S ∂ 2 u S ∂ t 2 − ∇ · σ S f S , in Ω S = 0 × (0 , T ) λ S � ∇ · u S � � u S � σ S I + 2 µ S ǫ = � ∇ u S � T � � u S � � 1 ∇ u S + ǫ = 2 u S = 0 , on Γ D × (0 , T ) σ S n S = 0 , on Γ N × (0 , T ) Γ D = [ AB ] ∪ [ CD ] , Γ N = [ DA ] .

Navier-Stokes equations We denote by v F the fluid velocity and by p F the fluid pressure. � ∂ v F � ∂ t + ( v F · ∇ ) v F ρ F − ∇ · σ F f F , ∀ t ∈ (0 , T ) , ∀ x ∈ Ω F = t ∇ · v F ∀ t ∈ (0 , T ) , ∀ x ∈ Ω F = 0 , t � v F � σ F − p F I + 2 µ F ǫ = σ F n F = on Σ 1 × (0 , T ) h in , σ F n F = h out , on Σ 3 × (0 , T ) v F = 0 , on Σ 2 × (0 , T )

Interface and initial conditions The interface Γ t is the image of Γ 0 via the map X → X + u S ( X , t ) . Interface conditions v F � � ∂ u S X + u S ( X , t ) , t = ∂ t ( X , t ) , ∀ ( X , t ) ∈ Γ 0 × (0 , T ) � σ F n F � � σ S n S � = − ( X , t ) , ∀ ( X , t ) ∈ Γ 0 × (0 , T ) ( X + u S ( X , t ) , t ) Initial conditions u S ( X , t = 0) u 0 ( X ) , in Ω S = 0 ∂ u S u 0 ( X ) , in Ω S ∂ t ( X , t = 0) = ˙ 0 v F ( x , t = 0) v 0 ( x ) , in Ω F = 0

Arbitrary Lagrangian Eulerian (ALE) framework. Notations Ω F = Ω F The reference fluid domain � n , the interface Γ n The velocity of the fluid mesh ϑ n = ( ϑ n 2 ) T is the solution of 1 , ϑ n ∆ ϑ n = 0 in Ω F ϑ n = 0 on ∂ Ω F ϑ n = v F , n on Γ n n , n \ Γ n , F n → R 2 The ALE map A t n +1 : Ω x 1 + ∆ t ϑ n x 2 + ∆ t ϑ n A t n +1 ( � x 1 , � x 2 ) = ( � 1 , � 2 ) . We define Ω F n +1 = A t n +1 (Ω F n ) and Γ n +1 = A t n +1 (Γ n ) v F , n +1 : Ω F n → R 2 and � p F , n +1 : Ω F We introduce � n → R defined by v F , n +1 ( � x ) = v F , n +1 ( x ) , p F , n +1 ( � x ) = p F , n +1 ( x ) , � � x ∈ Ω F x ) ∈ Ω F ∀ � n , x = A t n +1 ( � n +1

Time discretization of the fluid equations v F , n +1 : Ω F n → R 2 and � p F , n +1 : Ω F Find � n → R such that: �� � v F , n +1 − v F , n �� v F , n − ϑ n � � ρ F v F , n +1 + · ∇ � ∆ t � v F , n +1 � p F , n +1 = � − 2 µ F ∇ · ǫ f F , n +1 , in Ω F � + ∇ � n v F , n +1 = 0 , in Ω F ∇ · � n p F , n +1 ) · n F = h n +1 σ F ( � v F , n +1 , � , on Σ 1 in p F , n +1 ) · n F = h n +1 v F , n +1 , � σ F ( � out , on Σ 3 v F , n +1 = 0 , on Σ 2 � v F ( X , 0) = v 0 ( X ) , in Ω F 0 .

Time discretization of the structure equations 0 → R 2 such that: Find u S , n +1 , ˙ u S , n +1 , ¨ u S , n +1 : Ω S u S , n +1 − ∇ · σ S ( u S , n +1 ) = f S , n +1 , ρ S ¨ in Ω S 0 u S , n +1 = 0 , on Γ D σ S ( u S , n +1 ) n S = 0 , on Γ N u S ( X , 0) = u 0 ( X ) , in Ω S 0 � u S , n +1 � u S , n +1 = ˙ u S , n + ∆ t u S , n + δ ¨ ˙ (1 − δ )¨ u S , n + (∆ t ) 2 �� 1 � u S , n +1 � u S , n +1 = u S , n + ∆ t ˙ u S , n + θ ¨ 2 − θ ¨ . For δ = 1 2, the Newmark scheme is of second order in time.

Interface conditions We define T = A t n ◦ A t n − 1 · · · ◦ A t 1 . We have T (Γ 0 ) = Γ n . v F , n +1 ◦ T u S , n +1 , on Γ 0 × (0 , T ] � = ˙ ( σ F ( � v F , n +1 , � p F , n +1 ) n F ) ◦ T − σ S ( u S , n +1 ) n S , on Γ 0 × (0 , T ] = v F , n +1 , � p F , n +1 , u S , n +1 , Remark The global system of unknowns � u S , n +1 is implicit, but the fluid domain is computed u S , n +1 , ¨ ˙ explicitly as the image of Ω F n via the map x + ∆ t ϑ n ( � � x → � x ) . This is the meaning of the term “semi-implicit” of the title.

Weak formulation of the fluid equations � � w F ∈ ( H 1 (Ω F w F = 0 on Σ 2 W F � n )) 2 ; � Q F � n = L 2 (Ω F n = � n ) , v F , n +1 ∈ � p F , n +1 ∈ � W F Q F Find � n and � n such that: � � ρ F ��� v F , n − ϑ n � � v F , n +1 � v F , n +1 ρ F � w F + w F · � · ∇ � · � ∆ t Ω F Ω F � n n � � w F � � v F , n +1 � � w F � p F , n +1 + 2 µ F ǫ − ∇ · � � � : ǫ � Ω F Ω F n n � � σ F n F � w F = L F ( � w F ∈ � w F ) , W F − · � ∀ � n Γ n � v F , n +1 ) = 0 , q ∈ � Q F � q ( ∇ · � ∀ � n Ω F n

Weak formulation of the structure equations. Lagrangian coordinates � � W S = w S ∈ ( H 1 (Ω S 0 )) 2 ; w S = 0 on Γ D . u S , n +1 ∈ W S such that: Find ˙ � 2 ρ S u S , n +1 · w S + 2 θ ∆ t a S (˙ u S , n +1 , w S ) ∆ t ˙ Ω S 0 � ( σ S n S ) · w S = L S ( w S ) , ∀ w S ∈ W S , − Γ 0 where � � � λ S ( ∇ · u )( ∇ · w ) + 2 µ S ǫ ( u ) : ǫ ( w ) a S ( u , w ) = Ω S 0

Weak formulation of the structure equations. Eulerian coordinates � � W S � w ∈ ( H 1 (Ω S n )) 2 ; n = � w = 0 on Γ D � . � 2 ρ S v S , n +1 · � v S , n +1 , � ∆ t � w + 2 θ ∆ t ˜ a S ( � w ) Ω S � n w = ˜ w ∈ � ( σ S n S ) · � W S − L S ( � w ) , ∀ � n , Γ n where � � � λ S ( ∇ · u )( ∇ · w ) + 2 µ S ǫ ( u ) : ǫ ( w ) ˜ a S ( u , w ) = Ω S n

Global moving domain Ω n = Ω F n ∪ Ω S n Global velocity and pressure v n : Ω n → R 2 , p n : Ω n → R �� � � w ∈ ( H 1 (Ω n )) 2 ; � Q n = L 2 (Ω n ) � W n = w = 0 on Γ D ∪ Σ 2 , Characteristic functions related to the fluid domain χ Ω F t : Ω t → R and structure domain χ Ω S t : Ω t → R : � S 1 , on Ω t = t t = 1 − χ Ω S χ Ω S χ Ω F t . 0 , otherwise

Monolithic formulation for the fluid-structure equations p n +1 ) ∈ � W n × � v n +1 , � Find ( � Q n such that: � � n ρ F ��� n � � v n +1 � v n +1 n ρ F � v n − ˜ ∆ t · � w + · ∇ � · � χ Ω F χ Ω F ϑ w Ω n Ω n � � �� v n +1 � p n +1 + n 2 µ F ǫ − χ Ω F n ( ∇ · � w ) � χ Ω F : ǫ ( � w ) Ω n Ω n � 2 ρ S v n +1 · � v n +1 , � + χ Ω S ∆ t � w + 2 θ ∆ t ˜ a S ( � w ) n Ω n = ˜ w ) + ˜ w ∈ � L F ( � L S ( � w ) , ∀ � W n � v n +1 ) = 0 , q ∈ � n � q ( ∇ · � ∀ � χ Ω F Q n . Ω n

Finite element discretization Triangular P 1 + bubble for the velocity, P 1 for the pressure and P 0 for the characteristic functions. The velocity, the pressure as well as the test functions are continuous all over the global domain Ω n . Consequently, the continuity of velocity at the interface is automatically satisfied. The integrals over the interface do not appear explicitly in the global weak form due to the action and reaction principle. If the solution of the monolithic is sufficiently smooth, the continuity of stress at the interface holds in a weak sense.

Linear systems The linear system has not an unique solution, because the pressure p S can take any value.       B T A 0 v L  =    p F   B 0 0 0 p S 0 0 0 0 � p n +1 � We have added the term ǫ � q , then the bellow system has Ω n an unique solution and p S = 0 on Ω S n .       B T A 0 v L     =   ǫ M F p F B 0 0 ǫ M S p S 0 0 0 The matrix A is not symetric due to the convection term. We have used the GMRES algorithm for solving the linear system.

Time advancing schema from n to n + 1 We assume that we know Ω n , v n , u n , ¨ u n . n Step 1 : Compute ˜ ϑ v n +1 and pressure � p n +1 Step 2 : Solve the linear system and get � Step 3 : Compute the displacement and the acceleration u n + (∆ t ) 2 � 1 � u n + ∆ t (1 − 2 θ ) v n + 2 θ ∆ t � u n +1 v n +1 � = 2 − 2 θ ¨ � v n +1 − v n � 2 n +1 � u n ¨ = � − ¨ u ∆ t Step 4 : We define the map T n : Ω n → R 2 by: n ( � x + ∆ t ˜ u n +1 ( � x ) − u n ( � T n ( � x ) = � x ) χ Ω F n ( � x ) + ( � x )) χ Ω S n ( � x ) ϑ Step 5 : We set Ω n +1 = T n (Ω n ). We define v n +1 : Ω n +1 → R 2 and p n +1 : Ω n +1 → R 2 by: v n +1 ( x ) = � v n +1 ( � x ) , p n +1 ( x ) = � p n +1 ( � x ) , ∀ � x ∈ Ω n and x = T n ( � x ) .

Global fluid-structure mesh (top), the structure and fluid meshes (bottom) The global moving mesh is compatible with the interface: a triangle of the global mesh belongs either to the fluid region or to the structure region.