An arbitrary Lagrangian–Eulerian formulation for solving moving boundary problems with large displacements and rotations

Journal of Computational Physics 255 (2013) 660–679

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

An arbitrary Lagrangian–Eulerian formulation for solving moving boundary problems with large displacements and rotations Belkis Erzincanli, Mehmet Sahin ∗ Faculty of Aeronautics and Astronautics, Astronautical Engineering Department, Istanbul Technical University, Maslak, Istanbul 34469, Turkey

a r t i c l e

i n f o

Article history: Received 25 February 2013 Received in revised form 10 July 2013 Accepted 19 August 2013 Available online 2 September 2013 Keywords: Unstructured finite volume ALE methods Insect flight Incompressible viscous flow Coupled iterative solvers Large-scale computations

a b s t r a c t An Arbitrary Lagrangian–Eulerian (ALE) formulation based on the unstructured finite volume method is proposed for solving moving boundary problems with large displacements and rotations. The numerical method is based on the side-centered arrangement of the primitive variables that does not require any ad-hoc modifications in order to enhance pressure coupling. The continuity equation is satisfied within each element at machine precision and the summation of the continuity equations can be exactly reduced to the domain boundary, which is important for the global mass conservation. A special attention is given to construct an ALE algorithm obeying the discrete geometric conservation law (DGCL). The mesh deformation algorithm is based on the indirect Radial Basis Function (RBF) algorithm at each time level while avoiding remeshing in order to enhance numerical robustness. For the parallel solution of resulting large-scale algebraic equations in a fully coupled form, a matrix factorization is introduced similar to that of the projection method for the whole system and the parallel algebraic multigrid solver BoomerAMG is used for the scaled discrete Laplacian provided by the HYPRE library which we access through the PETSc library. The present numerical algorithm is initially validated for the decaying Taylor–Green vortex flow, the flow past an oscillating circular cylinder in a channel and the flow induced by an oscillating sphere in a cubic cavity. Then the numerical algorithm is applied to the numerical simulation of flow field around a pair of flapping Drosophila wings in hover flight. The time variation of the Eulerian coherent structures in the near wake is shown along with the aerodynamic loads. © 2013 Elsevier Inc. All rights reserved.

1. Introduction Many fluid dynamics problems involve moving/deforming boundaries that represent an interface between different fluids or between fluids and solids. Examples include a wide variety of fluid–structure interaction problems, a large class of free-surface problems, fluid–particle interactions, swimming/flying animals, etc. One of the most used techniques for the numerical simulation of moving boundary problems is the so-called Arbitrary Lagrangian–Eulerian (ALE) formulation [23]. The advantage of the technique is that the interface is sharply defined and the implementation of interface boundary conditions is straightforward. In the ALE method, the mesh follows the interface between the fluid and solid boundary and the governing equations are discretized on unstructured moving meshes. This differs from the standard Eulerian formulation in a way that the mesh movement has to fulfill special conditions in order to maintain the accuracy and the stability of the time

*

Corresponding author. E-mail address: [email protected] (M. Sahin).

0021-9991/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.08.038

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

661

integration scheme. This condition is satisfied by the enforcement of the so-called geometric conservation law (GCL) [53]. In the past decade there has been a considerable research effort into the development of stable and accurate time integration schemes for ALE schemes and identification of the role of the discrete geometric conservation law (DGCL) in terms of accuracy and stability [31,19]. The use of the ALE formulations for moving boundary problems with large displacements and rotations is rather challenging, in particular the moving boundary problem of rigid and flexible flapping wings due the flapping wing kinematics undergoing translation and rotation at large amplitudes which may lead to inadmissible elements. Ramamurti and Sandberg [40] applied an arbitrary Lagrangian Eulerian based Galerkin finite element procedure with linear tetrahedral elements to solve unsteady flow past a three-dimensional Drosophila wing. The approach allows the near-wing grid to move as the wing does, but remeshing is still required to eliminate badly distorted elements. Johnson [26] developed a dynamic mesh algorithm for complex moving boundary problems and tried to reproduce the experimental work of Dickinson et al. [17]. Bos [10] used the ALE based open-source toolbox OpenFOAM in order to investigate the effects of different wing kinematics. The objective of this study is to incorporate an arbitrary Lagrangian–Eulerian approach with the side-centered unstructured finite volume method given in [45] to establish a novel algorithm for the large-scale simulation of moving boundary problems in a fully coupled form. The numerical method is based on the side-centered finite volume method where the velocity vector components are defined at the mid-point of each cell face, while the pressure term is defined at element centroids. 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. This approach was initially used by Maliska and Raithby [35] to calculate the contravariant velocity components that enter the mass conservation constraint on boundaryfitted meshes and later by Hwang [24] and Rida et al. [42] for the solution of the incompressible Navier–Stokes equations on unstructured triangular meshes. Hwang [24] pointed out several important computational merits for the aforementioned grid arrangement. Rida et al. [42] called this scheme side-centered finite volume method and the authors reported superior convergence properties compared to the semi-staggered approach. The most appealing feature of this approach is the availability of efficient multigrid solvers which are sufficiently robust even on non-uniform and highly anisotropic meshes. Although the fully staggered approach with multigrid method also leads to a very robust numerical algorithm [54], obtaining the velocity components on unstructured staggered grids is not straightforward as well as the computation of inter grid transfer operators in multigrid. The use of all the velocity vector components significantly simplifies the numerical discretization of the governing equations on unstructured grids as well as the implementation of physical boundary conditions. The present arrangement of the primitive variables can be applied to any nonoverlapping convex polygon which is very important for the treatment of more complex configurations. As far as the authors’ knowledge goes, the present arrangement of the primitive variables is not considered within the ALE framework. Although the Immersed Boundary Method (IBM) [38,36] has been used extensively for the simulation of flow fields around moving/deforming bodies with complex geometrical shapes, the implementation of the physical boundary conditions on surfaces not aligned with the mesh is still a challenging task, in particular, at high Reynolds numbers as well as obtaining sufficiently smooth/accurate forces next to the solid body. In the case of the moving/deforming bodies with large amplitudes (as in flapping flight), the applications of the IBM requires either high mesh resolutions in large percentages of the fluid computational domain, or adaptive mesh refinement, in order to properly capture the viscous effects. An ALE algorithm requires a robust technique to model the motion and deformation of the underlying fluid mesh. In the literature, several mesh deformation algorithms have been proposed to compute the displacement of the internal points as the boundaries of a computational domain translate, rotate and deform in order to maintain mesh quality and validity. These approaches include the spring analogy [6], the elastic medium analogy [27], the edge swapping algorithm [16], the radial basis function (RBF) interpolation algorithm [9] and the remeshing algorithm [28]. The RBF interpolation algorithm leads to high quality meshes when domain boundaries exhibit large translations and rotations [9]. However, the RBF approach is too costly with its most straightforward implementation for large three-dimensional problems. An approximation algorithm for the RBF mesh deformation has been suggested by Rendall and Allen [41] in which the RBF is applied using a coarsened subset of surface mesh. A greedy algorithm is used to add points from the mesh that have the largest error. In the present work the RBF algorithm is implemented using a non-nested coarser mesh. However, the control points are not exactly located on the domain boundary as found in previous studies, but next to the boundary. Therefore, the mesh points between the RBF control points and the domain boundary are moved using the rigid body motion that will help to preserve the quality of boundary layer meshes. The rest of the points are deformed using the RBF interpolation. The present modification will ensure that the total volume of computation domain is conserved at machine precision. This is very important for the use of div-stable discretization of the incompressible Navier–Stokes with all Dirichlet boundary conditions where the total mass should be kept constant at all times. In addition, the indirect implementation of the RBF mesh deformation algorithm is used. In this approach, the wing rotation motion is split into several steps and the new grid nodes are computed starting from the initial mesh. This approach also ensures large mesh deformations while keeping mesh periodicity and quality. In the case of the deforming body, the control points can be moved using an algebraic equation based on the minimum distance function from the deforming solid boundary and the displacement vector at the nearest solid surface node. In a similar manner, the location of the grid points between the RBF control points and the solid domain boundary can also be computed. However, a special attention should be given to the deformation on the solid boundary in order to satisfy the global mass conservation in case all Dirichlet boundary conditions are employed. This may be achieved by solving the incompressible elasticity equations within the solid domain.

662

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

In the present work, we investigate the applicability of the fully coupled iterative solution of the momentum and continuity equations. The main reason is that employed mesh deformation algorithms may lead to extremely small elements for large deformations which significantly limit the allowable time step for decoupled approaches. In addition, the discretized governing equations may need to be coupled to the equation of motion of a body [47]. Multigrid techniques [21,56] are known to be the most efficient numerical techniques for solving large-scale problems that arise in numerical simulations of physical phenomena because of their computational costs and memory requirements that scale linearly with the degrees of freedom. The present two-level iterative solver is based on the multiplicative non-nested multigrid method with one V-cycle and it is described in [45] in detail. The combination of the present numerical method with the geometric nonnested multilevel preconditioner for the three-dimensional Stokes problem has enabled us to simulate the viscoelastic wake instabilities [46] for the first time. The implementation of the preconditioned Krylov subspace algorithm, matrix–matrix multiplication and the multilevel preconditioner were carried out using the PETSc [4] software package developed at the Argonne National Laboratories. However, the fully coupled multigrid techniques are not suitable for the time-dependent calculation of the incompressible Navier–Stokes equations with small time steps since the advection–diffusion operator is highly diagonally dominant and well conditioned. Therefore, a matrix factorization is introduced similar to that of the projection method [13] for the whole coupled system and we use the BoomerAMG solver for the scaled discrete Laplacian provided by the HYPRE [18] library, a high performance preconditioning package developed at Lawrence Livermore National Laboratory, which we access through the PETSc library. The computational domain is decomposed into a set of sub-domains or partitions using the METIS library [29]. The remainder of this paper is organized as follows: Section 2 provides some detail on the present ALE method with efficient iterative solvers and mesh deformation algorithm. In Section 3 the proposed method is initially validated for the decaying Taylor–Green vortex flow, the flow past an oscillating circular cylinder in a channel and the flow induced by an oscillating sphere in a cubic cavity. Then the numerical algorithm is applied to the numerical simulation of flow field around a pair of flapping Drosophila wings in hover flight. Concluding remarks are provided in Section 4. 2. Mathematical and numerical formulation The integral form of the incompressible Navier–Stokes equations that govern the motion of an arbitrary moving control volume Ω(t ) with boundary ∂Ω(t ) can be written in the Cartesian coordinate system in dimensionless form as follows: the continuity equation



−

n · u dS = 0

(1)

∂Ωe

the momentum equations



Re Ωd

∂u dV + Re ∂t



∂Ωd







n · (u − x˙ ) u dS + ∂Ωd

 np dS =

n · ∇ u dS

(2)

∂Ωd

In these equations, V is the control volume, S is the control volume surface area, n represents the outward normal vector, u represents the local fluid velocity vector, x˙ represents the grid velocity, p is the pressure and Re is the dimensionless Reynolds number. The local fluid velocity vector components are defined at the mid-point of each face meanwhile the pressure is defined at the element centroids. The above continuity equation is discretized within each hexahedral element Ωe , meanwhile the momentum equations are discretized over the dual finite volume Ωd . Fig. 1(a) illustrates the threedimensional hexahedral elements with a dual control volume.

Fig. 1. Three-dimensional unstructured mesh with a dual control volume (Ωd ) used for the velocity components (yellow volume) (a) and a covolume (Ωc ) used for the evaluation of gradient terms for the area vector A125 (red volume) (b). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

663

2.1. Three-dimensional numerical discretization The momentum equations along the x-, y- and z-directions are discretized over the dual finite volume shown in Fig. 1(a). The discrete contribution from the right cell shown in Fig. 1 is given below for each term of the momentum equation along the x-direction. The time derivative





n +1 3Re 3u 1

2

4t

i

+

uni +1



6 · 4t

n +1 V 12345





−

4Re 3un1 2

4t

+

i

uni

6 · 4t







n V 12345

+

n −1 Re 3u 1

2

4t

+

i

uni −1



6 · 4t

n −1 V 12345

i = 1, 2, 3, 4, 5, 10

with

(3)

The convective term due to fluid velocity

















+1 +1 n +1 +1 +1 n +1 +1 +1 n +1 +1 +1 n +1 Re un125 · An125 u 125 + Re un235 · An235 u 235 + Re un345 · An345 u 345 + Re un415 · An415 u 415

(4)

The convective term due to grid velocity

 +1 n +1  n +1  +1 n +1  n +1  +1 n +1  n +1  +1 n +1  n +1 −Re x˙ n125 · A125 u 125 − Re x˙ n235 · A235 u 235 − Re x˙ n345 · A345 u 345 − Re x˙ n415 · A415 u 415

(5)

The pressure term



p1 + p2 + p5

 +

n+1

3 p3 + p4 + p5 3

+1 An125

n+1

 ·i+

+1 An345

p2 + p3 + p5

 ·i+

n+1

3 p4 + p1 + p5 3

+1 An235 ·i

n+1

+1 An415 ·i

(6)

The viscous term

  n+1 n+1 n+1 ∂u ∂u ∂u +1 +1 +1 − An125 ·i+ An125 ·j+ An125 ·k ∂ x 125 ∂ y 125 ∂ z 125  n+1 n+1  n+1 ∂u ∂u ∂u n +1 n +1 n +1 A ·i+ A ·j+ A ·k − ∂ x 235 235 ∂ y 235 235 ∂ z 235 235  n+1 n+1  n+1 ∂u ∂u ∂u +1 +1 +1 An345 ·i+ An345 ·j+ An345 ·k − ∂ x 345 ∂ y 345 ∂ z 345  n+1 n+1  n+1 ∂u ∂u ∂u +1 +1 +1 An415 ·i+ An415 ·j+ An415 ·k − ∂ x 415 ∂ y 415 ∂ z 415

(7)

where V 12345 is the volume of the pyramid between the points x1 , x2 , x3 , x4 and x5 shown in Fig. 1, A125 , A235 , A345 and A415 are the area vectors of the dual volume triangular surfaces. The area vectors are computed from the cross product of vectors (as an example, A125 = 0.5(x2 − x5 ) × (x1 − x5 )). The values u125 , u235 , u345 and u415 are the velocity vectors defined at the mid-point of each dual volume triangular surfaces and p 1 , p 2 , p 3 , p 4 and p 5 are the pressure values at the points x1 , x2 , x3 , x4 and x5 , respectively. However, the pressure values are known only at the element centroids and the pressure values at x1 , x2 , x3 and x4 have to be computed. To compute the pressure at x1 , as an example, a second-order Taylor series expansion can be written as

   ∂ p  ∂ p  ∂ p  pi = p1 + (xc,i − x1 ) + ( y c ,i − y 1 ) + ( zc,i − z1 ) with i = 1, 2, . . . , m ∂ x x=x1 ∂ y x=x1 ∂ z x=x1

(8)

where m represents the number of the neighboring hexahedral elements connected to the point x1 and xc ,i corresponds to the neighboring element centroids. This overdetermined system of linear equations may be solved to compute the pressure value and the pressure gradient components in a least square sense using the normal equation approach, in which both sides are multiplied by the transpose. The modified system is solved using the singular value decomposition provided by LAPACK driver routines in order to avoid the numerical difficulties associated with solving linear systems with near rank deficiency. In a similar manner, the u-velocity component values at x1 , x2 , x3 and x4 are computed using the same approach. To compute the u-velocity component at x1 ,

ui = u1 +

   ∂ u  ∂ u  ∂ u  ( x − x ) + ( y − y ) + ( z f ,i − z1 ) with i = 1, 2, . . . , l 1 1 f ,i f ,i ∂ x x=x1 ∂ y x=x1 ∂ z x=x1

(9)

where l represents the number of the faces connected to the point x1 and x f ,i corresponds to the face mid-points. The overdetermined system of linear equations is also solved in a least square sense as before and the computed u-velocity

664

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

components are used to compute the velocity gradients defined at the mid-point of each dual volume triangular faces using the Green–Gauss theorem:

∇u =

∂u ∂u ∂u 1 i+ j+ k= ∂x ∂y ∂z Vc



u dA

(10)

∂Ωc

where V c covolume consists of two tetrahedral elements that share the same dual volume triangular surface area with their fourth vertices located at the mid-point of the hexahedral element faces. As an example, the covolume Ωc is shown in Fig. 1(b) for the triangular dual volume area vector A125 . The right-hand side of Eq. (10) is evaluated using the mid-point rule on each of the covolume faces. The convective velocity vector components u 125 , u 235 , u 345 and u 415 are computed at the mid-point of the dual volume triangular surfaces using the least square upwind interpolations [1,5]. As an example,









u 125 = β u 1 + ∇ u 1 · (x125 − x f ,1 ) + (1 − β) u 2 + ∇ u 2 · (x125 − x f ,2 )

(11)

where β is a weight factor determining the type of convection scheme used, ∇ u 1 and ∇ u 2 are the gradients of velocity components where the u 1 and u 2 velocity components are defined and x125 = (x1 + x2 + x5 )/3. For evaluating the gradient terms, ∇ u 1 and ∇ u 2 , a least square procedure is used in which the velocity data is assumed to behave linearly. Referring to Fig. 1 as an example, the following system can be constructed for the term ∇ u 1

⎡

x f ,2 − x f ,1 ⎢ x f ,3 − x f ,1 ⎢ ⎢ x f ,4 − x f ,1 ⎢ ⎢ x f ,5 − x f ,1 ⎢ ⎢ x f ,6 − x f ,1 ⎢ ⎢ x f ,7 − x f ,1 ⎣x −x f ,8 f ,1 x f ,9 − x f ,1

y f ,2 − y f ,3 − y f ,4 − y f ,5 − y f ,6 − y f ,7 − y f ,8 − y f ,9 −

y f ,1 y f ,1 y f ,1 y f ,1 y f ,1 y f ,1 y f ,1 y f ,1

⎤

z f ,2 − z f ,1 z f ,3 − z f ,1 ⎥ ⎡ ⎥ z f ,4 − z f ,1 ⎥ ⎥⎢ z f ,5 − z f ,1 ⎥ ⎢ ⎥⎢ z f ,6 − z f ,1 ⎥ ⎣ ⎥ z f ,7 − z f ,1 ⎥ z f ,8 − z f ,1 ⎦ z f ,9 − z f ,1

∂u ∂x ∂u ∂y ∂u ∂z

⎡u −u ⎤ 2 1 ⎤ ⎢ u3 − u1 ⎥ ⎢ u4 − u1 ⎥ ⎥ ⎥ ⎢ u5 − u1 ⎥ ⎥ ⎢ ⎥ ⎥= ⎢ ⎥ ⎦ ⎢ ⎢ u6 − u1 ⎥ ⎢ u7 − u1 ⎥ ⎦ ⎣

(12)

u8 − u1 u9 − u1

This overdetermined system of linear equations may be solved for ∇ u 1 in a least square sense using the same normal equation approach. The gradient term ∇ u 2 is also computed in a similar manner. To compute the fluxes due to mesh motion, a special attention is given to satisfy the discrete geometric conservation law (DGCL). The DGCL states that the volumetric increment of a moving element must be equal to the summation of the volumes swept by its surfaces that close the volume. Therefore, the mesh motion flux is evaluated as follows [19,37] +1 +1 x˙ n125 · An125 =

n −1 n  [An+1 + An125 ] 3  n +1 1  n −1  [A125 + A125 ] · x125 − xn125 · 125 − x125 − xn125 2t 2 2t 2

(13)

This approach will ensure that the DGCL is satisfied and the present ALE scheme preserves a uniform flow solution exactly independent of the mesh motion. However, Geuzaine et al. [19] showed that the compliance with the DGCL is neither a necessary nor a sufficient condition to preserve its order of time-accuracy established on fixed meshes. Because, the authors indicated by means of truncation error arguments that the linearization of the convective terms in Eq. (4) using the values at time level n will drop the accuracy of the numerical scheme to first-order on moving meshes. Hence, several sub-iterations have to be performed in order to maintain the second-order time accuracy. The discretization of the momentum equation along the y- and z-directions follows very closely the ideas presented here. It should be noted that the present dual volume surface integrals involve only triangular planar surfaces for the momentum equations which significantly simplify the three-dimensional numerical discretization. The continuity equation (1) is integrated within each hexahedral element Ωe and evaluated using the mid-point rule on each of the element faces

−

6  

u n +1 A x

i =1

 i

    + v n +1 A y i + w n +1 A z i = 0

(14)

where A = A x i + A y j + A z k 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 [7] of the form:

⎡

B 11 ⎢ 0 ⎢ ⎣ 0 B 41

0 B 22 0 B 42

0 0 B 33 B 43

⎤⎡

B 14 ⎢ B 24 ⎥ ⎥⎢ B 34 ⎦ ⎣ 0

⎤

⎡

⎤

u b1 ⎢ ⎥ v ⎥ ⎥ = ⎢ b2 ⎥ w ⎦ ⎣ b3 ⎦ p 0

(15)

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

665

where, B 11 , B 22 and B 33 are the convection diffusion operators, ( B 14 , B 24 , B 34 ) is the pressure gradient operator and ( B 41 , B 42 , B 43 ) is the divergence operator. It should also be noted that on a uniform Cartesian mesh, the multiplication of the matrices B 41 B 14 + B 42 B 24 + B 43 B 34 gives the classical five-point Laplace operator as in the MAC scheme [22] which is very important for the efficient implementation of present iterative solvers. 2.2. Iterative solver In practice, the solution of Eq. (15) does not converge very quickly and it is rather difficult to construct robust preconditioners for the whole coupled system because of the zero-block diagonal resulting from the divergence-free constraint. In the present paper, we use the following upper triangular right preconditioner which results in a scaled discrete Laplacian instead of a zero block in the original system. Then the modified system becomes

⎡

B 11 ⎢ 0 ⎢ ⎣ 0 B 41 where

⎡

⎤

0 B 22 0 B 42

0 0 B 33 B 43

⎡

⎤⎡

B 14 I 0 0 ⎢ B 24 ⎥ ⎥⎢ 0 I 0 B 34 ⎦ ⎣ 0 0 I 0 0 0 0

u I 0 0 ⎢ v ⎥ ⎢0 I 0 ⎢ ⎥=⎢ ⎣w ⎦ ⎣0 0 I p 0 0 0

⎤⎡

⎤⎡

⎤

⎡

⎤

B 14 q b1 ⎢ ⎥ ⎢ ⎥ B 24 ⎥ ⎥ ⎢ r ⎥ = ⎢ b2 ⎥ B 34 ⎦ ⎣ s ⎦ ⎣ b3 ⎦ p I 0

(16)

⎤

B 14 q ⎢r⎥ B 24 ⎥ ⎥⎢ ⎥ B 34 ⎦ ⎣ s ⎦ p I

(17)

and the zero block is replaced with B 41 B 14 + B 42 B 24 + B 43 B 34 , which is a scaled discrete Laplacian. Unfortunately, this leads to a significant increase in the number of non-zero elements due to the matrix–matrix multiplication. However, it is possible to replace the block matrices in the upper triangular right preconditioner with computationally less expensive matrices. The calculations indicate that the largest contribution for the pressure gradient in the momentum equations comes from the right and left elements that share the common face where the components of the velocity vector are discretized. Therefore, we will use the contribution from these two elements which leads to maximum three non-zero entries per row. Although, this approximation does not change the convergence rate of an iterative solver significantly, it leads to a significant reduction in the computing time and memory requirement. As an example, the two-dimensional Stokes flow in a lid-driven square cavity is solved on a uniform 201 × 201 mesh using an incomplete LU(ILU(4)) preconditioner. The original approach requires 184 iterations in order to reduce the relative residual norm to 10−8 , meanwhile the modified approach requires 192 iterations. The present one-level iterative solver is based on the restricted additive Schwarz method with the flexible GMRES(m) [44] algorithm. Because the zero block is removed, a block-incomplete factorization coupled with the reverse Cuthill–McKee ordering [14] can be used within each partitioned sub-domains. Multigrid methods [21,56] are known to be the most efficient numerical techniques for solving large-scale problems that arise in numerical simulations of physical phenomena because of their computational costs and memory requirements that scale linearly with the degrees of freedom. The two level iterative solver based on the multiplicative non-nested multigrid method with one V-cycle is described in [45] in detail and it has been successfully applied to the large-scale solution of the Stokes problems in three dimension. Although the fully coupled multigrid technique is shown to be very efficient [46], it is not suitable for the timedependent calculation of the incompressible Navier–Stokes equations with small time steps since the saddle point problem approaches to the following form:

⎡

I ⎢ 0 ⎢ ⎣ 0 B 41

0 I 0 B 42

0 0 I B 43

⎤

B 14 B 24 ⎥ ⎥ B 34 ⎦ 0

(18)

Therefore, a matrix factorization can be introduced similar to that of the projection method [13]

⎡

I ⎢ 0 ⎢ ⎣ 0 B 41

0 I 0 B 42

0 0 I B 43

⎤⎡

⎤⎡

0 I 0 0 0 I ⎢ ⎥⎢ 0 0⎥ 0 ⎥⎢ 0 I 0 ⎥⎢ ⎦⎣ 0 0 ⎦⎣ 0 0 I 0 0 0 0 − B 41 B 14 − B 42 B 24 − B 43 B 34 I 0

0 0 I 0 0 I 0 0

⎤

B 14 B 24 ⎥ ⎥ B 34 ⎦ I

(19)

The parallel algebraic multigrid solver BoomerAMG from the HYPRE library [18] is used for the inverse of the scaled discrete Laplacian. However, the Parallel Modified Independent Set (PMIS) algorithm within the BoomerAMG solver is used for coarsening scheme rather than the default Falgout algorithm due to its higher computational efficiency [50]. The implementation of the preconditioned Krylov subspace algorithm, matrix–matrix multiplication and the multilevel preconditioner are carried out using the PETSc [4] software package. The METIS library [29] is used to partition an unstructured mesh for a balanced

666

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

domain decomposition. The most appealing feature of the present fully coupled approach is that it allows us to avoid time step restriction due to the Courant–Friedrichs–Lewy condition. This is very important for large mesh deformations where the mesh deformation algorithms may lead to extremely small elements. In addition, the equation of motion of a deforming body has to be coupled with the incompressible Navier–Stokes equations for free-swimming/flying bodies [47]. 2.3. Mesh deformation algorithm There has been extensive research on mesh deformation and many different mesh deformation algorithms have been proposed in the literature in order to compute the displacement of the internal points as the boundaries of a computational domain translate, rotate and deform in order to maintain mesh quality and validity. These algorithms are based on the spring analogy [6], the elastic medium analogy [27], the edge swapping algorithm [16], the radial basis function (RBF) interpolation algorithm [9] and the remeshing algorithm [28]. The main advantage of the RBF interpolation algorithm is that it leads to high quality meshes when the domain boundaries exhibit large translations and rotations [9]. In the RBF method, an interpolation function is used to transfer the displacements known at the boundaries to the interior points.

x i =

N 





γ j φ xi − x j  + M (xi )

(20)

j =1

where N is the number of control points, γ j is the weight of control point x j , φ(xi − x j ) is the radial basis function, M (xi ) = β1 + β2 xi + β3 y i + β4 zi is a low degree polynomial and  is the Euclidean norm. There are four additional constraints due to translational and rotational rigid-body motion:

   

γj = 0

(21)

γjxj = 0

(22)

γj y j = 0

(23)

γjzj = 0

(24)

Therefore, the value of the weight coefficients for x can be evaluated by applying Eq. (20) for i = 1, 2, . . . , N using the control points:

⎡

φ11 ⎢ . ⎢ .. ⎢ ⎢ φ N1 ⎢ ⎢ 1 ⎢ ⎢ x1 ⎢ ⎣ y1 z1

· · · φ1N .. .. . . · · · φN N ··· 1 · · · xN · · · yN · · · zN

1

x1

y1

.. .

.. .

.. .

1 0 0 0 0

xN 0 0 0 0

yN 0 0 0 0

z1

⎤⎡

γ1

⎢ . .. ⎥ ⎢ . . ⎥ ⎥⎢ . ⎥ zN ⎥ ⎢ ⎢ γN ⎢ 0 ⎥ ⎥ ⎢ β1 ⎢ 0 ⎥ ⎥ ⎢ β2 0 ⎦ ⎣ β3 0 β4

⎤

⎡

x1 ⎥ ⎢ . ⎥ ⎢ .. ⎥ ⎢ ⎥ ⎢ x N ⎥ ⎢ ⎥=⎢ 0 ⎥ ⎢ ⎥ ⎢ 0 ⎥ ⎢ ⎦ ⎣ 0

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(25)

0

where φi , j = φ(xi − x j ) is the radial basis function. Once the weights γ j are computed, the deformation of internal points can be calculated from Eq. (20). However, the RBF approach is computationally expensive for large three-dimensional problems because it leads to large dense matrix systems. Rendall and Allen [41] addressed this issue by using fewer RBF control points based on an adaptive greedy algorithm. The authors indicated significant speeding up the RBF mesh deformation algorithm. In the current paper, the RBF control points are created using a relatively coarse mesh. However, the control points are not located on the fluid domain boundary as found in previous studies, but located close to the domain boundary. Therefore, the mesh points between the RBF control points and the domain boundary are moved using the analytical rigid body motion equations. This simple modification will ensure that the total volume of the fluid domain is conserved at machine precision which is important for the convergence of the employed div-stable discretization of the incompressible Navier–Stokes equations with all Dirichlet boundary conditions. In practice the implementation of the RBF mesh deformation algorithm is based on either the absolute or the relative implementations [10]. The absolute method constructs the linear system on an initial mesh and leads to very efficient implementation since the linear system needs to be factorized only once. However, the mesh quality is limited for large deformations. In the relative approach, the linear system is constructed from the node points at the previous time step. This technique was successfully used by Bos [10] for the flapping wing flight of insects. The problem with the second approach is that the mesh periodicity and quality may be lost for large time integrals. In the present paper, the indirect implementation of the RBF mesh deformation algorithm is proposed. In this approach, the rotation motion of the wing is split into several smaller substeps and the grid coordinates are always computed starting from the initial mesh. This method is also capable of handling large mesh deformations while maintaining mesh periodicity and quality. The implementation of the current algorithm is given in detail in Table 1.

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

667

Table 1 The indirect radial basis function algorithm. In here, N is the number of the control points, NP is the total number of vertices and NStep is the number of the indirect RBF iterations.

3. Numerical validations In this section, the proposed ALE algorithm is initially validated for the decaying Taylor–Green vortex flow, the flow past an oscillating circular cylinder in a channel, and the flow induced by an oscillating sphere in a cubic cavity. Then it is applied to the numerical simulation of flow field around a pair of flapping Drosophila wings.

668

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

3.1. The decaying Taylor–Green vortex flow The Taylor–Green vortex [52] is an analytical solution of the two-dimensional incompressible unsteady Navier–Stokes equations and it has been extensively used for testing and validation of the spatial and temporal order of convergence [30, 49,12]. The spatial domain is taken here as the unit square [0, 1] × [0, 1] and the analytical solution in non-dimensional form is given by 2 u (x, y , t ) = e −2π t /Re sin(π x) cos(π y )

v (x, y , t ) = −e

−2π 2 t /Re

cos(π x) sin(π y )   1 p (x, y , t ) = e −4π t /Re cos(2π x) + sin(2π y ) 4 2

(26) (27) (28)

For the present validation case the Reynolds number is taken to be Re = 10. The Reynolds number is based on the maximum velocity at the initial time t = 0, the size of the vortex and the kinematic viscosity of the fluid. The mesh deformation for the present calculations is taken to be

x(x, y , t ) = 0.1 sin(π t /0.8) sin(π x) sin(π y )

(29)

 y (x, y , t ) = 0.1 sin(π t /0.8) sin(π x) sin(π y )

(30)

The Dirichlet boundary conditions are applied for the velocity components at all time levels. The calculations are started at t = 0 by the exact solutions and are marched forward in time numerically until t = 0.4, at which time the numerical error is computed. The numerical error is taken to be:

Error =

u − u exact 2 √ Nu

(31)

where N u is the number of edges. The mesh space h on the unstructured quadrilateral element is defined as

1

h = √

Ne

(32)

where N e is the number of elements. In order to establish the spatial convergence of the method, an h-refinement study is performed on both uniform Cartesian meshes as well as unstructured quadrilateral meshes. For uniform Cartesian meshes, five different meshes are employed: mesh U1 with 21 × 21 node points, mesh U2 with 41 × 41 node points, mesh U3 with 81 × 81 node points, mesh U4 with 161 × 161 node points and mesh U5 with 321 × 321 node points. For unstructured quadrilateral meshes the following meshes are considered: mesh M1 with 432 node points and 391 elements, mesh M2 with 1729 node points and 1648 elements, mesh M3 with 6844 node points and 6683 elements, mesh M4 with 26 848 node points and 26 527 elements and mesh M5 with 105 453 node points and 104 812 elements. The successive meshes are generated using the mapping and paving algorithms provided within the CUBIT mesh generation environment [8]. In order to produce unstructured meshes with approximately uniform mesh size by the paving algorithm in a unit square, the computational domain is split into two by a circle of radius 0.35. The non-linear convective term in Eq. (4) is evaluated at time level n + 1 using two sub-iterations in order to guarantee the second-order convergence properties on moving meshes [19]. The convergence of numerical error with mesh spacing is shown in Fig. 2(a) with t = 10−4 and the numerical error decays at an algebraic rate as the mesh is refined. In here, it should be noted that the error measure is a function of mesh space h and time step t. For a sufficiently small t the numerical error is dominated by the spatial error and in a log–log scale the expected rate of convergence would appear as a straight line. The present ALE algorithm indicates an algebraic convergence rate of O (h2 ) for both structured and unstructured meshes. The convergence of numerical error with time step is also shown in Fig. 2(b) on the meshes M3, U3, M4 and U4. For a sufficiently small h the error is dominated by the temporal error for large time steps and the numerical error curves show that the ALE method has an algebraic convergence rate of O (t 2 ). However, the further reduction in time step t causes the convergence stall since the numerical error starts to be dominated by the spatial error. Then the only way to reduce the numerical error is to reduce the mesh size h. The convergence stall of the numerical error can also be seen in Fig. 2(a) if we continue to refine the mesh further with the same time step. 3.2. The flow past an oscillating circular cylinder in a channel The flow around a circular cylinder undergoing sinusoidal transverse oscillation in a channel with a specified amplitude and frequency is solved as a benchmark problem by Wan and Turek [55]. The computational domain has a size of [2.2 × 0.41]. The initial location of the cylinder center (x0 , y 0 ) is (1.1, 0.2) relative to the left bottom corner of the domain. The cylinder diameter is 0.1 and the cylinder center is oscillating sinusoidally such that the location of the cylinder center is given by x = x0 + A sin(2π f t ), where t is the time, A = 0.25 and f = 0.25 are amplitude and frequency of the oscillation,

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

√

669

√

Fig. 2. The spatial convergence of the numerical error (u − u exact 2 / N u ) with mesh refinement (h = 1/ N e ) (a) and temporal convergence of the numerical error with time step (b) for the decaying Taylor–Green flow at Re = 10.

Fig. 3. The instantaneous u-velocity component contours with streamtraces for an oscillating circular cylinder in a channel at several different time levels.

670

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

Fig. 4. The comparison of the drag and lift coefficients with the numerical results of Wan and Turek [55] for a circular cylinder undergoing sinusoidal transverse oscillations in a channel.

respectively. The kinematic viscosity ν is 1 × 10−3 m2 /s and the density is 1 kg/m3 . The calculations are started from the rest and the Dirichlet velocity boundary condition is imposed on all solid surfaces using their analytical values. The computational mesh consists of 191 014 vertices and 189 874 quadrilateral elements leading to 951 650 DOF. The mesh is highly clustered near the cylinder surface and on the lateral walls. Although there is a significant transverse translation for the cylinder, the mesh deformation algorithm is capable of handling such a large mesh deformation. The time step is set to 0.005 and the non-linear convective term in Eq. (4) is evaluated at time level n + 1 using two sub-iterations. The computed u-velocity component contours with the streamtraces are given in Fig. 3 at time levels t = 20 s, t = 21 s, t = 22 s and t = 23 s. The contour plot at t = 21 s is very similar to Fig. 5 of reference [55]. These pictures show relatively large velocities when the cylinder is passing through the channel center line. In addition, several very large separation bubbles are observed on the channel upper and lower walls. For a more accurate comparison, the computed drag coefficient 2 2 D and the lift coefficient C l = 2F y /ρ U max D are compared with the results of Wan and Turek [55] in Fig. 4. C d = 2F x /ρ U max From the comparison, we observe that the numerical results are indistinguishable from one another. 3.3. The flow due to a horizontally oscillating sphere in a cubic cavity The flow around a rigid sphere undergoing sinusoidal transverse oscillation in a cubic cavity with a specified amplitude and frequency is solved by Gilmanov and Sotiropoulos [20]. For this benchmark problem, a rigid sphere of diameter D is placed in the center of a cubic cavity with the edge length of L = 2D filled with viscous incompressible Newtonian fluid, which is initially at rest. Flow is induced by oscillating the sphere back and forth along the horizontal x-direction. The motion is initiated impulsively at t = 0 and the location of the sphere is prescribed as follows: x(t ) = h[1 − cos(2π t )] where h = 0.125D. The Reynolds number is based on the sphere diameter and the maximum sphere velocity and is taken to be Re = 20, which is the value used in the work of Gilmanov and Sotiropoulos [20]. The computational mesh consists of 492 596 vertices and 473 856 hexahedral elements leading to 4 794 264 DOF. The time step t is taken to be 0.0025 and two subiterations are used as before. The sequences of four snapshots are shown for the u-velocity contours with streamtraces in Fig. 5 at times t = 0, t = T /4, t = T /2 and t = 3T /4 where T is the period ( T = 1). At time t = T /4 and t = 3T /4 the sphere is stationary and the streamtraces are perfectly parallel to the solid surfaces. The computed streamtraces are relatively in good agreement with Fig. 6(b) in the study of Gilmanov and Sotiropoulos [20]. The time variation of the drag coefficient 2 π D 2 ) is given in Fig. 6. The drag coefficient is oscillating like sine function. The other components of C d = 4F x /(1/2ρ U max

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

671

Fig. 5. The instantaneous u-velocity component contours with streamtraces at several different time levels over a cycle for a rigid sphere oscillating in a cubic cavity at Re = 20: t = 0 (a), t = T /4 (b), t = T /2 (c) and t = 3T /4 (d).

Fig. 6. The time variation of drag coefficient for a rigid sphere oscillating in a cubic cavity at Re = 20.

the force coefficient are zero due to the flow symmetry. However, we cannot compare the maximum and minimum values of the drag coefficients since it is not available in the work of Gilmanov and Sotiropoulos [20]. Since the sphere is oscillating only in the right half part of y–z plane, the peak of drag coefficient is not symmetrical in one complete cycle.

672

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

3.4. The numerical simulation of flow field around the Drosophila wings The present ALE algorithm is applied to solve the flow field around a pair of flapping Drosophila wings in hover flight in order to further establish the reliability and accuracy of the proposed method. In the literature, the flow field around the Drosophila wings has been investigated extensively both experimentally [17,43,11,33] and numerically [51,40,26,3,2,34,10, 32,15]. The study of the insect flight reveals the sophisticated mechanism of generation of aerodynamic forces, well beyond values predicted by conventional wing, which can be considered as a combination of clap and fling, attached leading edge vortex, wake re-capture and rotational lift. We refer to the review paper by Sane [48] on the general aerodynamics of the insect flight. The use of the ALE formulations for the present problem is rather challenging due to extremely large mesh deformations which may lead to inadmissible elements. The three-dimensional wing shape and kinematics used in the present study are the same as those of robotic fly by Dickinson et al. [17]. The maximum wing chord length c is taken to be unity. The distance from the wing tip to the hinge location is 2.5 and the wing cross section has a constant thickness of 0.032. The hinge locations are 1.1 apart. The wing cross section area is computed to be 1.5578. The kinematics of the wing motion corresponds to the symmetrical rotation with respect to stroke reversal in the experiments [17] and it is approximated as follows:

α (t ) = φ(t ) =

αmax tanh C (t )

φmax sin

−1





tanh C (t ) sin(2π f t )

(0.97)





(33)



sin−1 0.97 sin 2π f (t − 0.25/ f )

(34)

where C (t ) = 1.6 + 1.6[sin(2π f t )]2 , f is the wing beat frequency, α (t ) is the wing angle of attack from the inertial frame x–z plane to the wing cross-section and φ(t ) is the wing stroke (azimuthal) angle from the inertial frame x– y plane to the wing axis of rotation. The maximum angle of attack αmax and stroke angle φmax are 40◦ and 80◦ , respectively. The maximum translational velocity at the wing tip is computed to be U max = 16.09 and the mean chord length is found to be c¯ = 0.7989. The non-dimensional Reynolds number is based on Re = U max c¯ /ν and taken to be 136 as in the experiments. The calculations are performed on a global domain of size [−5, 5] × [3, −17] × [−5, 5] as in the experiments. The computational mesh is shown in Fig. 7(a) and consists of 1 300 358 vertices and 1 276 666 hexahedral elements leading to 12 837 448 DOF. The mesh is stretched near the wing surface, in particular close to the leading and trailing edges, and it is symmetric according to the y–z plane. The mesh is created using the mapping, paving and sweeping algorithms available within the CUBIT mesh generation environment [8]. The mesh deformation is achieved by employing the indirect RBF interpolation as described in Section 2.3. The control points used for the indirect RBF method are shown in Fig. 7(b) and these points do not exactly lie on the computational domain boundary but next to the boundary. Because the use of the classical RBF method with coarsened grid points on the boundary does not guarantee that all the points on the boundary obey the rigid body motion at machine precision. However, the present simple modification ensures that all the points between the control points and the domain boundary obey the rigid body motion and the total volume of the computation domain is conserved. Therefore, the use of div-stable discretizations of the incompressible Navier–Stokes with all Dirichlet boundary conditions does not cause any convergence problem due to the incompatible boundary conditions. Although we have tested several radial basis functions such as the linear, cubic and thin plate spline, etc., the cubic radial basis function has produced the highest mesh quality. The effect of the number of steps used in the indirect radial basis function approach to mesh quality is shown in Fig. 8 with several different steps at the y = 0 plane. Apparently, the direct approach is not enough and leads to elements with negative volumes. Although there is an increase in the mesh quality with large number of steps, this increase is not significant after several steps and we employ only three steps. As boundary conditions, all Dirichlet boundary conditions are imposed on the solid walls. The magnitude of the velocity vector components on the wing surface is computed using the analytical values at the face mid-points. Meanwhile no-slip boundary conditions are used on the lateral walls of the rectangular prism. The present calculations are started impulsively with t = 0.005 and the flow becomes periodic after several cycles. The non-linear convective term in Eq. (4) is treated using two sub-iterations as before. The time variation of the Eulerian coherent structures in the near wake is analyzed with the λ2 -criterion developed by Jeong and Hussain [25] and the topological features of the flow structure is shown in Fig. 9 for the downstroke and in Fig. 10 for the upstroke. In these calculations we assume that the Drosophila body is aligned along the positive z-direction and the wing motion starts with the wing downstroke. The initial start motion of the Drosophila wings leads not only the leading and trailing edge vortices as in two dimension but also causes the formation of tip and root vortices. These vortices initially produce a U-type vortex ring around the wing edges. As the wing moves forward the tip and root vortices are continuously shed from the wing surface due to the axial flow within the leading edge vortex. Meanwhile, the leading edge vortex with its main component in the spanwise direction stays attached close to the upper leading edge and creates a low pressure region. Eventually the rotational wing motion leads to the finite C-type tip and root vortices which create the required downwash needed for lift. Due to local induced velocity, the tip and root vortices also move downward as they develop behind the wing. In addition, we observe relatively weak upward motion near the lateral solid walls since the computational domain is constrained with all Dirichlet (no-slip) boundary conditions. These observed vortices do not form a closed doughnut-shaped vortex ring as seen

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

673

Fig. 7. The computational mesh with 1 300 358 vertices and 1 276 666 hexahedral elements (12 837 448 DOF) (a) and the control points (1308 nodes) used for the RBF based mesh deformation algorithm (b) for the Drosophila wings.

in the numerical simulation of Aono et al. [2] where the authors presented the near- and far-field wake structure of a fully integrated wing–body combination using a different kinematic model. The computed wake structure is in very good agreement with the wake structures shown in Fig. 4 of Kweon and Choi [32] even though the authors’ simulation involves only a single wing. In addition, we observe topological similarities in the wake structure with the experimental flow visualization of Poelma et al. [39] around an impulsively-started wing at a constant angle of attack. The magnitude of the total forces and moments acting on the Drosophila left wing hinge location and the total power requirement for flapping are computed from the following surface integrals as a function of time



F(t ) =

(n · σ ) dS

(35)

∂Ω w



M(t ) =

r × (n · σ ) dS

(36)

n · (σ · u) dS

(37)

∂Ω w



P (t ) = ∂Ω w

where σ is the stress tensor including the pressure term, n · σ is the traction vector, r is the distance to the hinge location and ∂Ω w represents the Drosophila left wing surface. In the case of free-swimming/flying bodies, Eqs. (35) and (36) are needed to be written for the equation of motions and have to be coupled with the discrete incompressible Navier–Stokes 2 2 equations given in Eq. (15). The computed force coefficients (F/0.5ρ U max S ), moment coefficients (M/0.5ρ U max S c¯ ) and 3 power coefficient ( P /0.5ρ U max S ) are shown in Fig. 11 and the total lift coefficient is compared with the experimental result of Dickinson et al. [17] and the numerical result of Kweon andChoi [32].Due to the symmetry of the flow structure, the  following total force and moments are zero at each time level: F x = 0, M y = 0 and M z = 0. Although the computed lift coefficient is in a relatively good agreement with the numerical result of Kweon and Choi [32], where the authors

674

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

Fig. 8. The initial mesh (a) and deformed meshes at y = 0 plane using indirect radial basis function method with 1 iteration (direct method) (b), 2 iterations (c), 3 iterations (d), 4 iterations (e) and 5 iterations (f).

approximated the wing kinematics using smooth cubic spline functions, there is some difference with the experimental results of Dickinson et al. [17]. The slight difference with the numerical study of Kweon and Choi [32] is believed to be due to the lack of wing–wing interactions in their simulations. Although several researchers produced better agreements with the experiment, the variations in the precise kinematic patterns and wing morphologies make it very difficult to compare. As an example, Ramamurti and Sandberg [40] used similar wing kinematics with the same stroke and angle of attack

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

675

Fig. 9. The instantaneous downstroke wake structures (λ2 -criterion) around a pair of Drosophila wings at several different time levels: t = 0.00T (a), t = 0.10T (b), t = 0.20T (c), t = 0.30T (d), t = 0.40T (e) and t = 0.50T (f).

amplitudes, however, the wing stroke amplitude is between 90◦ and −70◦ rather than ±80◦ . Dai et al. [15] modified the ramp during stroke reversal in order to have a better agreement with the experiments. In the present simulations, it is quite interesting that the amplitude of the forces F x and F z is comparable to the lift force F y . The large sudden increase in lift after the stroke reversal is explained as rotational lift by Dickinson et al. [17]. The translational motion of the wing causes the fluid to accelerate towards to the wing, in particular, within the tip vortices. Even though the wing is decelerated and stopped at t = 0.5T , it still experiences relative velocity due to its wake resulting in non-zero positive drag and zero lift since the wing is parallel to y-axis. In addition, the energy dumped into the wake is partially recovered. During the

676

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

Fig. 10. The instantaneous upstroke wake structures (λ2 -criterion) around a pair of Drosophila wings at several different time levels: t = 0.50T (a), t = 0.60T (b), t = 0.70T (c), t = 0.80T (d), t = 0.90T (e) and t = 1.00T (f).

deceleration process, the bound vortex detaches from the wing surface and the wing rotation motion causes the formation of a trailing edge vortex which then will interact with the wing during the preceding stroke reversal as seen in Fig. 9(b) and Fig. 10(b). The resulting force due to the wing wake interaction highly depends on the initial position of previously created vortices which is very sensitive to wing kinematics. The calculations are also indicated that the lift force on the Drosophila wing surface is mainly produced in the regions close to the wing tips and the leading edges where the stable leading edge vortices attach themselves as shown in Fig. 12. The streamline directions are towards to the lower pressure regions on the wings. These finite C-type vortices create the downwash required for lift.

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

677

Fig. 11. The computed total force coefficients (a), the total moment coefficients (b), the power coefficient (c) for the Drosophila wing and the comparison of the lift coefficient with the experimental work of Dickinson et al. [17] and the numerical work of Kweon and Choi [32] (d).

Fig. 12. The computed instantaneous pressure coefficient contours during the downstroke motion at t = 0, t = 0.1T , t = 0.2T , t = 0.3T , t = 0.4T and t = 0.5T on the upper left wing surface (a) and the instantaneous v-velocity component isosurfaces with the streamtraces showing the stable leading edge vortex (b).

678

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

4. Conclusions A novel ALE approach based on the side-centered unstructured finite volume method has been presented for large-scale simulation of moving boundary problems in a fully coupled form. The numerical formulation uses the staggered arrangement of the primitive variables in order to avoid odd-even pressure decoupling or spurious pressure modes on unstructured meshes. The most appealing feature of the present approach is that it leads to very robust numerical algorithm when it is combined with multigrid methods [46] similar to that of the classical MAC method [54]. In the current discretization, the continuity equation is satisfied within each hexahedral elements at machine precision and the summation of the discrete equations can be exactly reduced to the domain boundary, which is important for the global mass conservation. In addition, a special attention is given to construct a second-order ALE algorithm obeying the discrete geometric conservation law (DGCL). A mesh deformation algorithm based on the indirect radial basis function method has been proposed in order to handle large mesh deformations caused by translations and rotations. The resulting system of linear algebraic equations are solved in a fully coupled manner using a matrix factorization similar to that of the projection method [13] and the BoomerAMG, a parallel algebraic-multigrid preconditioner, is employed for the scaled discrete Laplacian provided by the HYPRE library [18] that is called trough the PETSc library [4] interface. The present numerical algorithm is validated for the classical benchmark problems in the literature and then it is applied to the numerical simulation of flow field around a pair of flapping Drosophila wings in hover flight. The present fully-coupled ALE algorithm is sufficiently robust to deal with large mesh deformations seen in flapping wings and reveals highly detailed near wake topology which is very useful to study physics in biological flights and can also provide an effective tool for designing bio-inspired MAVs. Acknowledgements The authors acknowledge financial support from Turkish National Scientific and Technical Research Council (TUBITAK) through project number 111M332. The authors would like to thank Michael Dickinson and his graduate student Michael Elzinga for providing the experimental data. The authors are grateful for 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 of Turkey (UYBHM) under grant number 10752009 and the computing facilities at TUBITAK-ULAKBIM, High Performance and Grid Computing Center. References [1] W.K. Anderson, D.L. Bonhaus, An implicit upwind algorithm for computing turbulent flows on unstructured grids, Comput. Fluids 23 (1994) 1–21. [2] H. Aono, F. Liang, H. Liu, Near- and far-field aerodynamics in insect hovering flight: an integrated computational study, J. Exp. Biol. 211 (2008) 239–257. [3] P. Bai, E. Cui, F. Li, W. Zhou, B. Chen, A new bionic MAV’s flapping motion based on fruit fly hovering at low Reynolds number, Acta Mech. Sin. 23 (2007) 485–493. [4] S. Balay, K. Buschelman, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, B.F. Smith, H. Zhang, PETSc users manual, ANL-95/11, Mathematic and Computer Science Division, Argonne National Laboratory, 2004, http://www-unix.mcs.anl.gov/petsc/petsc-2/. [5] T.J. Barth, A 3-D upwind Euler solver for unstructured meshes, AIAA Paper 91-1548-CP, 1991. [6] J.T. Batina, Unsteady Euler airfoil solutions using unstructured dynamic meshes, AIAA J. 28 (1990) 1381–1388. [7] M. Benzi, G.H. Golub, J. Liesen, Numerical solution of saddle point problems, Acta Numer. 14 (2005) 1–37. [8] T.D. Blacker, S. Benzley, S. Jankovich, R. Kerr, J. Kraftcheck, R. Kerr, P. Knupp, R. Leland, D. Melander, R. Meyers, S. Mitchell, J. Shepard, T. Tautges, D. White, CUBIT Mesh Generation Environment Users Manual, vol. 1, Sandia National Laboratories, Albuquerque, NM, 1999. [9] A. de Boer, M.S. van der Schoot, H. Bijl, Mesh deformation based on radial basis function interpolation, Comput. Struct. 85 (2007) 784–795. [10] F.M. Bos, Numerical simulations of flapping foil and wing aerodynamics, PhD Thesis, Delft University of Technology, Delft, 2009. [11] J. Birch, M. Dickinson, The influence of wingwake interactions on the production of aerodynamic forces in flapping flight, J. Exp. Biol. 206 (2003) 2257–2272. [12] C.S. Chew, K.S. Yeo, C. Shu, A generalized finite-difference (GFD) ALE scheme for incompressible flows around moving solid bodies on hybrid meshfree– Cartesian grids, J. Comput. Phys. 218 (2006) 510–548. [13] A.J. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comput. 22 (1968) 745–762. [14] E. Cuthill, J. McKee, Reducing the bandwidth of sparse symmetric matrices, in: Twenty Fourth ACM National Conference, 1969, pp. 157–172. [15] H. Dai, H. Luo, J.F. Doyle, Dynamic pitching of an elastic rectangular wing in hovering motion, J. Fluid Mech. 693 (2012) 473–499. [16] M. Dai, D.P. Schmidt, Adaptive tetrahedral meshing in free-surface flow, J. Comput. Phys. 208 (2005) 228–252. [17] M.H. Dickinson, F.-O. Lehmann, S. Sane, Wing rotation and the aerodynamic basis of insect flight, Science 284 (1999) 1954–1960. [18] R. Falgout, A. Baker, E. Chow, V.E. Henson, E. Hill, J. Jones, T. Kolev, B. Lee, J. Painter, C. Tong, P. Vassilevski, U.M. Yang, Users manual, HYPRE high performance preconditioners, UCRL-MA-137155 DR, Center for Applied Scientific Computing, Lawrence Livermore National Laboratory, 2002, http://www.llnl.gov/CASC/hypre/. [19] P. Geuzaine, C. Grandmont, C. Farhat, Design and analysis of ALE schemes with provable second-order time-accuracy for inviscid and viscous flow simulations, J. Comput. Phys. 191 (2003) 206–227. [20] A. Gilmanov, F. Sotiropoulos, A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies, J. Comput. Phys. 207 (2005) 45–492. [21] W. Hackbusch, Multigrid Methods and Applications, Springer-Verlag, Heidelberg, 1985. [22] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, J. Comput. Phys. 8 (1965) 2182–2189. [23] C.W. Hirt, A.A. Amsden, J.L. Cook, An arbitrary Lagrangian Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227–253. [24] Y.H. Hwang, Calculations of incompressible flow on a staggered triangle grid, Part I: Mathematical formulation, Numer. Heat Transf., Part B, Fundam. 27 (1995) 323–1995. [25] J. Jeong, F. Hussain, On the identification of a vortex, J. Fluid Mech. 285 (1995) 69–94.

B. Erzincanli, M. Sahin / Journal of Computational Physics 255 (2013) 660–679

679

[26] A. Johnson, Dynamic-mesh CFD and its applications to flapping-wing micro-aerial vehicles, in: The 25th Army Science Conference, ASC, Orlando, FL, 27–30 November, 2006. [27] A. Johnson, T. Tezduyar, Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces, Comput. Methods Appl. Mech. Eng. 119 (1994). [28] A. Johnson, T. Tezduyar, Advanced mesh generation and update methods for 3D flow simulations, Comput. Mech. 23 (1999) 130–143. [29] G. Karypis, V. Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM J. Sci. Comput. 20 (1998) 359–392. [30] D. Kim, H. Choi, A second-order time-accurate finite volume method for unsteady incompressible flow on hybrid unstructured grids, J. Comput. Phys. 162 (2000) 411–428. [31] B. Koobus, C. Farhat, Second-order time-accurate and geometricaly conservative implicit schemes for flow computations on unstructured dynamic meshes, Comput. Methods Appl. Mech. Eng. 170 (1999) 103–129. [32] J. Kweon, H. Choi, Sectional lift coefficient of a flapping wing in hovering motion, Phys. Fluids 22 (2010) 071703. [33] F.O. Lehmann, S.P. Sane, M.H. Dickinson, The aerodynamic effects of wing–wing interaction in flapping insect wings, J. Exp. Biol. 208 (2005) 3075–3092. [34] Z. Liang, H. Dong, Unsteady aerodynamics and wing kinematics effect in hovering insect flight, in: Proceedings of 47th AIAA Aerospace Sciences Meeting and Exhibit with New Horizons Forum, AIAA, 2009, paper 2009-1299. [35] C.R. Maliska, G.D. Raithby, A method for computing three-dimensional flows using non-orthogonal boundary-fitted co-ordinates, Int. J. Numer. Methods Fluids 4 (1984) 519–537. [36] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. [37] A. Naderi, M. Darbandi, M. Taeibi-Rahni, Developing a unified FVE-ALE approach to solve unsteady fluid flow with moving boundaries, Int. J. Numer. Methods Fluids 63 (2010) 40–68. [38] C.S. Peskin, The fluid dynamics of heart valves: experimental, theoretical, and computational methods, Annu. Rev. Fluid Mech. 14 (1982) 235–259. [39] C. Poelma, W.B. Dickson, M.H. Dickinson, Time-resolved reconstruction of the full velocity field around a dynamically-scaled flapping wing, Exp. Fluids 41 (2006) 213–225. [40] R. Ramamurti, W.C. Sandberg, A three-dimensional computational study of the aerodynamic mechanisms of insect flight, J. Exp. Biol. 205 (2002) 1507–1518. [41] T. Rendall, C. Allen, Efficient mesh motion using radial basis functions with data reduction algorithms, J. Comput. Phys. 228 (2009) 6231–6249. [42] S. Rida, F. McKenty, F.L. Meng, M. Reggio, A staggered control volume scheme for unstructured triangular grids, Int. J. Numer. Methods Fluids 25 (1997) 697–717. [43] S.P. Sane, M.H. Dickinson, The control of flight force by a flapping wing: Lift and drag production, J. Exp. Biol. 204 (2001) 2607–2626. [44] Y. Saad, A flexible inner–outer product preconditioned GMRES algorithm, SIAM J. Sci. Stat. Comput. 14 (1993) 461–469. [45] M. Sahin, A stable unstructured finite volume method for parallel large-scale viscoelastic fluid flow calculations, J. Non-Newton. Fluid Mech. 166 (2011) 779–791. [46] M. Sahin, Parallel large-scale numerical simulations of purely-elastic instabilities behind a confined circular cylinder in a rectangular channel, J. NonNewton. Fluid Mech. 195 (2013) 46–56. [47] M. Sahin, K. Mohseni, An arbitrary Lagrangian–Eulerian formulation for the numerical simulation of flow patterns generated by the hydromedusa Aequorea victoria, J. Comput. Phys. 228 (2009) 4588–4605. [48] S.P. Sane, The aerodynamics of insect flight, J. Exp. Biol. 206 (2003) 4191–4208. [49] R.W. Smith, J.A. Wright, An implicit edge-based ALE method for the incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 43 (2003) 253–279. [50] H.D. Sterck, U.M. Yang, J.J. Heys, Reducing complexity in parallel algebraic multigrid preconditioners, SIAM J. Matrix Anal. Appl. 27 (2006) 1019–1039. [51] M. Sun, J. Tang, Unsteady aerodynamic force generation by a model fruit fly wing in flapping motion, J. Exp. Biol. 205 (2002) 55–70. [52] G.I. Taylor, On the decay of vortices in a viscous fluid, Philos. Mag. 46 (1923) 671–675. [53] P.D. Thomas, C.K. Lombard, Geometric conservation law and its application to flow computations on moving grids, AIAA J. 17 (1979) 1030–1037. [54] S.P. Vanka, Block-implicit multigrid solutions of Navier–Stokes equations in primitive variables, J. Comput. Phys. 65 (1986) 138–158. [55] D. Wan, S. Turek, Fictitious boundary and moving mesh methods for the numerical simulation of rigid particulate flows, J. Comput. Phys. 222 (2007) 28–56. [56] P. Wesseling, An Introduction to Multigrid Methods, John Wiley and Sons, New York, 1992.