Radial basis function based mesh deformation applied to simulation of flow around flapping wings

Computers & Fluids 79 (2013) 167–177

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Radial basis function based mesh deformation applied to simulation of flow around flapping wings Frank M. Bos, Bas W. van Oudheusden, Hester Bijl ⇑ Department of Aerospace Engineering, Delft University of Technology, Kluyverweg 2, P.O. Box 5058, Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 24 November 2011 Received in revised form 24 January 2013 Accepted 4 February 2013 Available online 7 March 2013 Keywords: Mesh deformation Radial basis functions OpenFOAMÒ Flapping wings

a b s t r a c t In this paper a mesh deformation technique based on radial basis function interpolation (introduced by the authors in [6]) is applied to flapping wings. The main difficulty of the simulation of flows around flapping wings is to maintain high mesh quality when the wing exhibits large translations and rotations. Standard mesh deformation techniques – based on solving the Laplace and solid body rotation stress equations – become inefficient or break down under these conditions. Using radial basis function interpolation, mesh quality can be preserved even for large translations and rotations. First, the performance of several radial basis functions is compared to that of standard meshmotion methods in 2D for a beam undergoing rotations and translations typical for flapping. After that, the best performing radial basis function – the thin plate spline – is applied to 3D flow simulation around a flapping wing as well as to a deforming wing. Efficiency of the mesh motion technique is improved using coarsening and smoothing of the radial basis function. The resulting method is efficient and provides superior mesh quality over standard Laplace based or solid body rotation stress mesh motion methods for flapping and deforming wings. Ó 2013 Published by Elsevier Ltd.

1. Introduction Flow around flapping wings at low Reynolds numbers is a challenging topic for numerical flow simulation methods. To begin with, the wing motion consists of large displacements and rotations, which is a major challenge for mesh motion solvers. In addition, it is crucial to correctly capture the Leading Edge Vortex, so mesh quality close to the wing is of utmost importance and subsequently posing additional challenges to the flow solver. Three main approaches to solve for flow around moving objects can be distinguished: 1. The immersed boundary method (IBM). 2. Overset grids. 3. The Arbitrary Lagrangian Eulerian (ALE) method. In the immersed boundary method [25] the equations are discretised on a permanent Cartesian background mesh, through which the body moves. The presence and movement of the body is represented by boundary conditions, that represent the flow-solid interface. While for non-moving objects, definition of the appropriate boundary conditions is already non-trivial, it is even harder when the object moves. See [24] for a discussion of several ⇑ Corresponding author. Tel.: +31 (0)15 27 85373; fax: +31 (0)15 27 87077. E-mail address: [email protected] (H. Bijl). 0045-7930/$ - see front matter Ó 2013 Published by Elsevier Ltd. http://dx.doi.org/10.1016/j.compfluid.2013.02.004

boundary treatments. In [31], a different boundary treatment using implicit force calculations is proposed. Compared to other numerical studies, the resulting forces show reasonable accuracy. Several researchers have applied the IBM method to solve for flow around flapping wings [38,20]. Especially for clap and fling, where wings touch each other IBM is the method of choice [24,23]. Disadvantages of the immersed boundary method are the difficulty to accurately capture the boundary layer and to meet the requirements for mass and momentum conservation, which becomes progressively important with increasing value of the Reynolds number. Therefore this method is not considered here. The Overset grid approach uses more than one grid to represent the flow field. For flapping wing simulations, a boundary fitted grid can be attached to the moving wing, while the background can be represented by a static Cartesian grid. As the grids overlap, interpolation must be used to exchange information between the grids. Especially in simulations with multiple moving wings, the method of Overset grids can be convenient. In this case each wing can have their own mesh. Examples of simulation of flow around rigid wings and a body by means of Overset grids are presented in [33,32,1]. A disadvantage of this method is the limited conservation due to interpolation between the grids. In the Arbitrary Lagrangian Eulerian (ALE) method the equations are discretised on a deforming or moving flow mesh. In the ALE method, a mesh deformation algorithm is needed to compute the motion of the internal points of the flow mesh, based on the

168

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

movement of the body. The ALE method is applied and discussed in [22,11,10,8]. Different mesh motion algorithms exist. With this method it is possible to obtain high mesh quality inside the boundary layer, in combination with conservation of mass and momentum. Potential critical issues are the computational time needed for solution of the mesh motion equations and the potential loss of mesh quality and validity during the simulation. In order to assess the quality of a mesh motion solver, different aspects are investigated, regarding in particular mesh quality and mesh motion efficiency. The quality of the resulting mesh is defined by the non-orthogonality and skewness of the finite volume cells. Efficiency is a measure of the used computation time to calculate the displacements of the mesh points at the new time-step. A variety of mesh deformation methods have been presented in literature using different approaches to calculate the motion of the mesh points. For structured meshes, there are efficient techniques available to deform the mesh, for example Transfinite Interpolation [35]. The displacements of the boundary points are in this method interpolated along grid lines to find the displacements of all interior mesh points. Using an additional mapping [34,21] the mesh quality can be strongly improved if the boundary is subjected to significant rotation and deformation. These methods are perfectly suitable for structured but unsuitable for unstructured grids. Since unstructured meshes are used for complex geometries, possibly in combination with mesh refinement, the focus is put on mesh deformation techniques which can be applied to unstructured meshes containing arbitrary polyhedral cells, used in the finite volume code of OpenFOAMÒ [36,17]. The most popular mesh deformation method, applicable to both structured and unstructured meshes, is called the spring analogy [3] where the point-to-point connection of every two neighbouring mesh points is represented by a linear spring. However, this method proved to lack robustness, especially on arbitrarily unstructured meshes, as was observed by [5], high resulting mesh quality was only achieved by specifying a problem specific spring stiffness. Additionally, [10,7] proposed a method to incorporate torsional springs to improve the robustness of this method. Other mesh deformation techniques involve solving a partial differential equation on the complete field of internal mesh displacements for given boundary point displacements. Concerning the governing partial differential equations, the Laplace and bi-harmonic operators [22,13] are often used in combination with a constant or variable distance-based diffusion coefficient to improve the mesh quality. Another choice of equations is made by [18] who used the pseudo-solid equation, which assumes static equilibrium for small deformations of a linear elastic solid (the mesh is treated as if it was a solid). The latter method is often used in the Arbitrary Lagrangian Eulerian formulation of finite element codes. [9] modified this method to incorporate rigid body rotation, which significantly improves the mesh quality for meshes subjected to large boundary translations and rotations. Since these considered mesh motion methods are using partial differential equations on the complete field of internal mesh points, existing iterative solvers can be used [16,15]. Therefore, the parallel implementation of these methods is fairly straightforward. Depending on the method, a variable diffusivity field needs to be defined, which acts as a stiffness of the system of equations and this influences the efficiency. One major drawback of these methods is that they all fail in maintaining high mesh quality when the boundary points move with high rotation angles. In this paper the mesh deformation methods based on the use of radial basis function (RBF) interpolation – developed in our group [6] – is further improved and applied to flapping wing problems. An RBF interpolation function is used to transfer prescribed boundary point displacements to the fluid boundary mesh. In [6] it was shown that radial basis function interpolation could improve

mesh quality considerably. With the RBF method even for large displacements and rotations mesh quality close to the object is preserved. Drawback of the method is that inversion of a full matrix with dimension equal to the number of boundary points. As a result it is computationally intensive. Recently two ideas to speed up the RBF algorithm were presented [14,27,26]. Since RBF mesh deformation is capable of dealing with large rotation angles, this technique is well suited for flapping wing simulations. In this paper, the suitability of the RBF mesh motion for this type of simulations is investigated. Hereto, the performance (resulting mesh quality, efficiency and robustness) of the RBF method is compared to that of two other popular mesh motion algorithms. In addition, a way to further improve the efficiency of the RBF method is presented. To demonstrate its capability, the resulting algorithm is applied to the simulation of flow around flapping wings. In this paper, two existing mesh deformation techniques, based on the Laplace equation with variable diffusivity and a modified pseudo-solid equation, are compared with radial basis function mesh motion for flapping cases. Both Laplace and pseudo-solid mesh motion techniques are commonly used within the OpenFOAMÒ community. These different mesh motion methods are described in Section 2. In order to assess the mesh quality, Section 2.5 discusses two different criteria, based on non-orthogonality and skewness. The resulting mesh quality of the different mesh motion methods is studied using a two-dimensional test case of a block which translates and rotates. The mesh quality is investigated using grid visualisation and histograms of the skewness and nonorthogonality criterion. This is the subject of Section 3. In addition to the simplified moving block, two more realistic test cases were considered, one using a three-dimensional flapping wing and the other involves a two-dimensional flexing airfoil. Since the radial basis function mesh motion method is computationally expensive, Section 2.4 deals with two techniques to increase its efficiency. Finally, the conclusions are drawn in Section 4. 2. Computational methods When a moving mesh problem is considered, the shape of the computational domain is varying in time. Therefore, a distinction can be made between the motion of the boundary points and the motion of the internal (fluid) points. For the displacement of the boundary, points can be considered to be given, where it is either externally defined, e.g. a prescribed rigid body motion, or it is part of the solution, as in the case of fluid–structure interaction problems. According to the given boundary point motion, the internal points need to be moved such as to maintain mesh quality and validity. The internal point motion influences the solution only through the discretization errors [11], provided that the ALE formulation is correctly implemented. The internal point motion can be calculated using different methods, as will be discussed in the next section. 2.1. Laplace equation with variable diffusivity One can think of a deforming computational domain as if it was a solid body undergoing internal stresses given by the Piola–Kirchhoff stress–strain equation [2]. That equation is non-linear and thus expensive to solve using standard numerical techniques. Therefore, other type of equations have been used, namely the Laplace equation and the solid body rotation stress (SBR Stress) equation, which is a variant of the linear stress equation [9]. A mesh motion method based on one of these equations is computationally cheap since the resulting matrix system is sparse, such that existing iterative solvers can be used efficiently.

169

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

When the mesh motion is governed by the Laplace equation, the given boundary point motion may be arbitrary and non-uniform. The nature of the Laplace equation is that the point displacements will be largest close to the moving boundary and small at large distance. Ideally, a user input is not desired, since that would decrease the robustness of the method. However, this method needs the specification of a variable diffusivity. This leads to the following definition of the Laplace equation:

$  ðcrxÞ ¼ 0; where x is the displacement field and c the diffusion coefficient, which decreases with the radius r from the deforming boundary as follows:

1 cðrÞ ¼ m : r

Solving the Laplace or the SBR Stress equation leads to a sparse system of equations, such that standard iterative techniques can be used, like the pre-conditioned Conjugate Gradient (PCG) method. 2.3. Radial basis function interpolation In the current work we use radial basis function interpolation to find the displacements of the internal fluid points for given boundary displacements. The interpolation function s(x) describing the displacement of all computational mesh points, is approximated by a sum of basis functions:

sðxÞ ¼

Nb X

cj /ðkx  xbj kÞ þ qðxÞ;

ð7Þ

j¼1

ð1Þ

The resulting mesh quality strongly depends on the chosen c(r) function, which depends on the distance from the moving boundary. This variable diffusion coefficient can be chosen such that a region next to the deforming or moving boundary closely moves with the boundary. The resulting mesh contains less cell quality deterioration next to the boundary. The current research uses a c(r) function like Eq. (1). In addition to the freedom of choosing a diffusion function, it is also possible to define c(r) for every internal mesh cell for all time-steps independently. This, however, appears to be problem dependent and thus optimisation of c(r) seems not cost effective. To maintain robustness, in the current work we use a quadratically, m = 2, decreasing diffusion coefficient, which was found to provide efficient in producing a smooth mesh motion [16]. Alternatively, one could also have used an exponentially decreasing diffusion coefficient or a diffusion coefficient related to the mesh deformation energy.

where the known boundary value displacements are given by xbj ¼ ½xbj ; ybj ; zbj , q(x) is a polynomial, Nb is the number of boundary points and / is a given basis function as a function of the Euclidean distance kxk. The minimal degree of polynomial q(x) depends on the particular choice of the basis function / [6]. A unique interpolant is obtained if the basis function is a conditionally positive definite function. If the basis functions are conditionally positive definite of order m 6 2, a linear polynomial for q(x) can be used [4]. We only applied basis functions that satisfy this criterion. A consequence of using a linear polynomial is that rigid body translations are exactly recovered. The polynomial q(x) is defined by the coefficients cj which can be defined by evaluating the interpolation function s(x) in the known boundary points:

sðxbj Þ ¼ Dxbj : Here Dxbj contains the known discrete values of the boundary point displacements. Together with the additional requirements: Nb X

2.2. Solid body rotation stress equation

cj pðxbj Þ ¼ 0;

j¼1

The second method to deform the mesh is based on the linear elasticity equation and is called the solid body rotation stress (SBR Stress) equation [9]. The equation of linear elasticity, valid for small displacements, may be written as

$  r ¼ 0;

ð2Þ

where r is the stress tensor. The stress tensor r is given in terms of the strain, which is given by the following constitutive relation:

r ¼ kTrðÞI þ 2l;

ð3Þ

where Tr indicates the trace and k and l are Lamé constants [2], which are properties of the elastic material. The following equation:

1 2

 ¼ ðrx þ rxT Þ;

ð4Þ

defines the relative change in length, where x is the position of an internal mesh point, which is treated as if it was a linear solid. Although Eq. (4) represents a solid body motion, it can be extended by inclusion of rotation represented by strain. In [9], an extra term was added to obtain the following strain relation:

1 2

 ¼ ðrx þ rxT þ rxT  rxÞ:

ð5Þ

Combining Eqs. (5), (3) and (2), together with k = l, the following solid body rotation stress equation is obtained:

$  ðcrxÞ þ $  ðcð$x  $xT Þ  kTrð$xÞIÞ ¼ 0;

ð6Þ

where c is a similar diffusion coefficient as in Eq. (1). Eq. (6) allows for rigid body motion and is still linear and therefore the computational costs are of the same order as the costs necessary to solve the Laplace equation.

which holds for all polynomials p(x) with a degree less or equal than that of polynomial q(x), the cj values can be determined [6]. The values for the coefficients cj and the linear polynomial can be obtained by solving the system:



  Ubb Dx b ¼ Q Tb 0

Qb

 

0

b

c

;

ð8Þ

where c is containing all coefficients cj, b the coefficients of the linear polynomial q(x), Ubb an nb  nb matrix contains the evaluation of the basis function /bi bj ¼ /ðkxbi  xbj kÞ and can be seen as a connectivity matrix connecting all boundary points with every other boundary point. Qb is an (nb  (d + 1)) matrix with row j given by   1 xbj , where d is the spatial dimension. In general, (8) leads to a dense matrix system, which is difficult to solve using standard iterative techniques. Therefore, it needs to be solved directly, by doing a LU decomposition. The possibilities of solving the system in a more efficient way are discussed in Section 2.4. After the coefficients in c and b having been obtained, they are used to calculate the values for the displacements of all internal fluid points Dxinj using the evaluation function (7),

Dxinj ¼ sðxinj Þ:

ð9Þ

The result of (9) is transferred to the mesh motion solver to update all internal points accordingly. The value of the interpolation function is equal to the actual displacement at the moving boundary and zero at the stationary outer boundaries. Every internal mesh point is moved based on its calculated displacement, such that no mesh connectivity is necessary. The size of the system of this method (8) is ((Nb + (d + 1))  (Nb + (d + 1))), where d is the spatial dimension, which is considerably smaller than for other techniques using

170

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

the mesh connectivity. The mesh connectivity techniques encounter systems of the order (Nin  Nin), with Nin the total number of mesh points, which is a dimension higher than the total number of boundary points. Inversion of the full (Nb  Nb) matrix takes approximately N3b operations, while the inversion of the sparse matrix (Nin  Nin) depends on the bandwidth b, Nin  b2. Assuming a 3D structured mesh, for a mesh with N3 elements the bandwidth will be N2, while for a 2D mesh with N2 elements, the bandwidth will be N. For relevant problems, the number of boundary points is only a fraction of the total number of internal points, such that it will be faster to invert the full (Nb  Nb) matrix compared to the sparse (Nin  Nin) matrix. Since the number of operations for the inversion is proportional to N3b , reducing the number of boundary points has a significant influence on the speed of the RBF method. Solving the system (8) provides the values of the coefficients c and b, which are then used for step two, the evaluation using Eq. (9). In contrast to the Laplace and SBR Stress methods, no partial differential equation needs to be solved and the evaluation of all internal boundary points is straightforward to implement in parallel, since no mesh connectivity is required. Concerning robustness, this method is not using a variable diffusion coefficient which has to be tuned by the user. Instead, the proper radial basis function needs to be chosen to satisfy the need for robustness. In general, the radial basis function should be smooth and global, but to what extent for a desired accuracy is subject of the present study.

2.3.1. Radial basis functions with compact support From literature, various radial basis functions are available, which are suitable for data interpolation. Two types of radial basis functions can be distinguished: functions with compact or global support. Functions with compact support have the following property:

 /ðx=rÞ ¼

f ðx=rÞ 0 6 x 6 r; 0

x > r;

where f(x/r) P 0 is scaled with a support radius r. When a support radius is used, only the internal mesh points inside a circle (twodimensional problem) or a sphere (three-dimensional problem) with radius r around a centre xj are influenced by the movement of the boundary points. When choosing r, it must be noted that larger values for the support radius lead to more accurate mesh motion. On the other hand, a larger support radius leads to a denser matrix system, while a low support radius results in a sparse system which can be solved efficiently. Additionally, when the change in position is given by an interpolated function with compact support, where the support radius does not reach the far-field, then it is not necessary to include the far-field boundary points, providing no polynomial terms are used. Ref. [6] compared the resulting mesh quality using a selection of radial basis functions. The best results were obtained using the continuous polynomial C2, the functions based on a continuous thin plate spline, C1, C 2a and C 2b and finally the globally supported thin plate spline, see Table 1. These radial basis functions are applied to our test problem and the best performing one is used throughout the current investigation.

2.3.2. Radial basis functions with global support In contrast to the functions with compact support, functions with global support are not equal to zero outside a certain radius, but cover the whole interpolation space. Radial basis functions with global support generally lead to dense matrix systems, which can be improved by multiplication with a smoothing function, as will be discussed in Section 2.4.

Table 1 Different used radial basis functions with compact (CP, CTPS) and global (TPS) support. Radial basis functions with compact support are non-zero within the range of the support radius r. Note that n = x/r. Taken from [37]. Radial basis functions with global support cover the whole interpolation space, i.e. the computational domain. RBF name

f(n)

2

(1  n)4(4n + 1)

CP C CTPS C1 CTPS

2 3 4 2 8 5 1 þ 80 3 n  40n þ 15n  3 n þ 20n logðnÞ 1  30n2  10n3 + 45n4  6n5  60n3log (n)

C 2a C 2b

1  20n2 + 80n3  45n4  16n5 + 60n4log (n)

CTPS Thin plate spline (TPS)

x2log (x)

2.3.3. Absolute and relative radial basis function interpolation In principle there are two different ways to implement this RBF mesh motion method, the absolute and the relative implementation. The absolute method performs a direct solve of the system (8) only once at the beginning of the simulation. This method is efficient since the direct matrix solve, which is more expensive than the evaluation, is only performed initially. On the other hand, the mesh quality is limited since the coefficients are not defined with respect to the previous time-step. Therefore, the relative method is used when large boundary displacements occur, like a 180 rotation. In this method the inverse is calculated at every time-step and the motion is defined as the increment with respect to the previous time-step. For reasonably small boundary displacements, it is much more efficient to use the absolute method. When applying the relative implementation it is important to use different techniques to decrease the number of boundary points, resulting in lower computation costs (see Section 2.4). Besides the large computational costs of the relative method, there is another potential drawback of this method. Since the relative method rebuilds the interpolation matrix every time-step, it is non-linear and therefore susceptible to hysteresis. This hysteresis behaviour is especially visible at larger time-steps, leading to a progressively increasing difference in mesh quality every motion period compared to the initial mesh. However, the flow solver stability usually puts a more strict constraint on the time-step size, such that this hysteresis behaviour is of lesser importance. 2.4. Improving computational efficiency When using direct methods, the computational costs of the new mesh motion solver based on radial basis function interpolation,    increases with an increased number of boundary O N 3b and internal points ðOðN b N in ÞÞ. The method consists of two computationally expensive steps: 1. Solving the system of Eqs. (8) for given boundary points xbi and corresponding displacements Dxbi to find the coefficient vectors c and b. The upper diagonal block matrix Ubb is in general a dense symmetric matrix of size (Nb  Nb). Therefore, standard direct solvers require OðN 3b Þ operations.

Table 2 Computational costs for solving the system and evaluate the RBF’s. Computational costs for solving the system and evaluation the RBF’s on all internal mesh points. Illustration of a two-dimensional and three-dimensional uniform square shaped Cartesian mesh with an equal number of boundary nodes Nb in x- and y-direction.

2D 3D

Internal points

All boundary points

Direct solve

RBF evaluation

N 2b N 3b

4Nb

64N 3b

4N 3b

6N 2b

216N 6b

6N 5b

171

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

104

70 Mesh motion Flowe quation

60

Total timestep

Computing time [s]

Computing time [s]

103

102

101

50 40 30 20 10

100 100

101

102

0

0

1

2

4

8

Number of grid cell s(x100)

Number of grid cells (x100,000)

(a)

(b)

16

Fig. 1. Computing times for different mesh motion solvers. (a) shows a comparison of computing times for the mesh motion around a rotating two-dimensional block based on the Laplace equation (), the SBR Stress equation (.), RBF absolute method (N), RBF relative method () and RBF relative method in combination with coarsening and smoothing techniques (). The computing times for the Laplace equation () and the SBR Stress equation (.) are nearly identical. (b) shows the computing times for one timestep of a three-dimensional flapping wing simulation. The times are subdivided into the contributions of solving the RBF mesh motion and the flow equations. The total time = flow solve time + mesh motion time + additional overhead. Meshes are used from 100k to 1600k.

2. Evaluation of the radial basis function summation, Eq. (7), at all Nin (of order OðN in Þ) internal mesh points, using the given boundary points xbi and the in step 1 computed coefficient vectors c and b. This evaluation leads to a computational cost of order OðN b N in Þ. For large two- and three-dimensional meshes both the system solve and the evaluation procedures may become computationally expensive, especially when direct methods are used. In order to illustrate the increasing costs with increasing number of mesh points, a two- and three-dimensional Cartesian grid is considered with a uniform mesh distribution. The computational domain is square shaped with an equal number of points on all edges of the boundary, equal to Nb. Table 2 shows the total number of internal and boundary points as a function of Nb, for both the two- and three-dimensional example. Additionally, the total number of operations for a direct solve and the RBF evaluation is given in the table, which shows that the computational costs for both a direct solve and evaluation scales with N 3b for the two-dimensional case. Concerning the three-dimensional case, the computational costs for the direct system solve are a factor Nb larger compared to the costs for the RBF evaluation. So when a complex threedimensional case, like a flapping wing, is considered, special treatment is necessary concerning the system solve, e.g. by reducing the number of boundary nodes or using advanced direct solver techniques. 2.4.1. Boundary point coarsening and smoothing From Table 2 it is seen that the total computation costs will decrease if a constraint is put on the number of mesh points. This can be achieved in two different ways. First, a major part of the computational cost is spent at solving the system, step 1, which is a factor Nb more expensive than the evaluation. It seems reasonable to reduce the number of moving boundary points by performing a proper selection procedure. Especially when the body displacement follows a rigid body motion, not all boundary points are necessary. Therefore, a coarsening technique was incorporated which selects a boundary point for every Nc points, where 1/Nc is the coarsening factor. However, reducing the number of boundary points introduces additional errors in mesh deformation for non-rigid body

motions such that the resulting mesh quality may not be optimal. In general, we found these errors negligible. More advanced coarsening techniques, based on greedy algorithms have been applied by [28]. Secondly, the efficiency of the RBF method is improved by the notion that all outer boundary points are fixed in general. Therefore, the outer boundary points can be neglected. This is achieved by specifying a smoothing function, such that the RBF contribution reduces to zero at the outer boundary, which is defined [14] as

8 ~ 6 0; x > < 1; 2 ~ ~ ~ ~ 6 1; wðxÞ ¼ 1  x ð3  2xÞ; 0 6 x > : ~ P 1; 0; x

ð10Þ

Fig. 2. The initial quadrilateral mesh with optimal quality around a moving block. The domain size is (25D  25D) around a block with size (5D  1D), at every unit spacing, a grid point is places such that an optimal quadrilateral mesh is obtained.

172

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

Fig. 3. The cell non-orthogonality is compared of the Laplace (a) the SBR Stress (b) and the RBF (absolute) (c) motion solvers, lower (blue) is better. The non-orthogonality is visualised for a two-dimensional block with a combined motion of translation and rotation. The boundary points translate over a distance of 2.5D in both X- and Y-direction and rotate 57.3° (1.0 rad) around the translating centre. The used basis function was the globally supported thin plate spline (TPS). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where ~ x is given by:

~x ¼

kxi  Rinner k : Router  Rinner

ð11Þ

xi represents the coordinate of the ith inner mesh point, to evaluate this smoothing function at that particular location in space. Rinner and Router are two radii, between those radii this smoothing function is decreasing from 1.0 to 0.0. When the RBF evaluation function ~ Þ the corrected RBF evaluation is obtained (9) is multiplied by /ðx by:

~inj ¼ sðxinj Þ  wðx ~Þ: Dx

ð12Þ

Within Rinner the contribution of the RBF evaluation remains unaltered, but outside Router the value of the RBF becomes zero, such that all fixed outer boundary points may be neglected. In principle, the inner radius is chosen to be multiple boundary lengths (wing chords) and the outer radius is chosen as the distance from the moving boundary to the outer boundary. Therefore, the system to be solved contains the control points on the moving boundary,

selected by the coarsening function. A more extensive method was recently employed by [29] using a simple decaying correction to remove the error at all of the unselected boundary points. In addition, it is interesting to address the computing times concerning the different mesh motion solvers. Fig. 1a shows the computing times for a rotating block with increasing grid resolution for the different mesh motion solvers (Laplace, SBR Stress and RBF). Concerning the RBF mesh motion solvers, the three described variants are used, the absolute implementation, relative implementation and the relative method in combination with the coarsening and smoothing techniques. It is clear that the mesh motion solvers based on solving a partial differential equation are fast, since standard iterative techniques can be used for these sparse systems. If the absolute and relative RBF methods are compared, it is observed that these methods require large computing times, at least an order of magnitude more. However, if the coarsening and smoothing techniques are applied the computing times are again of similar order compared with the fast Laplace mesh motion. The non-linear behaviour of the curve for the latter method is caused by the choice of the coarsening ratio to select the moving boundary points. While

173

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

the mesh strategy. It will be shown that the mesh quality is significantly better in case of the RBF mesh motion compared to the other methods considered. 2.5. Mesh quality measures

Fig. 4. The cell non-orthogonality is shown for the relative versions of the RBF mesh motion solver, lower (blue) is better. The non-orthogonality is visualised for a twodimensional block with a combined motion of translation and rotation. The boundary points translate a distance of 2.5D in both X- and Y-direction.The rotation is 180° (3.14 rad) and the used basis function was the globally supported thin plate spline (TPS). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

keeping the coarsening ratio fixed, the mesh resolution is increased, the number of boundary points used in the system solving is still growing non-linearly. For high resolution meshes, the order of computing times can be decreased further by increasing the coarsening ratio. Fig. 1b shows the computing times, per time-step, concerning a three-dimensional Navier–Stokes simulation of a harmonically flapping wing. The Reynolds number is defined to be 200, such that the resulting flow around the wing was laminar. The simulations are performed using the transient Navier–Stokes solver within the OpenFOAMÒ framework, including dynamic mesh capability. Second order spatial and temporal discretization schemes are used in combination with far field and symmetry boundary conditions. These three-dimensional simulations are performed for grid sizes of 100k, 200k, 400k, 800k and 1600k cells and needed a total computing time of about 8, 32, 48, 103, 190 h, respectively. All simulations were performed on four CPU cores of an AMD OpteronÒ 280 cluster. Fig. 1b shows how the computing time increases with increasing number of grid cells for both solving the flow equations and the RBF mesh motion. Validation showed that a mesh resolution of 800k provided accurate results for flapping wing aerodynamics. Additionally, Fig. 1b shows that the computing time used for RBF mesh motion is less than 10% of the computing time for the complete time-step across different mesh sizes. In this regard, the mesh quality appears to be a more critical factor while deciding

In order to compare the quality of the different meshes after mesh motion it is important to elaborate on how to interpret the mesh quality measures [19] of skewness and non-orthogonality. These mesh quality measures are related to the cell properties such as size, orientation, shape and skewness. To compare the mesh quality of different mesh motion solvers, a Cartesian grid around a square box (optimal initial mesh quality) is generated, as considered in more detail in Section 3. It is important that the ideal mesh motion solver maintains high quality in terms of skewness and non-orthogonality after mesh deformation. The mesh skewness is a measure of the distortion of an element. Additionally, the mesh non-orthogonality is defined as the average angle between any combination of two adjacent sides of the element. Both skewness and non-orthogonality are measures for the cell deformation, such that a lower value effectively means a higher mesh quality. Therefore, desirable mesh quality bounds (lower is better) are:

0:0 6 fskewness 6 1:0; 0 6 fnonortho 6 90 :

ð13Þ

When assessing the mesh quality it is important to analyse the lowest and average values of skewness and non-orthogonality. The lowest value provides an indication if the numerical simulation will be stable and converge at all. If the worst cell quality is too low, the simulation will diverge within a couple of iterations. On the other hand, the average value of the mesh quality measure will provide an indication of the average quality of the mesh. The higher the average quality of the mesh, the more stable, accurate and efficient the computation will be. In the next section, both the average and minimal value of the skewness and non-orthogonality mesh quality measures are used to compare the mesh motion solvers. 3. Results and discussion In Section 2, three different kinds of mesh motion solvers have been described, based on solving the Laplace equation, solving the solid body rotation (SBR Stress) equation, as well as methods based on interpolation using radial basis functions (RBF). This section addresses three numerical test cases which were performed to compare the mesh quality obtained with the different mesh motion solvers. Additionally, the effect of different radial basis functions is discussed. The first test case is a two-dimensional block which performs a combination of translation and rotation. The initial computational mesh is shown in Fig. 2. The domain size of the test problem is limited to 25D  25D and the size of the moving block is

Table 3 Comparison of mesh quality for different mesh motion solvers. The mean (over a half period) and maximum values of the skewness fs and non-orthogonality fo are compared at maximal displacement and rotation of the two-dimensional rectangular block. Results are shown for the Laplace, SBR Stress and RBF mesh motion solver, the latter using different RBFs. Method

(fs)max

Laplace SBR Stress CP C2 CTPS C1

0.95 0.96 0.65 0.86 0.81

(baseline) (+0.7%) (32%) (10%) (15%)

0.88

(7%)

0.11

(+22%)

76.0

(+5%)

32.8

(+63%)

0.52

(45%)

0.051

(41%)

52.5

(27%)

19.0

(6%)

CTPS C 2a CTPS C 2b TPS

(fs)ave 0.09 0.08 0.065 0.105 0.103

(fo)max (–) (4%) (28%) (+17%) (+14%)

72.1 75.1 59.6 74.5 73.8

(fo)ave (–) (+4%) (17%) (+3%) (+2%)

20.1 21.3 23.4 31.9 31.4

(–) (+6%) (+16%) (+59%) (+56%)

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

90

90

80

80

70

70

60

60

Number of cells

Number of cells

174

50 40 30

50 40 30

20

20

10

10

0

0

10

20

30

40

50

60

70

80

0

90

0

10

20

30

40

50

70

60

Non-orthogonality

Non-orthogonality

(a) Laplace

(b) RBF

80

90

Fig. 5. Cell non-orthogonality histograms for Laplace (a) and RBF mesh motion solvers (b) using the TPS.

5D  1D, the spacing corresponds to 1D in order to obtain a Cartesian grid as can be seen in the figure. The implication of using a larger computational domain is addressed in Section 2.4. Mesh motion simulations are performed using the different mesh motion solvers and different radial basis functions. 3.1. Translation and rotation of a two-dimensional block

500

500

450

450

400

400

350

350

Number of cells

Number of cells

The first test case considers a combined motion of translation and rotation to compare the mesh motion solvers under these conditions. The two-dimensional block is initially centred and translates 2.5D in both X-and Y-direction. In addition, the block is rotated around its translating centre with 57.3° (1.0 rad). The outer boundary points are kept fixed, such that the influence of the moving boundary points on all internal points could be studied independently. The resulting mesh quality, after the combined translation and rotation, is assessed using the fields of skewness and non-orthogonality. Using those fields, the maximal and average values are obtained as well as a complete visualisation of those mesh quality measures, combined with mesh quality histograms. Before proceeding to the results, some special settings, applicable to this test case, need to be described. First, when using the Laplace and SBR Stress mesh motion solver it is necessary to specify

the diffusivity coefficient which decreases proportionally to the distance from the moving boundary points. Following [16] and [15], a quadratically decreasing diffusivity coefficient, i.e. decreasing from the moving boundary, was chosen. Secondly, the boundary conditions need to be set for the motion solver. In both cases, Laplace and SBR Stress motion solvers, the boundary conditions on all outer boundary points are set to Dirichlet type with value 0, so the outer boundary points are kept fixed. Finally, the newly implemented RBF mesh motion solver is used with five different functions, CP C2, CTPS C1, CTPS C 2a , CTPS C 2b and TPS, based on the assessment performed by [6]. Fig. 3a and b shows the cell non-orthogonality at maximal mesh deformation for the Laplace and the SBR Stress motion solvers. As seen in the figure, with these standard methods the mesh quality near the moving boundary is low, especially near the leading and trailing edges of the block. The mesh motion method which solves the Laplace equation with a quadratically decreasing diffusion coefficient, is simply not robust enough to maintain high mesh quality when the boundary rotates. As can be seen in Fig. 3a, the mesh deformation is largest near the boundary, which is not desirable. In cases where large rotation angles occur, it is best to apply a mesh motion solver, which leads to the motion of all internal mesh points, coping with the boundary deformation. The gain in mesh

300 250 200

300 250 200

150

150

100

100

50

50

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0 0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

Skewness

Skewness

(a) Laplace

(b) RBF

Fig. 6. Cell skewness histograms for Laplace (a) and RBF mesh motion solvers (b) using the TPS.

0.7

0.8

0.9

1

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

175

Fig. 7. Cell non-orthogonality for a three-dimensional wing with an ellipsoidal planform and harmonic kinematics. Results obtained with the RBF mesh motion solver with the globally supported thin plate spline (TPS), lower (blue) is better. (a) shows the mesh quality at t = 0.25T and (b) at t = 0.75T, respectively during the downstroke and upstroke. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

quality by using the SBR Stress method is marginal. From Fig. 3b it can be seen that the cells with a high non-orthogonality occur at a distance from the moving boundary. Still, the mesh quality near the body surface, especially near the leading and trailing edges, needs improvement. Improvement can be obtained by specifying a constraint to an inner mesh region, such that it moves according to the body motion. Fig. 3c shows the mesh obtained using RBF interpolation using a thin plate spline function. When comparing the cell non-orthogonality with Fig. 3a and b it may be seen that most of the cell deformation, with the RBF mesh motion, occurs in the outer regions of the mesh and all cells are dealing with the boundary displacement. Additionally, Fig. 4 shows the resulting mesh, which is obtained using the relative implementation of the RBF mesh motion after a rotation of 180 deg. It shows that even for this excessive rotation, the mesh remains valid. This method is robust, yet computationally expensive as will be shown in Section 2.4.

Table 3 shows a quantitative comparison of the resulting mesh quality obtained with different mesh motion solvers, including different RBF’s. The maximum and averaged (over a half flapping period) values of both skewness and non-orthogonality are compared. The mesh quality calculated by the Laplace method is low, since the maximal skewness and maximal non-orthogonality are large, respectively fsmax ¼ 0:95 and fomax ¼ 72:1. The results of the Laplace and SBR Stress mesh motion solver are similar, within 6.1%. Concerning the RBF method, four RBF’s with compact support are used (CP C2, CTPS C1, CTPS C 2a and CTPS C 2b ) and one with global support, the thin plate spline TPS. In order to eliminate the effect of the support radius, it is set to r = 75 which is about three times the domain size such that the results in Table 3 are independent of the support radius for r > 75. The RBF method provides high mesh quality, maximal and averaged, for both C2 and TPS compared to the other functions. The C2 and TPS RBF’s result in a maximum skewness of respectively 32% and 45% compared to the Laplace motion

Fig. 8. Deforming mesh around a flexible rectangular block which is translating, rotating and flexing using the RBF mesh motion solver with the globally supported thin plate spline (TPS). The flexing is defined by simple harmonic functions. (a) shows the mesh deformation at t = 0.25T and (b) at t = 0.75T.

176

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

solver, while the difference of the other RBF’s is only about 10%. Similar results are shown in the table for the average skewness and maximum orthogonality. Concerning the average orthogonality, the C2 RBF is outperformed by the TPS RBF, which is the only mesh motion solver resulting in a lower value compared to the Laplace method. Overall, the basis function thin plate spline provides the highest mesh quality for both skewness and non-orthogonality, such that this globally supported function was used for the current investigations. Finally, in addition to the non-orthogonality visualisations, Figs. 5 and 6 show the histograms of cell non-orthogonality and skewness for the Laplace and RBF mesh motion solver, which are considered to result in the lowest and highest mesh quality, respectively. For both non-orthogonality and skewness, the RBF mesh motion results in a gradual shape, due to the fact that more distant internal cells are coping with the mesh motion compared to the Laplace mesh motion.

and df represents the direction vector of the flexing. In this model problem of a moving two-dimensional block, the amplitudes of both translation and rotation were fixed to At = (2.5, 2.5, 0.0) and Aa = (0.0, 0.0, 1.0), respectively. The flexing was defined such that the main flexing direction is perpendicular to the flexible boundary surface, df = (0.0, 1.0, 0.0) with a flexing amplitude of Af = 0.5 which is about 10% of the flexible boundary length, cf = 5.0. Fig. 8 shows the resulting mesh deformation at t = 0.25T and t = 0.75T where T = 1/f is the motion period. It can be seen that the whole mesh is deformed by the RBF mesh motion solver, like was the case with a rigid airfoil. Still, some high non-orthogonality can be observed within a region of 1  2 block lengths, which is mainly caused by the fixed points on the outer boundary and the small computational domain. A larger domain will undoubtedly lead to high quality meshes when using RBF mesh motion.

4. Conclusions 3.2. Flapping of a three-dimensional wing It was shown that high mesh quality was obtained for a simplified two-dimensional test case, by using radial basis function interpolation. Especially, the globally supported thin plate spline (TPS) provided high quality and robust mesh deformation. To show the three-dimensional capabilities of the RBF mesh motion solver, the mesh around a flapping wing is shown in Fig. 7 at mid-stroke. The planform of the wing was ellipsoidal and the wing kinematics was chosen as harmonic. Fig. 7a shows the non-orthogonality during the downstroke, while Fig. 7b presents the mesh quality halfway of the upstroke. The TPS was used as radial basis function, without any user input, since it has global support. From the figure it is clear that a large part of the near wake is deformed in order to deal with the three-dimensional wing motion. Concerning the RBF mesh motion solver, the cells close to the wing take a larger part of the deformation compared to the Laplace method. It should be noted that the numerical computation of the Navier–Stokes equations did not converge using the Laplace method, due to the limited mesh accuracy close to the wing. When comparing Fig. 7a and b it is seen that the mesh during both the upstroke and downstroke is symmetric. This is caused by the fact that the radial basis function interpolation determines the new internal mesh points with respect to the initial mesh, such that the initial mesh is recovered after every flapping period. 3.3. Flexing of a two-dimensional block In extension to a rigid wing approach, a flexible wing was incorporated by prescribing the wing flexing using harmonic functions, to mimic realistic insect wing deformation [30]. In order to show that the RBF mesh motion method is able to deal with a flexing boundary, the model problem of a moving block is used. The motion of the two-dimensional block can be decomposed into translation, rotation and flexing, all defined with respect to the initial configuration. The translation and rotation are defined by:

xt ¼ At  sinð2pftÞ;

a ¼ Aa  sinð2pftÞ; where, At and Aa represent the translation and rotation amplitude vectors, f the frequency and t the time. The flexing of the boundary is defined by combining two harmonic functions like:



x  dc  sinð2pftÞ  df ; xf ¼ Af  cos 2p 2  cf where Af is the flexing amplitude vector, dc is the direction vector in-plane of the flexing surface, cf the length of the flexing surface

In this paper, a comparison was made of different mesh motion techniques for flapping wing simulations. Firstly, two methods were described which are commonly used within the code of OpenFOAMÒ, both based on solving a partial differential equation. The first method solves the Laplace equation with a variable diffusion coefficient, which is used to control the final mesh quality. In the second method, the linear stress equation was modified to include rigid body rotations in order to cope with the severe mesh deformation present in flapping wing simulations. As with the Laplace equation, the solid body rotation stress (SBR Stress) mesh motion uses the diffusivity, acting as a stiffness, to control the quality of the mesh. The diffusivity, in both cases, is defined to decrease quadratically with the distance from the moving boundary. Besides solving a partial differential equation, the motion of mesh points can be defined using interpolation techniques. A new mesh motion solver is incorporated in OpenFOAMÒ, which uses the interpolation of radial basis functions (RBF). For given boundary point displacements the internal mesh displacements are obtained by solving a system of equations to obtain an array of interpolation coefficients. Using those coefficients, the internal point displacements are obtained by evaluating the radial basis functions. This new mesh motion technique does not need any information about the mesh connectivity and can be applied to arbitrary unstructured meshes containing polyhedral cells, which is the way OpenFOAMÒ deals with the finite volume implementation. The three mesh motion solvers are tested on a two-dimensional rectangular block which moves through a Cartesian mesh. The cell non-orthogonality and skewness are compared. Additionally, different radial basis functions, concerning the RBF mesh motion, are compared. The RBF mesh quality provides superior mesh quality over the Laplace and SBR Stress mesh motion solvers. Especially when using the thin plate spline (TPS) or the continuous polynomial C2 as radial basis functions, the mesh quality is high in terms of low skewness and non-orthogonality. Note that the German Aerospace Centre (DLR) uses the TPS as well in the TAU code [12]. The TPS has global support, whereas the C2 basis function has compact support. The RBF mesh motion was successfully tested on simple test problems and for a three-dimensional flapping wing with the possibility to incorporate a flexing moving boundary. Since the RBF mesh motion technique encounters a dense system of equations, different methods are implemented to increase its efficiency. First of all, a subset of the moving boundary points was selected, because not all points are necessary when the body performs a rigid body motion. So a coarsening algorithm selects those control points. Secondly, a smoothing function is used to

F.M. Bos et al. / Computers & Fluids 79 (2013) 167–177

decrease the RBF contribution to zero at the outer domain boundaries. Therefore, it is justified to neglect the outer (fixed) boundary points, which reduces the system of equations considerably. From this discussion it is concluded that, concerning the three-dimensional wing simulations, the globally supported TPS should be used in combination with the coarsening and smoothing techniques to increase the efficiency of the RBF mesh motion method. Acknowledgement This research has been supported by the Dutch Organisation for Scientific Research, NWO-ALW Grant 814.02.019. References [1] Aono H, Liang F, Liu H. Near- and far-field aerodynamics in insect hovering flight: an integrated computational study. J Exp Biol 2008;211:239–57. [2] Baruh H. Analytical dynamics. McGraw-Hill Inc.; 1999. [3] Batina JT. Unsteady euler airfoil solutions using unstructured dynamic meshes. AIAA J 1990;28(8):1381–8. [4] Beckert A, Wendland H. Multivariate interpolation for fluid–structureinteraction problems using radial basis functions. Aerospace Sci Technol 2001;5(2):125–34. [5] Blom FJ. Considerations on the spring analogy. Int J Numer Methods Fluids 2000;32:647–69. [6] de Boer A, van der Schoot MS, Bijl H. Mesh deformation based on radial basis function interpolation. Comput Struct 2007;85:784–95. [7] Degand C, Farhat C. A three-dimensional torsional spring analogy method for unstructured dynamic meshes. Comput Struct 2002;80:305–16. [8] Donea J. An arbitrary Lagrangian–Eulerian finite element method for transient fluid–structure interactions. Comput Methods Appl Mech Eng 1982:689–723. [9] Dwight RP. Robust mesh deformation using the linear elasticity equations. In: Deconinck H, Dick E, editors. Computational fluid dynamics 2006. Springer; 2004. [10] Farhat C, Degand C, Koobus B, Lesoinne M. Torsional springs for twodimensional dynamic unstructured meshes. Comput Methods Appl Mech Eng 1998;163:231–45. [11] Ferziger JH, Peric M. Computational methods for fluid dynamics. 3rd ed. Berlin: Springer-Verlag; 2002. [12] Heinrich R, Kroll N. Fluid–structure coupling for aerodynamic analysis and design: a dlr perspective. In: 46th AIAA aerospace sciences meeting and exhibit, Reno, 2008-561; 2008. [13] Helenbrook BT. Mesh deformation using the biharmonic operator. Int J Numer Methods Eng 2003;56(7):1007–21. [14] Jakobsson S, Amoignon O. Mesh deformation using radial basis functions for gradient-based aerodynamics shape optimization. Comput Fluids 2007;36: 1119–36. [15] Jasak H. Dynamic mesh handling in openfoam. In: 47th AIAA aerospace sciences meeting, Orlando, 2009-341; 2009.

177

[16] Jasak H, Tukovic´ Z. Automatic mesh motion for the unstructured finite volume method. Transactions of FAMENA 2007; 30(2). [17] Jasak H, Weller HG, Nordin N. In-cylinder CFD simulation using a c++ objectoriented toolkit. SAE technical paper 2004-01-0110; 2004. [18] Johnson AA, Tezduyar TE. Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces. Comput Methods Appl Mech Eng 1994;119:73–94. [19] Knupp PM. Algebraic mesh quality metrics for unstructured initial meshes. Finite Elem Anal Des 2003;39:217–41. [20] Lee J, Kim J, Choi H, K. Y. Sources of spurious force oscillations from an immersed boundary method for moving-body problems. J Comput Phys 2011;230:2677–95. [21] Liu X, Qin N, Xia H. Fast dynamic grid deformation based on delaunay graph mapping. J Comput Phys 2006;211:405–23. [22] Löhner R, Yang C. Improved ale mesh velocities for moving bodies. Commun Numer Methods Eng 1996;12(10):599–608. [23] Miller L, Peskin C. Flexible clap and fling in tiny insect flight. J Exp Biol 2009;212:3076–90. [24] Mittal R, Iaccarino G. Immersed boundary methods. Annu Rev Fluid Mech 2005;37:239–61. [25] Peskin CS. The immersed boundary method. Acta Numer 2002:479–517. [26] Rendall TCS, Allen CB. Improved radial basis function fluid–structure coupling via efficient localised implementation. Int J Numer Methods Eng 2008;78(10): 1188–208. [27] Rendall TCS, Allen CB. Unified fluid–structure interpolation and mesh motion using radial basis functions. Int J Numer Methods Eng 2008; 74(10):1519–59. [28] Rendall TCS, Allen CB. Efficient mesh motion using radial basis functions with data reduction algorithms. J Comput Phys 2009;228(17):6231–49. [29] Rendall TCS, Allen CB. Parallel efficient mesh motion using radial basis functions with application to multi-bladed rotors. Int J Numer Methods Eng 2010;81(1):89–105. [30] Shyy W, Lian Y, Tang J, Liu H, Trizila P, Stanford B, et al. Computational aerodynamics of low reynolds number plunging, pitching and flexible wings for mav applications. Acta Mech Sin 2008;24:351–73. [31] Su SW, Lai MC, Lin CA. An immersed boundary technique for simulating complex flows with rigid boundary. Comput Fluids 2007;36:313–24. [32] Sun M, Lan S. A computational study of the aerodynamic force and power requirements of dragonfly hovering. J Exp Biol 2004;207:1887–901. [33] Sun M, Yu M. Flows around two airfoils performing fling and subsequent translation and subsequent clap. Acta Mech Sin 2003;19(2):103–17. [34] Wang ZJ. Vortex shedding and frequency selection in flapping flight. J Fluid Mech 2000;410:323–41. [35] Wang ZJ, Przekwas AJ. Unsteady flow computation using moving grid with mesh enrichment. Tech. rep., AIAA-94-0285; 1994. [36] Weller HG, Tabor G, Jasak H, Fureby C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput Phys 1998;12(6). [37] Wendland H. Konstruktion und untersuchung radialer basisfunktionen mit kompaktem träger. Tech. rep.; 1996 [38] Zhu L, Peskin CS. Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method. J Comput Phys 2002;179:452–68.