Numerical Simulation of the Wake Structure and Thrust/Lift Generation of a Pitching Airfoil at Low Reynolds Number Via an Immersed Boundary Method

ABSTRACT:

In this study, an accurate computational algorithm in the context of immersed boundary methods is developed and used to analyze an incompressible flow around a pitching symmetric airfoil at Reynolds number (Re = 255). The boundary conditions are accurately implemented by an iterative procedure applied at each time step, and the pressure is also updated simultaneously. Flow phenomena, observed at different oscillation frequencies and amplitudes, are numerically modeled, and the physics behind the associated vortex dynamics is explained. It is shown that there are four flow regimes associated with four wake structures. These include three symmetric flow regimes, with adverse, favorable and no vortex effects, and an asymmetric flow regime. The phenomena associated with these flow regimes are discussed, and the critical or transitional values of the Strouhal (St) and normalized amplitude (A_D) numbers are presented. It is shown that, at the fixed pitching amplitude, A_D = 0.71, the transition from adverse (drag generation) to favorable (thrust generation) symmetric flow regime occurs at St = 0.23. Moreover, at this particular amplitude, transition from symmetric to asymmetric regime occurs at St = 0.48. It is also shown that, at St = 0.22, the wake is always deflected and the flow is asymmetric for large enough amplitudes A_D > 2. The dipole vortices and lift generation are two characteristics of asymmetric vortex street. This numerical study also reveals that the initial phase angle has a dominant effect on the appearance of dipole vortices and vortex sheet deflection direction. Numerical results are in good agreement with the available experimental data.

KEYWORDS:
Immersed boundary method; Flapping airfoil; Karman vortex street; Thrust generation; Lift generation

INTRODUCTION

Birds, insects and marine creatures use, and benefit from, the phenomena associated with flapping wings, tails and fins. Inspired by the nature, scientists and engineers have focused on the performance of oscillating wings and airfoils and have found that flapping mechanism is more efficient compared to the classical fixed wing designs at low Reynolds numbers (Re) (Davis 2007).

The flow field around a flapping airfoil is very complex due to the discontinuous production of vortices and transient interactions between them. Studies on pitching airfoils have shown that there are different flow regimes associated with different wake patterns. For example, Kármán vortex street (KVS) and reverse Kármán vortex street (RKVS) are two different flow patterns that result in drag enhancement and thrust generation, respectively. In addition, the vortices in the wake region of a pitching airfoil can be symmetric or deflected, and this also affects the momentum exchange between the airfoil body and the surrounding fluid (Ashraf 2010). It is, therefore, very important to understand the development of these flow regimes, the transition from one regime to another and to quantify the favorable and adverse effects associated with flapping airfoils.

The transient flow field around flapping airfoils has been studied theoretically, experimentally and numerically for further understanding of the effects of parameters such as frequency, amplitude and flapping mode on the lift/thrust generation and momentum loss.

On the theoretical side, Jones et al. (1998) have used the effective angle of attack to explain the thrust generation by a flapping airfoil. The role of the vortex dynamics on thrust and drag generation, or momentum surfeit or deficit in the wake region, is explained by von Kármán and Burgers (1934). Also, Weis-Fogh (1973) explains how insects fly using the unsteady lift generation mechanism.

Many experimental studies have also been carried out. Anderson et al. (1998) have provided experimental results related to thrust generation through a combined pitching and plunging motion of a NACA0012 airfoil in a water tunnel. Based on this study, the optimum thrust is obtained at a Strouhal number (St) between 0.25 and 0.4. Taylor et al. (2003) have experimentally studied wing frequencies of 42 different birds, bats and insects and concluded that the Strouhal numbers associated with these natural flyers are between 0.19 and 0.41. These results are in agreement with the results reported by Anderson et al. (1998). Lentik et al. 2007 have experimentally studied the vortex dynamics of a flapping airfoil in a soap film tunnel at Re = 1,000 and reported that, at high amplitude oscillations, asymmetric wakes are generated when the frequency is increased. Godoy-Diana et al. (2008) have experimentally studied the vortex structure around a pitching airfoil and presented the results in the form of a parameter map for different wake types. Bohl and Koochesfahani (2009) have studied the wake patterns behind a sinusoidally pitching symmetric airfoil at various oscillation frequencies and concluded that the wake pattern can be controlled by adjusting the frequency, amplitude and the oscillation mode (shape of the oscillation wave form). They have also reported the transition point at which the KVS regime turns to the RKVS one. Zanotti et al. (2014) studied, both experimentally and numerically, a pitching airfoil in deep dynamic stall condition. They reported that, during upstroke motion, 3-D numerical models are in better agreement with the experiments as compared to 2-D models.

Considering the complexity of the vortex dynamics around a flapping airfoil and the continuous progress in computational algorithms and computer software and hardware, numerical simulation is becoming more and more popular. Wang (2000) has numerically calculated the flow around a 2-D symmetric airfoil undergoing pure plunging motion by an incompressible Navier-Stokes solver at Re = 1,000 and reported that the flow separates from both the leading and trailing edges at this particular Re. Lewin and Haj-Hariri (2003), as well as Young and Lai (2007), have provided numerical results on wake structures and flow regimes around plunging airfoils. According to the latter study, aerodynamic forces are more strongly affected by the leading edge vortices while the wake structure is mainly controlled by trailing edge vortices. Schnipper et al. (2009) have provided information regarding the effects of the wake pattern on aerodynamic forces for a flapping airfoil. He et al. (2012) have studied transition and symmetry-breaking in the wake of a flapping airfoil by an immersed boundary method (IBM). They have presented results that show flow patterns from the KVS regime to the RKVS one. Yu et al. (2012) investigated the asymmetric wake vortex structure around an oscillating airfoil both numerically and experimentally. They reported that the deflected wake appears around St = 0.31 for pitching amplitude equal to 5º. They also reported that, as the St increases, the dipole mode of the vortex pair becomes more prevalent. Dipole vortex is necessary for the formation of asymmetric wake. Ren et al. (2013) reported that the direction of wake asymmetry changes at St = 0.37 for a fixed amplitude corresponding to α = 5º. Khalid et al. (2014) have numerically simulated the equivalence of pitching and plunging motions found in a flapping NACA0012 airfoil. They have reported that wake deflection is observed, though not dominant, at low Strouhal numbers. Bai et al. (2014) used an immersed boundary method to simulate the turbulent flow around a horizontal axis turbine operating under free surface waves. The IBM which is implemented in this investigation uses a 3-D finite volume solver. Also, Wei and Zheng (2014) have implemented an IBM to study the wake downstream of a 2-D heaving airfoil.

In this study, the physics behind the wake structure, the transition from one flow regime to another and the thrust/lift performance of a sinusoidally flapping symmetric airfoil over a wide range of Strouhal numbers and oscillation amplitudes is studied. Furthermore, the effect of the initial phase angle on the deflected vortex street, not sufficiently discussed in the relevant literature, is studied. To carry out this investigation, an accurate computer code based on the IBM is developed and used to simulate the flow. IBM is particularly well suited for moving boundary problems including the problem of flow around a flapping airfoil.

The proposed solution method is a finite difference iterative algorithm, in which both pressure and forces at the Lagrangian (boundary) points are updated simultaneously. Therefore, boundary conditions are implemented by updating the body force and pressure iteratively in an inner loop at each time step. Two test cases are used to validate the computational results.

The present numerical study is carried out at a constant Reynolds number (Re = 255). The airfoil is also selected such that the experimental results reported by Godoy-Diana et al. (2008) are numerically modeled. Furthermore, detailed information regarding the vortex shedding, time-averaged velocity field, and aerodynamic coefficients at different oscillation frequencies and amplitudes is reported. The phenomenon of transition from drag-generating wake to thrust-generating wake is also discussed.

IMMERSED BOUNDARY METHOD

GENERAL VIEW

In this section, some background information relevant to the commonly used IBMs is presented.

IBMs use a background uniform Cartesian grid, which does not generally conform to the boundary of the flow domain. Boundary conditions are then imposed by source terms in the governing equations discretized on the background Cartesian mesh. In moving boundary problems, the displacement of boundary nodes corresponds to the shift of source terms between fixed grid points.

In IBMs, one has to distinguish between two different types of nodal points. Fixed Cartesian grid points are called the Eulerian points and the possibly moving boundary points are called the Lagrangian points. To numerically solve an incompressible flow around a moving object via IBM, it is assumed that the boundary imposes force on the fluid and the fluid imposes an equal and opposite force on the boundary. In classical Peskin's method, the immersed boundary is represented by a set of elastic fibers whose locations are tracked by Lagrangian points (Peskin 2002). Figure 1 shows the coordinates of Lagrangian and Eulerian points and the zone of influence for the force term at Lagrangian point K.

Figure 1
(a) Lagrangian (k (Lima e Silva et al. 2003).) points; (b) zone of influence for the force term at Lagrangian point ) and Eulerian (

The force terms on Lagrangian points, represented by , are related to the boundary force as follows:, are then related by the Hook's law or a similar relation to the displacement of these points. These force terms are distributed on Eulerian grid points by a proper delta-type function. Therefore, boundary force at each Lagrangian point is distributed over a group of cells around that Lagrangian point. The forces exerted at the background grid points,

(1) $f\left(\overrightarrow{x},t\right)={\int }_{{\Gamma }_b}\overrightarrow{\mathbf{F}}\left(s,t\right)\delta \left(\overrightarrow{x}-\overrightarrow{\mathbf{X}}\left(s,t\right)\right)\mathrm{ds}$

where:

t is time; Γ_b is the boundary of the solid domain; s is the body-fitted coordinate along the boundary; ∆ is a smooth Dirac delta function; = (x, y) are the Cartesian coordinates of an Eulerian point; = (X, Y) are the Cartesian coordinates of a Lagrangian point.

The distributed forces on Eulerian points are used in the discrete governing equations as source terms. These equations are then solved to calculate fluid velocities at Eulerian grid points, represented by . The velocities at boundary Lagrangian points, U, are calculated afterwards as follows:

(2) $\overrightarrow{\mathbf{U}}\left(\overrightarrow{\mathbf{X}}\left(s,t\right)\right)={\int }_{\overrightarrow{x}\varepsilon {\Omega }_f+{\Omega }_S}\overrightarrow{u}\left(\overrightarrow{x},t\right)\delta \left(\overrightarrow{x}-\overrightarrow{\mathbf{X}}\left(s,t\right)\right)d\overrightarrow{x}$

where:

Ω_s is the solid domain and Ω_f is the fluid domain.

Note that discrete governing equations are solved at all Eulerian grid points regardless of the positions of those points. In other words, non-physical velocities are assigned to the Eulerian grid points outside of the flow domain. These numerical results are some by-products of the IBM computations.

The method just described is called the continuous forcing method and is often used to solve fluid-structure interaction problems with elastic boundary at low Reynolds numbers.

The direct forcing method, developed by Mohd-Yusof (1997) and further refined by Tseng and Ferziger (2003), Balaras (2004) and Mittal et al. (2008), is another approach in the family of immersed boundary techniques. This method is better suited for high Re flows. The force is implemented into the momentum equation directly by substituting the regularized no-slip condition near immersed boundary.

Exact implementation of the boundary conditions as well as the satisfaction of conservation laws are two major concerns in all IBMs. To satisfy mass and momentum conservation near the boundary, the cut cell method (Chung 2006) has also been developed. In this method, finite volumes are defined so that immersed boundary coincides with boundary cell faces. This method and some later improvements, e.g. Ls-STAG (Cheny and Botella 2010), are rather more complex compared to the two previously mentioned methods.

In this paper, an iterative-direct forcing IBM is employed. Variables at Lagrangian and Eulerian grid points are linked through a simple discrete delta function. Details of the proposed method are given in the next section.

THE PROPOSED COMPUTATIONAL ALGORITHM

Formulation and Force Term Calculations

In the formulation presented here, which has been proposed by Wang et al. (2008), a first-order time discretization similar to Uhlmann's method (Uhlmann 2005) is used to implement the Lagrangian force terms. The governing equations for incompressible flow in the entire computational domain (Ω_f + Ω_g) are:

(3) $\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}·▽\overrightarrow{u}=-\frac{1}{\mathit{\rho }}▽P+v{▽}^2\overrightarrow{u}+\overrightarrow{f}$

(4) $▽·\overrightarrow{u}=0$

(5) $\overrightarrow{u}={\overrightarrow{U}}_{\Gamma }\mathrm{on}{\Gamma }_b$

where:

ρ is the fluid density; ν is the kinematic viscosity, defined as ν = μ / ρ, where μ is the fluid dynamic viscosity; V and V² are the regular centered difference approximations for the gradient and Laplace operators, respectively. The source term Eq. 3 is a fraction of the relevant external boundary force _Γ at the immersed boundary. This is done to assure that the velocities at Lagrangian points satisfy the no-slip boundary condition. is the velocity vector on the Eulerian points; P is the pressure; in that acts at Eulerian grid point. The force term is imposed at the Lagrangian point to enforce the desired velocity

An estimate of the velocity can be obtained by the following equation:

(6) $\frac{\overrightarrow{\hat{u}}-{\overrightarrow{u}}^n}{\Delta t}+{\left(\overrightarrow{u}·▽\overrightarrow{u}\right)}^n=-\frac{1}{\mathit{\rho }}▽P^n+v{\left({▽}^2\overrightarrow{u}\right)}^n$

It is important to note that is an intermediate velocity vector at Eulerian points that satisfies the momentum equation without considering the effects of the boundary conditions. The superscript n indicates the last time step, and ∆t is the time step. To calculate nodal velocity values that satisfy both the boundary conditions and governing equations, the effect of the force term should be added as follows:

(7) ${\overrightarrow{u}}^{n+1}=\overrightarrow{\hat{u}}+\overrightarrow{f}\Delta t$

Rewriting Eq. 7 at a Lagrangian point, it results in:

(8) ${\overrightarrow{U}}^{n+1}=\overrightarrow{\hat{U}}+\overrightarrow{F}\Delta t$

where:

ⁿ⁺¹ is the velocity at the Lagrangian point at time level n + 1; U is an intermediate Lagrangian point velocity calculated as follows:

(9) $\overrightarrow{\hat{U}}\left(\overrightarrow{X}\left(s,t\right)\right)=\int \overrightarrow{\hat{u}}\left(\overrightarrow{x,t}\right)\mathit{\delta }\left(\overrightarrow{x}-\overrightarrow{X}\left(s,t\right)\right)d\overrightarrow{x}$

The velocity term ⁿ⁺¹ in Eq. 8 can now be replaced by _Γb). This results in the following expression:

(10) $\overrightarrow{F}=\frac{{\overrightarrow{U}}_{{\Gamma }_b}-\overrightarrow{\hat{U}}}{\Delta t}$

Force terms at Eulerian points can now be calculated using the discrete delta function:

(11) $\overrightarrow{f}\left(\overrightarrow{x,t}\right)=\sum _{k=1}^N\overrightarrow{F}\left({\overrightarrow{X}}_k\right)·\mathit{\delta }\left(\overrightarrow{x}-{\overrightarrow{X}}_k\right)\Delta s_k$

where:

N is the number of Lagrangian points; ∆S_k is the distance between k^th and (k - 1)^th Lagrangian points.

Numerical Scheme

A finite difference method on a staggered grid is used in this study. The immersed boundary is discretized by N Lagrangian markers, _k = (X_k, Y_k), where the marker spacing is equal to ∆S = Γ_b / N. The following discrete Dirac delta function is used to transfer the information between Lagrangian and Eulerian points:

(12) $\mathit{\delta }\left(\overrightarrow{x}-\overrightarrow{X}\right)=d_h\left(x-X_k\right)d_h\left(y-Y_k\right)$

where:

x and y are Cartesian coordinates at Eulerian point.

The hat function, d_h, which is equivalent to the quadratic interpolation scheme (Peskin 2002), is defined as follows:

(13) $d_h\left(r\right)=\left\{\begin{array}{cc}\frac{1}{8}\left(3-2\frac{\left|r\right|}{h}\right.+\sqrt{1+4\frac{\left|r\right|}{h}-4{\left(\frac{\left|r\right|}{h}\right)}^2}\hfill & \left|r\right|\le h\hfill \\ \frac{1}{8}\left(5-2\frac{\left|r\right|}{h}\right.+\sqrt{-7+12\frac{\left|r\right|}{h}-4{\left(\frac{\left|r\right|}{h}\right)}^2}\hfill & h\le \left|r\right|\le 2h\hfill \\ 0\hfill & \left|r\right|\ge 2h\hfill \end{array}\right\}$

where:

h is the grid size and r is a variable like (x - X_k) or (y - Y_k) in Eq. 12.

To solve the Navier-Stokes equations, the fractional step method proposed by Kim and Moin (1985) is implemented. The non-linear convection term is treated by second-order Adams-Bashforth method, and Crank-Nicolson method is used for the discretization of the diffusion terms. All the space derivatives are approximated by the second-order central difference method. The time advancement and calculation of flow field at time level n + 1 is carried out through the following steps. The first Eq. 14. (intermediate velocity at Eulerian points) is calculated by

(14) $\begin{array}{c}\frac{\overrightarrow{\tilde{u}}-{\overrightarrow{u}}^n}{\Delta t}=-\frac{3}{2}{\left(\overrightarrow{u}·▽\overrightarrow{u}\right)}^n+\frac{1}{2}{\left(\overrightarrow{u}·▽\overrightarrow{u}\right)}^{n-1}\hfill \\ \frac{1}{2}v{▽}^2{\left(\overrightarrow{u}\right)}^n+\frac{1}{2}v{▽}^2{\left(\overrightarrow{u}\right)}^{n-1}\hfill \end{array}$

Then, an inner loop is used to update pressure and impose velocity boundary condition in an iterative procedure. The following equations are used at the m^th inner loop iteration:

(15) $\overrightarrow{\hat{u}}=\overrightarrow{\tilde{u}}-\frac{\Delta t}{\mathit{\rho }}▽P^{n,m-1}$

(16) $\overrightarrow{\hat{U}}\left({\overrightarrow{X}}_k\right)=\sum _{\overrightarrow{X}\epsilon \Omega }\overrightarrow{\hat{u}}\mathit{\delta }\left(\overrightarrow{x}-{\overrightarrow{X}}_k\right)h^2$

(17) ${\overrightarrow{F}}^{n,m}\left({\overrightarrow{X}}_k\right)=\frac{{\overrightarrow{U}}_{\Gamma }-\overrightarrow{\hat{U}}\left({\overrightarrow{X}}_k\right)}{\Delta t}$

(18) ${\overrightarrow{f}}^{n,m}\left(\overrightarrow{x}\right)=\sum {\overrightarrow{F}}^{n,m}\left({\overrightarrow{X}}_k\right)\mathit{\delta }\left(\overrightarrow{x}-{\overrightarrow{X}}_k\right)\Delta s_k$

(19) ${\overrightarrow{u}}^{*}\left(\overrightarrow{x}\right)=\overrightarrow{\hat{u}}\left(\overrightarrow{x}\right)+{\overrightarrow{f}}^{n,m}\left(\overrightarrow{x}\right)\Delta t$

(20) $\overrightarrow{\tilde{u}}={\overrightarrow{u}}^{*}\left(\overrightarrow{x}\right)+\frac{\Delta t}{\mathit{\rho }}▽P^{n,m-1}$

(21) ${▽}^2{P}^{n,m}=\frac{\mathit{\rho }▽\overrightarrow{\tilde{u}}\left(\overrightarrow{x}\right)}{\Delta t}$

(22) $\max \left|{\overrightarrow{F}}^{n,m}\left({\overrightarrow{X}}_k\right)\right|\le \epsilon$

where:

()are Eulerian points., ,

Equation 22, with error tolerance ε = 0.001, is checked in each inner loop to impose the velocity boundary condition exactly. At the end of inner loop, the velocity and pressure are updated as follows:

(23) $P^{n+1}=P^{n,m}$

(24) ${\overrightarrow{u}}^{n+1}\left(\overrightarrow{x}\right)=\overrightarrow{\tilde{u}}\left(\overrightarrow{x}\right)-\frac{\Delta t}{\mathit{\rho }}▽P^{n+1}$

Vector quantities, , are the intermediate velocities between the time levels n and n + 1.* and ,

The computational algorithm can now be summarized as follows:

1. Obtain the intermediate velocity at grid points, Eq. 14.), by solving (

   Start the inner loop:

2. Obtain the intermediate velocity at grid points, Eq. 15.), by solving (

3. Calculate the intermediate velocities at Lagrangian points, _k), using Eq. 16. (

4. Obtain the force terms at Lagrangian points, ^n,m (_k), by Eq. 17.

5. Distribute the Lagrangian forces among the Eulerian points using Eq. 18.

6. Correct the intermediate velocities and obtain Eq. 19.* () by

7. Obtain u by Eq. 20.

8. Update pressure by solving Poisson equation - Eq. 21.

9. Check Eq. 22 for convergence; if the convergence is achieved, then go to step 10; otherwise, go to step 2.

   End of inner loop:

10. Update pressure and velocity by Eqs. 23 and 24.

11. Go to the next time step.

VALIDATION STUDIES

Two test cases are used to validate the proposed numerical method. The first test case is a laminar flow around a fixed circular cylinder which is steady at Re = 40 and unsteady at Re = 100. The calculated drag coefficients at various Reynolds numbers are compared to the values reported by Park et al. (1998) and Lima e Silva et al. (2003), shown in Table 1.

Table 1
Comparison between the drag coefficients computed by the present method and two reported results.

Re	Present method	Lima e Silva et al. (2003)	Park et al. (1998)
10	2.79	2.81	2.78
20	2.02	2.04	2.01
40	1.52	1.54	1.51
80	1.37	1.40	1.35
100	1.35	1.39	1.33

Also, the distributions of the pressure coefficient, C_P = (P - P_∞)/(1/2ρU_∞²), along the surface of circular cylinder at Re = 40 and Re = 100, are investigated. Here, U_∞ is the far-field free stream pressure which is assumed zero and is the free stream velocity. Figure 2 shows that the computational results in this study agree well with the results of Park et al. (1998). In the case of unsteady flow (at Re = 100), pressure coefficients are averaged in time.

Figure 2
Comparison of distribution of the pressure coefficient Cp along the cylinder surface for flow past a stationary circular cylinder at Re = 40 and 100.

A comparison between the streamlines near the cylinder is shown in Fig. 3. In some of the conventional direct forcing methods, e.g. the method proposed by Su et al. (2007), the streamlines cross the boundary. This indicates that the boundary condition is not fully satisfied. In the present method, iterativedirect forcing is employed, and the boundary conditions are fully satisfied by updating the pressure and body force simultaneously in an iterative procedure.

Figure 3
Streamline around a circular cylinder. (a) Present numerical method at Re = 40; (b) Conventional method (like Su et al. 2007) at Re = 40; (c) Present numerical method at Re = 100; (d) Conventional method (like Su et al. 2007) at Re = 100.

The second test case is a laminar flow around an inline (back and forth in the x direction) oscillating a circular cylinder in a stagnant fluid as described by Dütsch et al. (1998). The periodic oscillation is given by the harmonic function, x = -Asin (2πf_et) , in which A denotes the oscillation amplitude of the cylinder and f_e is the frequency of the oscillation. Keulegan-Carpenter and Reynolds numbers are defined as KC = U_m /f_ed and Re = U_m d/ν, where, Um is the maximum velocity of cylinder during oscillation, d is the cylinder diameter and ν is the kinematic viscosity. The computation is performed at KC = 5 and Re = 100 at which the experimental and numerical results by Dütsch et al. (1998) are available. Figure 4 shows the velocity profiles along the transverse (y_a) axis at x_a = -0.6d for two different phase positions compared to the experimental results of Dütsch et al. (1998). Velocity profiles (u_a) agree relatively well with the experimental results.

Figure 4
The x velocity component in the transverse (y_a) direction at x_a = –0.6d for two different phase positions (φ = 2πft). (a) φ = 180º; (b) φ = 210º.

Also, pressure and vorticity isolines have been compared with the numerical study of Guilmineau and Queutey (2002) in Figs. 5 and 6. Good agreement is evident.

Figure 5
Vorticity isolines for two different phase positions (φ = 2πft). (a) Guilmineau and Queutey 2002, at φ = 0º; (b) Present method, at φ = 0º; (c) Guilmineau and Queutey 2002, at φ = 95º; (d) Present method, at φ = 95º.

Figure 6
Pressure isolines for two different phase positions (φ = 2πft). (a) Guilmineau and Queutey 2002, at φ = 0º; (b) Present method, at φ = 0º; (c) Guilmineau and Queutey 2002, at φ = 95º; (d) Present method, at φ = 95º.

THE FLAPPING AIRFOIL MODEL

The computational domain and the flapping airfoil are shown in Fig. 7. The chord of the airfoil C is 23 mm and the diameter of the leading edge semi-circle D is 5 mm. These dimensions are chosen to simulate the experimental study carried out by Godoy-Diana et al. (2008).

Figure 7
Computational domain in this study.

Three non-dimensional parameters, i.e. Reynolds number, Re = U_∞D/v, normalized amplitude, A_D = A/D, and the Strouhal number, St = fD/U_∞, are important in this study, where f is the oscillation frequency. The oscillatory motion of the airfoil is described by the following equation:

(25) $\mathit{\theta }={\mathit{\theta }}_0+{\mathit{\theta }}_A\sin (2\Pi \mathrm{ft})$

In Eq. 25, θ₀ is the initial angle of attack and θ_A is the maximum angle of attack associated with the amplitude A_D, as shown in Fig. 7. A 1,200 x 720 uniform grid is used as the background Cartesian mesh, and the ∆t depends on the oscillation frequency. The time step is chosen so that each cycle of the oscillation corresponds to 2,000 time steps. For example, when the Strouhal number is St = 0.22, the time step is ∆t = 1.14 x 10^-4 s. Relevant flow parameters are taken from the experiment of Godoy-Diana et al. (2008), i.e. Re = 255, U_∞ = 0.1 m/s, and ρ = 1 kg/m³. The viscosity is tuned in each numerical test so that Re = 255 is kept fixed in all cases.

A grid refinement study was carried out at A_D = 0.71 and St = 0.22. The results are shown in Table 2 for 4 grid sizes. In Table 2, C is the chord of the airfoil and CD is the drag coefficient computed as follows:

Table 2
Grid refinement study for drag coefficient.

Grid size	Grid resolution	CD
0.0625C	400 x 240	0.781
0.0416C	600 x 360	0.762
0.03125C	800 x 480	0.752
0.025C	1,000 x 600	0.748
0.02083C	1,200 x 720	0.745
0.0166C	1,500 x 900	0.745

(26) $C_D=\frac{\sum _{k=1}^N{\left(\overrightarrow{F}\left({\overrightarrow{X}}_k\right)\right)}_x}{\frac{1}{2}\mathit{\rho }{U}_{\infty }^{2}h}$

In Eq. 26, (_k))_x terms are the x-components of the forces at Lagrangian points. (

All results presented in this paper are obtained using the 1,200 x 720 grid. Figure 8 shows the grid convergence study. The results for the finest grid are considered accurate, and the L₂ and L_∞ norms of the error obtained on the coarser grids are calculated and shown in Fig. 8. The results demonstrate the second-order accuracy of the method. The norms L₂ and L_∞ are defined as follows:

Figure 8
Grid convergence study. (a) L2 norm; (b) X norm.

(27) $L_2=\frac{1}{n}{\left(\sum _{i=1}^n\left({\left|{\Phi }_{\mathrm{coarser}}-{\Phi }_{\mathrm{finest}}\right|}^2\right)i\right)}^{\frac{1}{2}}$

(28) $L_{\infty }=\max \left|{\Phi }_{\mathrm{coarser}}-{\Phi }_{\mathrm{finest}}\right|$

where:

n is the number of grid points; φ is the velocity component.

SIMULATION RESULTS

Before presenting the computational results, it is appropriate to briefly discuss some phenomena associated with flapping airfoils. As explained next, these phenomena are simulated and observed in this numerical study and have also been reported by Yu et al. (2012) and Ren et al. (2013).

An airfoil flapping in a fluid generates two main vortices near the leading edge. Defining the positive vortex as a counter clockwise vortex, as shown in Fig. 9, a negative vortex is generated at the upper surface of the foil and a positive vortex is generated at its lower surface. The strengths of these leading edge vortices and their subsequent positions behind the foil depend on the flapping amplitude and frequency for a given Re. When these vortices shed away from the airfoil, they alter the flow pattern in the wake region. Three different scenarios regarding the orientation of the vortex pair in the wake region are shown in Fig. 9. For low values of oscillation frequency/amplitude, the upper (negative) and lower (positive) vortices keep their initial relative positions with respect to the symmetry line of the airfoil. In this mode of flapping motion, the momentum of the flow in the wake region is reduced due to the adverse effect of vortex system on the main flow (Fig. 9a). At higher frequencies/amplitudes, the vortex pair rolls into the wake region in an aligned situation shown in Fig. 9b. At still higher frequency/amplitude, the leading edge vortices cross the symmetry line behind the foil. In other words, the negative vortex appears below the symmetry line and the positive one appears above that line in this case (Fig. 9c). The flow symmetry breaks down at more intensified flapping, and a diverted jet of energized flow is observed.

Figure 9
(a) Kármán vortex street (drag-generating vortex); (b) Aligned vortex; (c) Reverse Kármán vortex street (thrust-generating vortex). The black arrow shows the main flow direction and the gray one shows the flow direction generated by vortices.

The vortex pair shown in Fig. 9a, results in momentum loss, and the flow pattern is known as the Kármán (or Benard-von Kármán) vortex street. These vortices are clearly drag generating vortices. In contrast, the vortex pair shown in Fig. 9c enhances the wake flow momentum, and the pattern is known as the reverse Kármán vortex street. The RKVS is, therefore, a thrust-generating vortex street. When the two vortices are aligned, as shown in Fig. 9b, they have no favorable or adverse effect on the main flow, shown by the thick black arrow.

The asymmetric or deflected jet flow regime is a thrust and lift-producing regime and is here called the asymmetric reverse Kármán vortex street regime (ARKVS) for the sake of briefness. In other words, in addition to the thrust production, side force is also generated due to the asymmetric nature of the flow in the ARKVS regime.

If the flapping airfoil is used as a propulsive device, it is obviously necessary to realize how and when the flow regime changes. For a specified Re, a performance map can be developed for both analysis and design purposes.

In this paper, two groups of numerical studies are carried out and reported to quantitatively investigate the phenomena just described. First, the St is fixed, and the effects of oscillation amplitude are considered. Then, the oscillation amplitude is fixed, and the effects of the frequency variation or St are studied. Finally, a map which shows the effects of the variation of both frequency and amplitude is presented.

FLAPPING AT FIXED FREQUENCY (St = 0.22)

Vortex structures and calculated average fluid velocities at grid points are shown in Fig. 10 for St = 0.22, as well as A_D = 0.36; 0.71; 1.07; 1.77 and 2.14. The calculated velocities are average values in an oscillation cycle, and the contour of the velocity field is also shown in each case for further clarification.

Figure 10
(a) Vortex shedding and (b) Time-averaged velocity at St = 0.22 and A_D = 0.36; 0.71; 1.07; 1.77 and 2.14.

As Fig. 10 shows two main leading edge vortices formed and shed away from the airfoil as a result of the flapping motion. It is clearly seen that boundary layer separation from the top and bottom sides of the flapping airfoil alters the pattern of the vortex street. At A_D = 0.36, momentum loss in the wake region is observed in the corresponding velocity-time averaged contour diagram. This is a drag generating KVS flow regime. By increasing the amplitude, negative and positive vortices are positioned along the symmetry line behind the airfoil at A_D = 0.71. In this aligned vortex regime, some loss in fluid momentum exists due to viscous effects. However, the vortex system does not play a significant role in momentum transfer downstream. Based on the numerical results, transition from the KVS to RKVS regime occurs at A_D ∼ 0.72. Further increase in the flapping amplitude, represented by A_D = 1.07 and A_D = 1.77, results in the intensification of the symmetric jet flow and the propulsive force. This intensification of the momentum transfer is clearly seen in the corresponding velocity-time averaged contour diagrams. At a higher flapping amplitude, i.e. A_D = 2.14, the symmetry in the flow field breaks down, and the asymmetric flow regime is observed. This corresponds to the ARKVS regime in which reactive side forces are generated.

The physics behind the asymmetric flow regime is not understood quite well, but the phenomenon may possibly be attributed to the intensive interactions between leading and trailing edge vortices (Yu et al. 2012). It is known that in the ARKVS regime a strong dipole vortex is formed on one side of the symmetry line, and a weaker single vortex is shed on the other side. Dipole vortex is made up of two vortices with different signs that are shed in each pitching cycle. In this vortex structure, corresponding to A_D = 2.14, the jet-like flow is bent towards the side of the strong dipole vortex as shown in Fig. 10.

The schematics in Fig. 11 help to better understand why the KVS flow regime changes to the RKVS regime. At low oscillation amplitudes, the induced momentum generated by airfoil pitching motion is not strong enough to push the vortex generated at one side of the foil to move to the other side. By increasing the oscillation amplitude, the push becomes strong enough to switch the positions of negative and positive vortices and to convert the KVS regime to the thrust-producing RKVS regime.

Figure 11
(a) The induced momentum pushes the negative vortex at the upper surface downwards; (b) The induced momentum pushes the positive vortex at the lower surface upwards.

To further investigate and quantify the intensified jet flow behind the foil, time averaged horizontal velocities along a-b and c-d lines, shown in Fig. 7, are calculated and shown in Figs. 12 and 13, respectively. It is seen that the velocity behind the airfoil at A_D = 0.36 is decreased due to the drag generating KVS regime. At A_D = 0.71, no net momentum is transferred to the main flow by the aligned vortex system. However, viscous effects reduce the velocity behind the airfoil slightly. At A_D = 1.07 and 1.77, velocity and momentum behind the airfoil increase due to the thrust-generating RKVS regime. Figure 9 makes it clear that there is a jet-like flow at the middle of c-d line. This corresponds to the transition from KVS to RKVS.

Figure 12
Time averaged horizontal velocity along a-b line, shown in Fig. 7, at St = 0.22 and 4 different amplitudes. A non-dimensional coordinate along a-b line is employed.

Figure 13
Time averaged horizontal velocity along c-d line, shown in Fig. 7, at St = 0.22 and 4 different amplitudes. A non-dimensional coordinate along c-d line is employed.

Simulation results also provide detailed information regarding the transition from symmetric RKVS to ARKVS flow regime. By considering the results shown in Fig. 14, it is concluded that the transition occurs at normalized amplitude very close to A_D = 2 for this particular flapping airfoil.

Figure 14
Asymmetric KVS at (a) A_D = 2, (b) 2.04 and (c) 2.08.

At A_D = 2, the dipole vortex is not recognizable but the deflection of the vortex sheet has just been initiated. At higher A_D values, the jet deflection and dipole vortex structure are clearly seen.

The calculated transient lift coefficients corresponding to A_D = 0.71; 1.77 and 2.14 are shown in Fig. 15. Here, the lift coefficient is calculated as follows:

Figure 15
Transient lift coefficient for St = 0.22 and A_D = 0.71, 1.77 and 2.14.

(29) $C_L=\frac{\sum _{k=1}^{N_b}{\left(\overrightarrow{F}\left({\overrightarrow{X}}_k\right)\right)}_y}{\frac{1}{2}\mathit{\rho }{U}_{\infty }^{2}h}$

where:

(_k))_y values are the y-components of the force at boundary (Lagrangian) points. Note that C_L curves corresponding to A_D = 0.71 and 1.77 are symmetric with respect to the C_L = 0 line, meaning that no side (lift) force is generated at these amplitudes. In contrast, the C_L curve corresponding to A_D = 2.14 is not symmetric with respect to the C_L = 0 line, which means that there is a net side force (lift) in this case. The average lift coefficients at various flapping amplitudes and St = 0.22 are shown in Fig. 16. It is clearly seen that lift generation starts close to A_D = 2. As shown in this figure, when the deflected vortex is generated behind the airfoil, a side force (lift) is produced. Therefore, the side force is an indication of the formation of a deflected vortex behind the airfoil. (

Figure 16
The effect of oscillation amplitude on average lift coefficient for St = 0.22.

To understand what exactly happens near the transition point, which results in transition from symmetric to deflected vortex, the vortex patterns are investigated at 4 time steps, t = 2T; t = 4T; t = 6T and t = 8T, being T the period of oscillation. As shown in Fig. 17, the first dipole vortex separates from the trailing edge at t = 2T. This dipole vortex is strong enough to bend the flow direction. After that, second and third dipole vortices separate from the trailing edge and follow the first dipole vortex. These dipole vortices are then stretched along the vortex street and result in symmetry breaking. It is clear that the first dipole vortex is very important in the sense that it determines the deflected vortex direction. Numerical analysis shows that, when initial phase angle alters, the deflected vortex direction alters too. This is clearly seen in Fig. 18.

Figure 17
Asymmetric KVS at (a) t = 2T, (b) 4T, (c) 6T and (d) 8T.

Figure 18
Shedding vortices at 2 initial phase angles with 90º difference. (a) Initial phase = 0º; (b) Initial phase = 90º.

Table 3 summarizes the vortex types and flow field effects for various AD values at St = 0.22.

Table 3
Vortex types and flow field effects for various A_D values at St = 0.22.

No.	Type of vortex	Effect on the flow	Range of oscillation amplitude
1	KVS (drag generating vortex)	Decreasing momentum behind the airfoil/ increasing the drag force	AD < 0.70
2	Aligned vortex street (neutral type vortex)	Neither decreasing nor increasing momentum behind the airfoil	0.70 < AD < 0.72
3	RKVS (thrust generating vortex)	Increasing momentum behind the airfoil/ decreasing the drag force and then generating a thrust force	0.72 < AD < 2
4	Deflected vortex (thrust/lift generating vortex)	Simultaneous thrust and lift force generation	2 < AD

FLAPPING AT FIXED AMPLITUDE (A _D = 0.71)

The effect of St on the vortex structure and thrust performance is now studied for fixed oscillation amplitude A_D = 0.71 and different values of St: 0.1; 0.22; 0.3; 0.4 and 0.5. Vortex shedding and average velocity contours in this case are shown in Fig. 19.

Figure 19
(a) Vortex shedding; (b) Time-averaged velocity at A_D = 0.71, and St = 0.1; 0.22; 0.3; 0.4 and 0.5.

Once again, to further explore the flow features, average horizontal velocities on a-b and c-d lines are calculated and shown in Figs. 20 and 21, respectively. The same trends explained before regarding Figs. 12 and 13 are observed here again.

Figure 20
Average horizontal velocities at grid points along a-b line at A_D = 0.71 and St = 0.1; 0.22; 0.3 and 0.4. A non-dimensional coordinate along a-b line is employed.

Figure 21
Average horizontal velocities at grid points along c-d line at A_D = 0.71 and St = 0.1; 0.22; 0.3 and 0.4. A non-dimensional coordinate along c-d line is employed.

Figure 22 shows the time history of lift coefficient for three different Strouhal numbers, i.e. St = 0.22; 0.3 and 0.5, and Fig. 23 shows the effect of the St on average lift coefficient for A_D = 0.71. The same trends explained before regarding Figs. 15 and 16 are observed here again. It is clear that lift is generated for St > 0.48 in this case.

Figure 22
Transient lift coefficient for A_D = 0.71 and St = 0.22; 0.3 and 0.5.

Figure 23
The effect of St on average lift coefficient for A_D = 0.71.

Table 4 summarizes the vortex types and flow field effects for various St values at A_D = 0.71.

Table 4
Vortex types and flow field effects for various St values at A_D = 0.22.

No.	Type of vortex	Effect on the flow	Range of Strouhal number
1	KVS (drag generating vortex)	Decreasing momentum behind the airfoil/ increasing the drag force	St < 0.22
2	Aligned vortex street (neutral type vortex)	Neither decreasing nor increasing momentum behind the airfoil	0.21 < St < 0.23
3	RKVS (thrust generating vortex)	Increasing momentum behind the airfoil/ decreasing the drag force and then generating a thrust force	0.23 < St < 0.48
4	Deflected vortex (thrust/lift generating vortex)	Simultaneous thrust and lift force generation	0.48 < St

PERFORMANCE MAP FOR A RANGE OF ST NORMALIZED AMPLITUDE

The effects of both amplitude and frequency on the drag/thrust coefficient are shown in Fig. 24 and compared to the results reported by He et al. (2012). The drag coefficient for the corresponding motionless airfoil is C_D0 = 0.86. As mentioned before, in the RKVS regime, the vortex system transfers momentum to the flow and generates thrust force. The situation C_D = 0 in Fig. 24 corresponds to a situation in which the RKVS regime generates just enough thrust to overcome the drag. The reduction of drag coefficient for high values of amplitude/frequency corresponds to more thrust generation in the RKVS regime.

Figure 24
Drag coefficients in terms of the normalized amplitude at St = 0.1; 0.22 and 0.3.

If the flapping airfoil is used for the purpose of generating propulsion or lift, it is obviously necessary to realize how and when the flow regime changes. For a fixed Re, a performance map can be developed for both analysis and design purposes. Such a map, also called the phase space, is experimentally obtained by Godoy-Diana et al. (2008) for the flapping airfoil just described in "The Flapping Airfoil Model" section. Here, the map is reproduced and shown in Fig. 25 based on the numerical results obtained in this study. Different regions in Fig. 25 represent different types of vortex systems and flow regimes described previously.

Figure 25
A flapping foil performance map. Zone 1: Kármán vortex street; zone 2: aligned vortex street; zone 3: reverse Kármán vortex street; zone 4: asymmetric reverse Kármán vortex street.

The KVS regime, which is characterized by zone number 1 in Fig. 25, is observed when at least one of A_D or St values is sufficiently low. The lower-left zone of the map corresponds to the KVS. At moderately low A_D and St values, characterized by zone number 2, the aligned vortex system is developed. This corresponds to a flow regime between KVS and RKVS. Zone 3 in Fig. 25 represents the RKVS regime. As shown in Fig. 25, the RKVS is formed at the middle of the map. Note that the zero drag (C_D = 0) line lies in this region of the map as well. At the left-hand side of the C_D = 0 line in zone 3, the generated thrust is still less than the drag. Propulsive force is obtained at oscillations corresponding to the right-hand side of the C_D = 0 line in zone 3. Finally, the thrust and lift generating regime, i.e. the ARKVS regime, is represented by zone 4 in Fig. 25. This regime is obtained when both frequency and amplitude are sufficiently high.

CONCLUSIONS

In this paper, an IBM is employed to simulate the flow around a flapping airfoil. The boundary conditions are accurately implemented by an iterative procedure applied at each time step. Vortex and wake patterns as well as lift and drag coefficients are studied at different oscillation amplitudes and frequencies. Four distinguished flow regimes, controlled by the vortex system generated due to the flapping motion, are observed. These flow regimes are the KVS, the Aligned Vortex Street (AVS), the RKVS and finally the ARKVS. The KVS regime generates drag; the AVS is a neutral regime with no net momentum transfer to the main flow; the RKVS is a thrust producing regime; and the ARKVS is a thrust and lift generating flow regime. For the particular symmetric airfoil used in this study, two groups of flapping scenarios are studied at a fixed Reynolds number (Re = 255). First, the frequency is fixed at St = 0.22, and the amplitude is changed. Then the amplitude is fixed at A_D = 0.71, and the oscillation frequency is altered. For the fixed frequency case, St = 0.22, the transition from KVS to RKVS occurs at A_D = 0.72, and the ARKVS is developed around A_D = 2. For the fixed amplitude oscillation, A_D = 0.71, the transition from KVS to RKVS occurs at St = 0.23, and the ARKVS is developed around St = 0.48. Vorticity and velocity contours are presented, which helps to understand the physics behind the flow. The effects of the dipole vortex system in deflecting the wake flow are also discussed. The computational results compared to the available experimental data show that the iterative-direct forcing immersed boundary method is a reliable tool for studying complex transient flows such as the flow around a flapping airfoil.

REFERENCES

• Anderson JM, Streitlien K, Barrett DS, Triantafyllou MS (1998) Oscillating foils of high propulsive efficiency. J Fluid Mech 360:41-72. doi: 10.1017/S0022112097008392
  » https://doi.org/10.1017/S0022112097008392
• Ashraf MA (2010) Numerical simulation of the flow over flapping airfoils in propulsion and power extraction regimes (PhD thesis). Australia: School of Engineering and Information Technology.
• Bai X, Avital EJ, Munjiza A, Williams JJR (2014) Numerical simulation of a marine current turbine in free surface flow. Renew Energ 63:715-723. doi: 10.1016/j.renene.2013.09.042
  » https://doi.org/10.1016/j.renene.2013.09.042
• Balaras E (2004) Modeling complex boundaries using an external force field on fixed cartesian grids in large eddy simulations. Comput Fluids 33(3):375-404. doi: 10.1016/S0045-7930(03)00058-6
  » https://doi.org/10.1016/S0045-7930(03)00058-6
• Bohl DG, Koochesfahani MM (2009) MTV measurements of the vertical field in the wake of an airfoil oscillating at high reduced frequency. J Fluid Mech 620:63-88. doi: http://dx.doi.org/10.1017/S0022112008004734
  » http://dx.doi.org/10.1017/S0022112008004734
• Cheny Y, Botella O (2010) The LS-STAG method: a new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties. J Comput Phys 229(4):1043-1076. doi: 10.1016/j.jcp.2009.10.007
  » https://doi.org/10.1016/j.jcp.2009.10.007
• Chung MH (2006) Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape. Comput Fluids 35(6):607-623. doi: 10.1016/j.compfluid.2005.04.005
  » https://doi.org/10.1016/j.compfluid.2005.04.005
• Davis WA (2007) Nano air vehicle: a technology forecast. Technical Report. Montgomery (AL): Blue Horizons Paper, Center for Strategy and Technology, Air War College.
• Dütsch H, Durst F, Becker S, Lienhart H (1998) Low-Reynoldsnumber flow around an oscillating circular cylinder at low Keulegan-Carpenter numbers. J Fluid Mech 360:249-264. doi: 10.1017/S002211209800860X
  » https://doi.org/10.1017/S002211209800860X
• Godoy-Diana R, Aider JL, Wesfreid JE (2008) Transitions in the wake of a flapping foil. Phys Rev E 77:207-221. doi: 10.1103/PhysRevE.77.016308
  » https://doi.org/10.1103/PhysRevE.77.016308
• Guilmineau E, Queutey P (2002) A numerical simulation of vortex shedding from an oscillating circular cylinder. J Fluid Struct 16(6):773-794. doi: 10.1006/jfls.2002.0449
  » https://doi.org/10.1006/jfls.2002.0449
• He GY, Wang Q, Zhang X, Zhang SG (2012) Numerical analysis on transitions and symmetry-breaking in the wake of a flapping foil. Acta Mech Sinica 28(6):1551-1556. doi: 10.1007/s10409-012-0158-8
  » https://doi.org/10.1007/s10409-012-0158-8
• Jones KD, Dohring CM, Platzer MF (1998) Experimental and computational investigation of the Knoller-Betz effect. AIAA Journal 36(7):1240-1246. doi: 10.2514/2.505
  » https://doi.org/10.2514/2.505
• Khalid MSU, Akhtar I, Durrani NI (2014) Analysis of Strouhal number based equivalence of pitching and plunging airfoils and wake deflection. J Aerospace Eng 229(8). doi: 10.1177/0954410014551847
  » https://doi.org/10.1177/0954410014551847
• Kim J, Moin P (1985) Application of a fractional-step method to incompressible Navier-Stokes equations. J Comput Phys 59(2):308-323. doi: 10.1016/0021-9991(85)90148-2
  » https://doi.org/10.1016/0021-9991(85)90148-2
• Lewin GC, Haj-Hariri H (2003) Modeling thrust generation of a two dimensional heaving airfoil in a viscous flow. J Fluid Mech 492:339-362. doi: 10.1017/S0022112003005743
  » https://doi.org/10.1017/S0022112003005743
• Lima e Silva ALF, Silveira-Neto A, Damasceno JJR (2003) Numerical simulation of two-dimensional flows over a circular cylinder using the immersed boundary method. J Comput Phys 189(2):351-370. doi: 10.1016/S0021-9991(03)00214-6
  » https://doi.org/10.1016/S0021-9991(03)00214-6
• Mittal R, Dong H, Bozkurttas M, Najjar FM, Vargas A, Von Loebbecke A (2008) A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J Comput Phys 227(10):4825-4852. doi: 10.1016/j.jcp.2008.01.028
  » https://doi.org/10.1016/j.jcp.2008.01.028
• Mohd-Yusof J (1997) Combined immersed-boundary/B-spline methods for simulations of flow in complex geometries. Center for Turbulence Research Annual Research Briefs 161(1):317-327.
• Park J, Kwon K, Choi H (1998) Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160. J Mech Sci Technol 12(6):1200-1205. doi: 10.1007/BF02942594
  » https://doi.org/10.1007/BF02942594
• Peskin C (2002) The immersed boundary method. Acta Numer 11:479-517. doi: 10.1017/S0962492902000077
  » https://doi.org/10.1017/S0962492902000077
• Ren W, Hu H, Liu H, Wu JC (2013) An experimental investigation on the asymmetric wake formation of an oscillating airfoil. Proceedings of the 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Paper. 0794; Grapevine, USA.
• Schnipper T, Andersen A, Bohr T (2009) Vortex wakes of a flapping foil. J Fluid Mech 633:411-423. doi: 10.1017/S0022112009007964
  » https://doi.org/10.1017/S0022112009007964
• Su SW, Lai MC, Lin CA (2007) An immersed boundary technique for simulating complex flows with rigid boundary. Comput Fluids 36(2):313-324. doi: 10.1016/j.compfluid.2005.09.004
  » https://doi.org/10.1016/j.compfluid.2005.09.004
• Taylor GK, Nudds RL, Thomas ALR (2003) Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency. Nature 425:707-711. doi: 10.1038/nature02000
  » https://doi.org/10.1038/nature02000
• Tseng YH, Ferziger JH (2003) A ghost-cell immersed boundary method for flow in complex geometry. J Comput Phys 192(2):593-623. doi: 10.1016/j.jcp.2003.07.024
  » https://doi.org/10.1016/j.jcp.2003.07.024
• Uhlmann M (2005) An immersed boundary method with direct forcing for the simulation of particulate flows. J Comput Phys 209(2):448-476. doi: 10.1016/j.jcp.2005.03.017
  » https://doi.org/10.1016/j.jcp.2005.03.017
• Von Kármán T, Burgers MJ (1934) Aerodynamic theory. Vol. 2. Berlin: Springer.
• Wang Z, Fan J, Luo K (2008) Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles. Int J Multiphas Flow 34(3):283-302. doi: 10.1016/j.ijmultiphaseflow.2007.10.004
  » https://doi.org/10.1016/j.ijmultiphaseflow.2007.10.004
• Wang ZJ (2000) Vortex shedding and frequency selection in flapping flight. J Fluid Mech 410:323-341. doi: 10.1017/S0022112099008071
  » https://doi.org/10.1017/S0022112099008071
• Wei Z, Zheng ZC (2014) Mechanisms of wake deflection angle change behind a heaving airfoil. J Fluid Struct 481-13. doi: 10.1016/j.jfluidstructs.2014.02.010
  » https://doi.org/10.1016/j.jfluidstructs.2014.02.010
• Weis-Fogh T (1973) Quick estimates of flight fitness in hovering animals, including novel mechanism for lift production. J Exp Biol 59:169-230.
• Young J, Lai JCS (2007) Mechanisms influencing the efficiency of oscillating airfoil propulsion. AIAA Journal 45(7):1695-1702.
• Yu ML, Hu H, Wang ZJ (2012) Experimental and numerical investigations on the asymmetric wake vortex structures around an oscillating airfoil. Proceedings of the 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Paper. 0299; Nashville, Tennessee.
• Zanotti A, Melone S, Nilifard R, D'Andrea A (2014) Experimental-numerical investigation of a pitching airfoil in deep dynamic stall. J Aerospace Eng 228(4):557-566. doi: 10.1177/0954410013475954
  » https://doi.org/10.1177/0954410013475954

Publication Dates

• Publication in this collection
  Jul-Sep 2015

History

• Received
  03 Mar 2015
• Accepted
  05 Dec 2015