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 (AD) numbers are presented. It is shown that, at the fixed pitching amplitude, AD = 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 AD > 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.


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, Numerical Simulation of the Wake Structure and Thrust/Lift Generation of a Pitching Airfoil at Low Reynolds Number Via an Immersed Boundary Method 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.

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.
The force terms on Lagrangian points, represented by F , 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, f, are related to the boundary force as follows: where: x = (x, y) are the Cartesian coordinates of an Eulerian point; 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 = (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 u.The velocities at boundary Lagrangian points, U, are calculated afterwards as follows: 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: It is important to note that u 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: where: u is the velocity vector on the Eulerian points; P is the pressure; ρ is the fluid density; ν is the kinematic viscosity, defined as ν = μ / ρ , where μ is the fluid dynamic viscosity; V and V 2 are the regular centered difference approximations for the gradient and Laplace operators, respectively.The source term f in Eq. 3 is a fraction of the relevant external boundary force F that acts at Eulerian grid point.The force term F is imposed at the Lagrangian point to enforce the desired velocity U Γ at the immersed boundary.This is done to assure that the velocities at Lagrangian points satisfy the no-slip boundary condition.
An estimate of the velocity can be obtained by the following equation: Rewriting Eq. 7 at a Lagrangian point, it results in: where: U n+1 is the velocity at the Lagrangian point at time level n + 1; U is an intermediate Lagrangian point velocity calculated as follows: The velocity term U n+1 in Eq. 8 can now be replaced by the known velocity at the boundary (U Γ ).This results in the following expression: Force terms at Eulerian points can now be calculated using the discrete delta function: 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, X 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: 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: where: (u, û, u*)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: 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 u (intermediate velocity at Eulerian points) is calculated by Eq. 14.
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: Vector quantities, u, u* and u* , are the intermediate velocities between the time levels n and n + 1.
The computational algorithm can now be summarized as follows: 1. Obtain the intermediate velocity at grid points, u (x), by solving Eq. 14. Start the inner loop: 2. Obtain the intermediate velocity at grid points, û (x), by solving Eq. 15. 3. Calculate the intermediate velocities at Lagrangian points, U (X k ) , using Eq.16. 4. Obtain the force terms at Lagrangian points, F n,m (X k ), by Eq. 17. 5. Distribute the Lagrangian forces among the Eulerian points using Eq.18.

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  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.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.
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.

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 In Eq. 26, (F (X 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 2 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 2 and L ∞ are defined as follows: where: n is the number of grid points; ϕ is the velocity component.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.
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.
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.
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 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.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.clear that there is a jet-like flow at the middle of c-d line.This corresponds to the transition from KVS to RKVS.
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.
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.where: (F (X 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.
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 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.

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.
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

Figure 1 .
Figure 1.(a) Lagrangian (X) and Eulerian (x) points; (b) zone of influence for the force term at Lagrangian point k (Lima e Silva et al. 2003).
the Wake Structure and Thrust/Lift Generation of a Pitching Airfoil at Low Reynolds Number Via an Immersed Boundary Method (20) ˜˜ˆ˜ Numerical Simulation of the Wake Structure and Thrust/Lift Generation of a Pitching Airfoil at Low Reynolds Number Via an Immersed Boundary Method 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.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 e t) , 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 e d and Re = U m d/ν, where, U m 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.6dfor 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 2 .Figure 3 .
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.

Figure 7 .
Figure 7. Computational domain in this study.
simulate the experimental study carried out by Godoy-Diana et al. (2008).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: Numerical Simulation of the Wake Structure and Thrust/Lift Generation of a Pitching Airfoil at Low Reynolds Number Via an Immersed Boundary MethodIn Eq. 25, θ 0 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 ofGodoy-Diana et al. (2008), i.e.Re = 255, U ∞ = 0.1 m/s, and ρ = 1 kg/m 3 .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

Figure 9 .
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 Wake Structure and Thrust/Lift Generation of a Pitching Airfoil at Low Reynolds Number Via an Immersed Boundary Method 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.

Figure 11 .
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.

Figure 12 .
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 .
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.

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

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.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.
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.

Figure 21 .
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 20 .
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 24 .
Figure 24.Drag coefficients in terms of the normalized amplitude at St = 0.1; 0.22 and 0.3.

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

Table 1 .
Comparison between the drag coefficients computed by the present method and two reported results.
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.

Table 3 .
Vortex types and flow field effects for various A D values at St = 0.22.
Table 3 summarizes the vortex types and flow field effects for various AD values at St = 0.22.
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.