Modeling and numerical simulations of fish like swimming Michel - PowerPoint PPT Presentation

Modeling and numerical simulations of fish like swimming Michel Bergmann, Angelo Iollo INRIA Bordeaux Sud-Ouest, ´ equipe MC2 Institut de Math´ ematiques Appliqu´ ees de Bordeaux 33405 TALENCE cedex, France Workshop Maratea, may 13 2010 – p. 1

Context and objectives ◮ Context : ANR CARPEiNETER Cartesian grids, penalization and level set for the simulation and optimisation of complex flows ◮ Objectives : ֒ → Model and simulate moving bodies S (translation, rotation, deformation, ..) ֒ → Couple Fluid and Structures ֒ → Cartesian meshes Avoid remeshing ֒ → Penalization of the equations To take into account the bodies ֒ → Level Set To track interfaces (fluid/fluid, fluid/structures) Workshop Maratea, may 13 2010 – p. 2

Outline Flow modeling Numerical approach Method: discretization / body motion Validation Applications: 2D fish swimming Parametrization Classification: BCF On the power spent to swim Maneuvers and turns Fish school (3 fishes) 3D locomotion Conclusions and future works Workshop Maratea, may 13 2010 – p. 3

Flow modeling ∂ Ω ∂ Ω 1 Ω 1 u 1 u 2 ∂ Ω 2 Ω f Ω 2 Ω i : Domain "body" i Ω f : Domain "fluid" Ω = Ω ∪ Ω i : Entire domain Workshop Maratea, may 13 2010 – p. 4

Flow modeling ◮ Classical model: Navier-Stokes equations (incompressible): � ∂ u � ρ ∂t + ( u · ∇ ) u = −∇ p + µ ∆ u + ρ g dans Ω f , (1a) ∇ · u = 0 dans Ω f , (1b) u = u i sur ∂ Ω i (1c) u = u f sur ∂ Ω (1d) Numerical resolution Need of meshes that fit the body geometries ֒ → Costly remeshing for moving and deformable bodies!! Workshop Maratea, may 13 2010 – p. 5

Flow modeling ◮ Penalization model: penalized Navier-Stokes equations (incompressible): � ∂ u � N s � ρ ∂t + ( u · ∇ ) u = −∇ p + µ ∆ u + ρ g + λρ χ i ( u i − u ) dans Ω , (2a) 1=1 ∇ · u = 0 Ω , dans (2b) u = u f sur ∂ Ω . (2c) λ ≫ 1 penalization factor → Solution eqs (2) tends to solution eqs (1) w.r.t. ε = 1 /λ → 0 . χ i characteristic function: χ i ( x ) = 1 if x ∈ Ω i , (3a) χ i ( x ) = 0 else if . (3b) Numerical resolution No need of meshes that fit the body geometries ֒ → Cartesian meshes Workshop Maratea, may 13 2010 – p. 6

Flow modeling ◮ Transport of the characteristic function for moving bodies ∂χ i ∂t + ( u i · ∇ ) χ i = 0 . (4) Other choice: χ i = H ( φ i ) where H is Heaviside function and φ i the signed distance function ( φ i ( x ) > 0 if x ∈ Ω i , φ i ( x ) = 0 si x ∈ ∂ Ω i , φ i ( x ) < 0 else if). ∂φ i ∂t + ( u i · ∇ ) φ i = 0 . (5) ◮ Density � � � N s � N s � ρ = ρ f 1 − χ i + ρ i χ i . (6) i =1 i =1 Workshop Maratea, may 13 2010 – p. 7

Flow modeling ◮ Dimensionless equations with U ∞ , D , ρ f , Re = ρU ∞ D : µ � N s ∂ u ∂t + ( u · ∇ ) u = −∇ p + 1 Re ∆ u + g + λ χ i ( u i − u ) dans Ω , (7a) 1=1 ∇ · u = 0 dans Ω , (7b) u = u f sur ∂ Ω (7c) ◮ Body velocity u i : u i = u i + Û u i + � u i (8) with: u i translation velocity Û u i rotation velocity � u deformation velocity (imposed for the swim) Workshop Maratea, may 13 2010 – p. 8

Numerical approach | Method ◮ Space: Cartesian meshes, collocation with compact "non oscillating" scheme, Centered FD 2nd order and upwind 3rd order for convective terms ◮ Time: 1 st order explicit euler, implicit penalization (larger λ ) u ( n +1) − u ( n ) + ( u ( n ) · ∇ ) u ( n ) = − ∇ p ( n +1) + 1 Re ∆ u ( n +1) + g ∆ t � N s χ i ( n +1) ( u i ( n +1) − u ( n +1) ) , + λ 1=1 ∇ · u ( n +1) =0 ⇒ Problems ֒ → Pressure uncoupled → The function χ i ( n +1) and velocity u i ( n +1) are not known ֒ ⇒ Solutions ֒ → Chorin scheme (predictor/corrector) ֒ → Fractional step method (2 steps) Workshop Maratea, may 13 2010 – p. 9

Numerical approach | Method ◮ Fractional steps method u ( n +1) − u ( n ) + ( u ( n ) · ∇ ) u ( n ) = − ∇ p ( ∗ ) + 1 Re ∆ u ( n +1) + g ∆ t + � ∇ p ( ∗ ) − ∇ p ( n +1) � � N s χ i ( n +1) ( u i ( n +1) − u ( n +1) ) , + λ 1=1 ∇ · u ( n +1) =0 u ( n +1) = f ( u ( n +1) , p ( n +1) ) i Step 1: ⇒ u ( ∗ ) , p ( ∗ ) Step 2 : ⇒ � u ( n +1) , � p ( n +1) Step 3 : ⇒ u ( n +1) = � f ( � u ( n +1) , � p ( n +1) ) i Step 4 : ⇒ u ( n +1) , p ( n +1) Workshop Maratea, may 13 2010 – p. 10

Numerical approach | Method ◮ Step 1 : prediction u ( ∗ ) − u ( n ) + ( u ( n ) · ∇ ) u ( n ) = −∇ p ( ∗ ) + 1 Re ∆ u ( ∗ ) + g ∆ t ◮ Step 2 : correction u ( n +1) − u ( ∗ ) � = ∇ p ( ∗ ) − ∇ p ( n +1) ∆ t u ( n +1) =0 ∇ · � with ψ = ∇ p ( ∗ ) − ∇ p ( n +1) , on a ∆ ψ = ∇ · u ( ∗ ) u n +1 = � u ∗ − ∇ ψ � p ∗ + ψ p n +1 = � � ∆ t Workshop Maratea, may 13 2010 – p. 11

Numerical approach | Method ◮ Etape 3 : body motion Compute forces F i and torques M i m d u i = F i + m g , translation velocity , m mass u i (14a) d t d J Ω i = M i , angular velocity , J inertia matrix Ω i (14b) d t Rotation velocity Û u i = Ω i × r G with r G = x − x G ( x G center of mass). Re ( ∇ u + ∇ u T ) et n outward normal unit vector at s i : 1 Stress tensor T ( u , p ) = − p I + � F i = − T ( u , p ) n d x , (15a) ∂ Ω i � M i = − T ( u , p ) n × r G d x . (15b) ∂ Ω i Evaluation of forces and torques Cartesian mesh: no direct acces to ∂ Ω i ֒ → Not easy evaluation .... Workshop Maratea, may 13 2010 – p. 12

Numerical approach | Method Definition : Arbitrarily domain Ω f i ( t ) surrounding body i . Forces: � � F i = − d u d V + ( T + ( u − u i ) ⊗ u ) n d S d t Ω fi ( t ) ∂ Ω fi ( t ) � (16a) + (( u − u i ) ⊗ u ) n d S. ∂ Ω i ( t ) Torques: � � M i = − d u × r G d V + ( T + ( u − u i ) ⊗ u ) n × r G d S d t Ω fi ( t ) ∂ Ω fi ( t ) � (16b) + (( u − u i ) ⊗ u ) n × r G d S. ∂ Ω i ( t ) Evaluation of forces and torques The term onto ∂ Ω i vanishes in our case (no transpiration) ֒ → Easy evaluation! Workshop Maratea, may 13 2010 – p. 13

Numerical approach | Method ◮ Step 4 : Update velocity using implicit penalization � N s u ( n +1) − � u ( n +1) χ i ( n +1) ( u i ( n +1) − u ( n +1) ) = λ ∆ t 1=1 ◮ Summary: ⊲ Solve Navier-Stokes equation without penalization ⇒ � u ( n +1) , � p ( n +1) ⊲ Compute body motion ⇒ u ( n +1) , χ ( n +1) i i ⊲ Correct solution with penalization ⇒ u ( n +1) , p ( n +1) ◮ Remark: ⊲ Step 4 can be implemented in step 1using explicit body velocity (time order is 1). Workshop Maratea, may 13 2010 – p. 14

Numerical approach | Validation ◮ Improvement of the penalization order ֒ → Test case: 2D Green-Taylor vortex with analytical solution ( 0 ≤ x, y ≤ π , Re = 100 ) € � u ( t, x ) = sin( x ) cos( y ) exp( − 2 t/Re ) , v ( t, x ) = − cos( x ) sin( y ) exp( − 2 t/Re ) , u ( T f , x ) − u ( T f , x )) 2 d x. ( � E = p ( t, x ) = 1 4 (cos(2 x ) + cos(2 y )) exp( − 4 t/Re ) . Ω 1 χ = 0 0.9 Fluid 0.8 0.7 0.6 χ = 1 0.5 " Body " 0.4 0.3 0.2 0.1 0 ֒ → "Non intrusive" body ⇒ penalization velocity depends on space and time Workshop Maratea, may 13 2010 – p. 15

Numerical approach | Validation χ = 0 χ = 0 " Fluid " " Fluid " 1 - No penalization ֒ → use analytical boundary conditions → Numerical scheme order, (∆ x ) 2 ⇒ 2 nd order x i − 1 x i +1 x i − 2 x i 0 No penalization → order 2 1.2E-07 -2 1.1E-07 -4 1.0E-07 9.1E-08 -6 7.9E-08 log E -8 6.8E-08 5.7E-08 -10 4.5E-08 -12 3.4E-08 2.3E-08 -14 1.1E-08 -16 -4.5 -4 -3.5 -3 -2.5 log ∆ x Workshop Maratea, may 13 2010 – p. 16

Numerical approach | Validation χ = 1 χ = 0 " Body " " Fluid " 2 - "Exact" penalization: ֒ → use analytical penalization values u n i = u n i ⇒ 2 nd order x i − 1 x i +1 x i − 2 x i 0 No penalization → order 2 "Exact" pen. → order 2 1.2E-07 -2 1.1E-07 -4 1.0E-07 9.0E-08 -6 7.8E-08 log E -8 6.7E-08 5.6E-08 -10 4.5E-08 -12 3.4E-08 2.2E-08 -14 1.1E-08 -16 -4.5 -4 -3.5 -3 -2.5 -2 log ∆ x Workshop Maratea, may 13 2010 – p. 16

Numerical approach | Validation χ = 1 χ = 0 " Body " " Fluid " 3 - "Standard" penalization: ֒ → use only boundary velocity u n i = u n φ =0 ⇒ 1 nd order x i − 1 x i +1 x i − 2 x i 0 No penalization → order 2 "Exact" pen. → order 2 4.7E-03 -2 "Classic" pen. → order 1 4.3E-03 -4 3.8E-03 3.4E-03 -6 3.0E-03 log E -8 2.6E-03 2.1E-03 -10 1.7E-03 -12 1.3E-03 8.5E-04 -14 4.3E-04 -16 -4.5 -4 -3.5 -3 -2.5 log ∆ x Workshop Maratea, may 13 2010 – p. 16

Numerical approach | Validation χ = 1 χ = 0 " Body " " Fluid " 4 - "Improved" penalization: ֒ → use Level Set informations φ =0 − φ i ( ∂ u i /∂ n ) n − 1 u n i = u n ⇒ 2 nd order x i − 1 x i +1 x i − 2 x i 0 No penalization → order 2 "Exact" pen. → order 2 1.4E-05 -2 "Classic" pen. → order 1 1.3E-05 -4 "Improved" pen. → order 2 1.2E-05 1.0E-05 -6 9.0E-06 log E -8 7.8E-06 6.5E-06 -10 5.2E-06 -12 3.9E-06 2.6E-06 -14 1.3E-06 -16 -4.5 -4 -3.5 -3 -2.5 -2 log ∆ x Workshop Maratea, may 13 2010 – p. 16