Numerical Simulation of Drosophila Flight Based on Arbitrary - - PowerPoint PPT Presentation

Numerical Simulation of Drosophila Flight Based on Arbitrary Lagrangian Eulerian (ALE) Method

Belkis ERZINCANLI & Mehmet SAHIN 65th Annual Meeting of the APS Division of Fluid Dynamics Aeronautical Engineering Department, November 18-20, 2012 – San Diego,CA. Faculty of Aeronautics and Astronautics, Istanbul Technical University, 34469, Maslak/Isatanbul, TURKEY

ITU/FAA

Faculty of Aeronautics and Astronautics

Contents

• Motivations
• ALE Formulation
• Governing equations
• Discrete Geometric Conservation Law (GCL)
• 3D Numerical discretization
• Iterative solver
• Numerical Validations
• Kovasznay flow
• The flow past an oscillating circular cylinder in a channel
• The flow induced due to a horizontally oscillating sphere
• The Numerical Simulation of Flow Field Around Drosophila
• Conclusions and Future Work

APS DFD 2012 2

ITU/FAA

Faculty of Aeronautics and Astronautics

The Fruit Fly, Drosophila

APS DFD 2012 3

ITU/FAA

Faculty of Aeronautics and Astronautics

Debat et al., Evolution, 60(12), 2006, pp. 2529–2538

Insect Flight Aerodynamics

APS DFD 2012 4

ITU/FAA

Faculty of Aeronautics and Astronautics

Clap & Fling A comparison of 2-D linear translation vs 3-D flapping translation Stable attachment of the leading edge vortex

Sane, S.P.,”The Aerodynamics of Insect Flight”, The Journal of Experimental Biology 206,4191-4208, (2003).

Numerical Aprroaches Used for Moving Boundary Flow Problems

• Immersed Boundary Method (IBM): This approach solves Navier-Stokes

equations on a Cartesian mesh with a moving immersed body (Miller & Peskin, 2004; Gilmanov & Sotiropoulos, 2005; Mittal et al., 2006).

• Fictious Domain Method: This approach uses Eulerian mesh and the body

motion is represented using the distributed Lagrangian multiplier method.

• Overset (Chimera) Grid Method: This method uses several body fitted

meshes inside the fluid domain and an interpolation algorithm is used between the neighbouring meshes (Liu and Aono, 2009).

• Arbitrary Lagrangian-Eulerian (ALE) Method: The mesh follows the

interface between the fluid and solid boundary and the governing equations are discretized on a moving/deforming mesh (Johnson, 2006; Ramamurti & Sandberg, 2002).

APS DFD 2012 5

ITU/FAA

Faculty of Aeronautics and Astronautics

A Fully Coupled Parallel Unstructured Finite Volume Method based on ALE Formulation

• In the Arbitrary Lagrangian-Eulerian (ALE) method, the interface between

the solid and fluid is sharply represented.

• The mesh movement has to satisfy a special condition called Geometric

Conservation Law (GCL) in order to maintain the accuracy and the stability

• f the time integration scheme.
• An unstructured ALE method based on the side-centered finite volume

technique has been developed for parallel large-scale moving boundary problems.

• The present arrangement of the primitive variables leads to a stable

numerical scheme and it does not require any ad-hoc modifications in order to enhance the pressure-velocity coupling.

• The most appealing feature of present primitive variable arrangement is the

availability of very efficient multigrid solvers.

• A fully coupled implicit iterative solver is used for large-scale calculations

since the CFL restriction may lead to inadmissible time steps on deforming meshes and does not introduce any splitting error when the advection- diffusion operator is treated impicitly.

APS DFD 2012 6

ITU/FAA

Faculty of Aeronautics and Astronautics

Mathematical and Numerical Formulation

The incompressible Navier-Stokes equations that govern the viscous fluid flow of an arbitrary moving Eulerian Lagrangian control volume Ω(t) with boundary ∂Ω(t) can be written in Cartesian coordinate system in dimensionless form as follows: The continuity equation: The momentum equation:

Three-dimensional dual volume

APS DFD 2012 7

ITU/FAA

Faculty of Aeronautics and Astronautics

Grid velocity

Numerical Discretization

Side-centered finite volume method

• A stable numerical scheme with exact mass conservation
• No ad-hoc modifications for pressure-velocity coupling
• Very efficient multigrid solvers are available

(a) Two-dimensional dual volume (b) Three-dimensional dual volume

APS DFD 2012 8

ITU/FAA

Faculty of Aeronautics and Astronautics

Numerical Discretization Cont’d

The contributions for the x-momentum equation for the right element: The time derivation is the volume of the pyramid between points

APS DFD 2012 9

ITU/FAA

Faculty of Aeronautics and Astronautics

Numerical Discretization Cont’d

The convective term due to fluid motion are the area vectors of the dual volume triangular

• surfaces. are the velocity vectors defined at the mid-

point of each dual volume surface area.

APS DFD 2012 10

ITU/FAA

Faculty of Aeronautics and Astronautics

Numerical Discretization Cont’d

The convective term due to mesh motion First-order backward difference is used for grid velocity

APS DFD 2012 11

ITU/FAA

Faculty of Aeronautics and Astronautics

Geometric conservation law satisfied.

Numerical Discretization Cont’d

The pressure term 2nd order Taylor series expansion is used for the pressure values at .

APS DFD 2012 12

ITU/FAA

Faculty of Aeronautics and Astronautics

Numerical Discretization Cont’d

The viscous term Gauss-Green theorem is used for the gradient terms:

APS DFD 2012 13

ITU/FAA

Faculty of Aeronautics and Astronautics

Numerical Discretization Cont’d

The continuity equation is integrated within each hexahedral elements and evaluated using the mid-point rule on each of the element faces where is the hexahedral element surface area vector and u, v and w are the velocity vector components defined at the mid-point of each hexahedral element face. The discretization of above equations leads to a saddle point problem of the form: where, are the convection diffusion operators, is the pressure gradient operator and is the divergence operator.

APS DFD 2012 14

ITU/FAA

Faculty of Aeronautics and Astronautics

Iterative Method

The preconditioner matrix is Where .For the inverse of the scaled Laplacian S, we use two-cycle AMG solver provided by the HYPRE library, a high performance preconditioning package developed at Lawrence Livermore National Laboratory, which we access through the PETSC library.

APS DFD 2012 15

ITU/FAA

Faculty of Aeronautics and Astronautics

TEST CASE I: 2D Kovasznay Flow

APS DFD 2012 16

ITU/FAA

Faculty of Aeronautics and Astronautics

The spatial convergence of the error with mesh refinement .

TEST CASE II: An Oscillating Circular Cylinder in a Channel

• The location of the cylinder center
• Amplitude and frequency of the oscillation
• Time step
• Reynolds number
• 70,667 quadrilateral elements
• 71,349 vertices
• 354,699 total DOF

APS DFD 2012 17

ITU/FAA

Faculty of Aeronautics and Astronautics

TEST CASE II: An Oscillating Circular Cylinder in a Channel

The computed u-velocity vector component contours with streamtraces.

APS DFD 2012 18

ITU/FAA

Faculty of Aeronautics and Astronautics

t=20.0s t=21.0s t=22.0s t=23.0s

TEST CASE II: An Oscillating Circular Cylinder in a Channel

The comparion of cd and cl plots with results from Wan and Turek (JCP, 2007).

APS DFD 2012 19

ITU/FAA

Faculty of Aeronautics and Astronautics

ITU/FAA

Faculty of Aeronautics and Astronautics

TEST CASE III: An Oscillating Sphere in a Cubic Cavity

• The location of the sphere
• Time step
• Reynolds number (Gilmanov & Sotiropoulos , 2005)

The Reynolds number is based on the sphere diameter and the maximum sphere velocity.

• 473,856 hexahedral elements
• 492,596 vertices
• 4,794,264 total DOF

APS DFD 2012 20

ITU/FAA

Faculty of Aeronautics and Astronautics

TEST CASE III: An Oscillating Sphere in a Cubic Cavity

• The location of the sphere x(t) = h [1- cos(2฀t)] h = 0.125D

APS DFD 2012 21

The instantaneous streamtraces with u-velocity component contours at different phases over a cycle for a rigid sphere oscillating in a cubic cavity at Re = 20: t=0, t=T/4, t=T/2 and t=3T/4. The time variation of drag coefficient for a rigid sphere oscillating in a cubic cavity at Re = 20.

The Numerical Simulation of Flow Field Around Drosophila

APS DFD 2012 22

ITU/FAA

Faculty of Aeronautics and Astronautics

Robofly experiment of M. Dickinson (CalTech / ICB), Science (1999).

The Numerical Simulation of Flow Field Around Drosophila

APS DFD 2012 23

ITU/FAA

Faculty of Aeronautics and Astronautics

Downstroke motion at times t=0, 0.1, 0.2, 0.3, 0.4 and 0.5. Upstroke motion at times t=0.5, 0.6, 0.7, 0.8, 0.9 and 1. The approximated kinematics of the Drosophila wing motion:

The Numerical Simulation of Flow Field Around Drosophila

APS DFD 2012 24

ITU/FAA

Faculty of Aeronautics and Astronautics

The initial computational fine mesh [a] and the control points mesh used for RBF based mesh deformation algorithm [b] for the Drosophila wing. The computational mesh consists of 1,300,358 vertices and 1,276,666 hexahedral elements (12,837,448 DOF).

[a] [b]

Radial Basis Function (RBF) Method for Mesh Deformation

where

APS DFD 2012 25

ITU/FAA

Faculty of Aeronautics and Astronautics

is the weight of center is the radial basis function is a low-degree polynomial

M(x)

is the Eucledian norm (Thin Plate Spline)

RBF - The Linear System Ax = b

where

APS DFD 2012 26

ITU/FAA

Faculty of Aeronautics and Astronautics

The resulting algebraic equations for the displacement in the x-direction: The constraints due to rigid body motion where is the radial basis function.

APS DFD 2012 27

ITU/FAA

Faculty of Aeronautics and Astronautics

The initial mesh [a] and deformed meshes at y=0 plane using indirect Radial Basis Function (RBF) method with 1 iteration (direct method) [b], 2 iterations [c] and 3 iterations [d]. [a] [b] [c] [d]

The Mesh Deformation with RBF Method

APS DFD 2012 28

ITU/FAA

Faculty of Aeronautics and Astronautics

The instantaneous downstroke wake structures (lambda2-criterion) around the Drosophila wing at several different time levels: t = 0.0T [a], t = 0.1T [b], t = 0.2T [c], t = 0.3T [d], t = 0.4T [e] and t = 0.5T [f ].

[a] [b] [c] [d] [e] [f]

The Numerical Simulation of Flow Field Around Drosophila

APS DFD 2012 29

ITU/FAA

Faculty of Aeronautics and Astronautics

The instantaneous upstroke wake structures (lambda2-criterion) around the Drosophila wing at several different time levels: t = 0.5T [a], t = 0.6T [b], t = 0.7T [c], t = 0.8T [d], t = 0.9T [e] and t = 1.0T [f ].

[a] [b] [c] [d] [e] [f]

The Numerical Simulation of Flow Field Around Drosophila

The Numerical Simulation of Flow Field Around Drosophila

APS DFD 2012 30

ITU/FAA

Faculty of Aeronautics and Astronautics

The instantaneous wake structures (Q-criterion) around the Drosophila wing.

The Numerical Simulation of Flow Field Around Drosophila

APS DFD 2012 31

ITU/FAA

Faculty of Aeronautics and Astronautics

The computed instantaneous pressure contours during the downstroke motion at t=0, t=0.1, t=0.2, t=0.3, t=0.4 and t=0.5 [a] and the instantaneous v-velocity component isosurfaces with the streamtraces showing the stable leading edge vortex [b] . [b] [a]

The Numerical Simulation of Flow Field Around Drosophila

APS DFD 2012 32

ITU/FAA

Faculty of Aeronautics and Astronautics

The computed total force coefficients for a pair of Drosophila wings [a] and the comparison of the lift coefficient with the experimental (Dickinson et. al, 1999) and numerical results (Kweon and Choi, 2010) [b]. [a] [b]

Conclusions and Future Work

• A parallel fully-coupled unstructured ALE method based on the side-

centered unstructured finite volume has been developed and successfuly tested for 2D and 3D benchmark problems.

• Boddy fitted unstructured meshes are allowed on the solid surface.
• The continuity equation is satisfied within each element at machine

precission.

• A special attention is given to satisfy the Geometric Conservation Law.
• The current mesh deformation is achieved by using the radial basis function

(RBF) approach. Large mesh deformations with high mesh quality and periodicity can be achieved.

• A fully-coupled implicit iterative solver is used at each time level.
• The numerical method is applied to the numerical solution of the flow field

around a pair of flapping Drosophila wing in hover.

APS DFD 2012 33

ITU/FAA

Faculty of Aeronautics and Astronautics

Acknowledgement

The authors gratefully acknowledge the financial support from Scientific and Technical Research Council of Turkey (TUBITAK ) under project number

• 111M332. The authors also would like to acknowledge the use of the Chimera

machine at the Faculty of Aeronautics and Astronautics at ITU, the computing resources provided by the National Center for High Performance Computing

• f Turkey (UYBHM) under grant number 10752009 and the computing

facilities at TUBITAK ULAKBIM, High Performance and Grid Computing Center. The authors also would like to thank Michael Dickinson and his graduate student Michael Elzinga for providing the experimental data.

APS DFD 2012 34

ITU/FAA

Faculty of Aeronautics and Astronautics

THANK YOU FOR YOUR ATTENTION!

APS DFD 2012 35

ITU/FAA

Faculty of Aeronautics and Astronautics