Vortex Lattice Simulations of Attached and Separated Flows around Flapping Wings

Flapping flight is an increasingly popular area of research, with applications to micro-unmanned air vehicles and animal flight biomechanics. Fast, but accurate methods for predicting the aerodynamic loads acting on flapping wings are of interest for designing such aircraft and optimizing thrust production. In this work, the unsteady vortex lattice method is used in conjunction with three load estimation techniques in order to predict the aerodynamic lift and drag time histories produced by flapping rectangular wings. The load estimation approaches are the Katz, Joukowski and simplified Leishman–Beddoes techniques. The simulations’ predictions are compared to experimental measurements from wind tunnel tests of a flapping and pitching wing. Three types of kinematics are investigated, pitch-leading, pure flapping and pitch lagging. It is found that pitch-leading tests can be simulated quite accurately using either the Katz or Joukowski approaches as no measurable flow separation occurs. For the pure flapping tests, the Katz and Joukowski techniques are accurate as long as the static pitch angle is greater than zero. For zero or negative static pitch angles, these methods underestimate the amplitude of the drag. The Leishman–Beddoes approach yields better drag amplitudes, but can introduce a constant negative drag offset. Finally, for the pitch-lagging tests the Leishman–Beddoes technique is again more representative of the experimental results, as long as flow separation is not too extensive. Considering the complexity of the phenomena involved, in the vast majority of cases, the lift time history is predicted with reasonable accuracy. The drag (or thrust) time history is more challenging.

. Flapping mechanism the entire assembly and to limit the inertial loads caused by the motion.    The loads measured by the aerodynamic balance included both aerodynamic and inertial 105 contributions. In order to separate the two, wind-off flapping and pitching tests were carried out for 106 each tested configuration. The force measurements obtained during these tests included only inertial 107 contributions and were to be subtracted from the load values measured during the wind-on tests, 108 which also included aerodynamic loads. The wind-off tests were carried out after replacing the wings 109 with metal bars that had the same inertial characteristics in flap and generated negligible aerodynamic 110 loads. Extension springs were also installed on top of the flapping arms in order to mimic the effect 111 of the lift acting on the wings at wind-on conditions. Despite these precautions, the flap and pitch 112 angle signals at wind-off and wind-on conditions were not identical, as demonstrated in figure 5 113 for a particular test case. All measured signals had to be cycle averaged before the subtraction of 114 the wind-off loads from the wind-on loads could be performed. Figure 5 shows the result of the between the two pitch response. Nevertheless, it was assumed that pitching has a very small effect 118 on the inertial loads compared to flapping and therefore the subtraction could be carried out. The 119 cycle averaging procedure was also applied simultaneously to the measured lift and drag signals,  Finally, the parasite drag due to skin friction and interference from the wing supports was also 125 removed from the drag measurements, since the vortex lattice method cannot model such effects.
Static wind-on tests were carried out with the wings set to zero flap and pitch angles and the measured simulations are carried out with a single wake shed from the trailing edge.

142
The two vorticity surfaces (bound and trailing edge wake) are exemplified in figure 6. Note that 143 x, y and z are global coordinates but x is aligned with the chord and y is aligned with the span when 144 γ = θ = 0 • . Once a vortex ring is shed into the wake, it is allowed to travel at the local flow velocity, 145 so that the wakes are force-free.

146
According to the vortex model by Vatistas et al. [37], the velocity induced at a general field point by a straight line vortex segment, whose ends with respect to the field point are given by vectors r 1 and r 2 and whose strength is Γ, is given by where K is defined as K = 1 if the vortex segment belongs to a bound vortex ring and if the segment belongs to a wake vortex ring. In this latest expression where ν is the kinematic viscosity of air, δ = 1 + 2 × 10 −4 √ Γ 2 /ν and t v is the time elapsed since the 147 vortex segment was shed into the wake.

148
The impermeability boundary condition is imposed by calculating the velocities induced by all the vortex rings normal to the control points of the geometric panels of the wing from equation 1 and setting the total normal velocity to zero. The boundary condition becomes where Γ b and Γ TE are vectors containing the vortex strengths of the bound and trailing edge wake gives rise to a starting vortex. This means that the first flapping cycle is not representative of the fully 159 developed flapping flow and at least two cycles must be simulated.

160
Once the strengths of the bound vortex rings have been calculated, the forces acting on the wing can be calculated using either the Joukowski approach or the Katz approximation [34,35]. For the former, the force generated by the i, jth vortex ring is given by where i is the chordwise panel index, j is the spanwise panel index and F st m is the steady contribution of the mth vortex segment of the ring, given by Here, U m is the local flow velocity evaluated at the centre of the segment, l m is the vortex segment vector and Γ m its vortex strength. The latter is generally equal to Γ ij for all the vortex segments of each vortex ring except for the trailing edge spanwise segments (i.e. the bound segments that are in contact with the wake segments), where Γ m = 0. Finally, F unst ij in equation 4 is the unsteady force contribution of the i, jth vortex ring, given by where A ij is the area of the corresponding wing panel and n ij a unit vector normal to the panel. The

161
derivative ∂Γ ij /∂t is calculated using a first order backwards difference formula.

162
The Katz approximation requires the calculation of three different velocities on the ijth panel:

163
• the velocity due to the wing's motion U m ij ,

164
• the velocity induced by the vorticity in the wake U w ij and

165
• the velocity induced by only the bound chordwise vortex segments, U bc ij 166 Then the lift acting on the i, jth panel is given by where τ c ij , τ s ij are chordwise and spanwise unit vectors tangential to the panel while ∆c ij and ∆b ij are the spanwise and chordwise lengths of the panel. The panel's angle of attack due to the shape and motion of the wing is The drag acting on the panel is given by |U m ij | 2 and P ij n ij is an orthogonal projection of the normal to the panel in the direction of the flow. The total aerodynamic force induced by the i, jth panel is then given by Note that, for the Katz method, L ij P ij n ij is always perpendicular to the local airspeed due to the 167 motion of the wing while D ij U m ij /|U m ij | is always parallel to it. Consider a static, symmetric and 168 rectangular wing at an angle of attack α to a horizontal free stream. The contribution L ij P ij n ij to the 169 total drag will be zero for any value of α, as the lift is always perpendicular to the free stream. The   183 An alternative to the Katz and Joukowski aerodynamic load estimation techniques is to use the 2D unsteady drag model presented by Leishman [38]. For a 2D airfoil oscillating in pitch in an airflow, the unsteady drag coefficient can be written as

Estimation of separated flow aerodynamic loads using the Leishman-Beddoes model
where α(t) is the instantaneous pitch angle, c n the normal force coefficient and c c the chordwise force coefficient. The latter is also known as leading edge suction as it contributes thrust. For attached flow, the normal force coefficient can be estimated from where α E is an unsteady effective angle of attack and the lift curve slope is assumed to be equal to 2π. The factor η usually takes values between 0.95 and 0.97 and represents the fact that airfoils in viscous flows cannot attain the full amount of leading edge suction. Leishman and Crouse [39] select the effective angle of attack as where c n is the unsteady normal force coefficient caused by circulatory contributions only.

184
In the context of a VLM simulation, the unsteady normal force coefficient can be obtained from either the Katz or Joukowski methods. For example, from equation 5 and for the jth spanwise section This force can be used to obtain an estimate of the local unsteady angle of attack α Ej This expression is an overestimation, as c n j contains both circulatory and non-circulatory 185 contributions. This point will be discussed later on.

186
The Leishman-Beddoes model can also be used to estimate the amount of local flow separation and the ensuing changes to the aerodynamic loads. The chordwise position of the separation point, f s , can be estimated from Kirchoff theory. For the jth spanwise section, the separation point lies at where α 1 is the value of the static angle of attack when f s = 0.7, i.e. the separation point lies at 187 0.7c and S 1 , S 2 are constants to be determined for each airfoil and Reynolds number. Furthermore, 188 α * Ej = (c n j − c n 0 )/2/pi, were c n 0 is the normal coefficient offset required to make the airfoil's lift curve 189 symmetrical (see next paragraph). Consequently, the normal and chordwise force coefficients at the 190 jth spanwise station are given by Finally, the total lift and drag can be approximated from The aerodynamic load coefficients of equations 12 and 13 were calculated at each time instance

209
As mentioned earlier, equation 9 is an overestimation, since the Katz approach does not differentiate between circulatory and non-circulatory aerodynamic loads. A simple and approximate solution to this problem is adopted in the present work. Equation 10 is re-written as and the correction factor in both equations 14 and 11 is set to η = 0.75. This modification implies 210 that around 20% of the inviscid aerodynamic normal force over the entire span is non-circulatory in       histories and concern cases where neither method is very successful (i.e. pure flapping with negative 365 θ 0 and pitch lagging). The Leishman-Beddoes approach can provide better drag predictions than the 366 other two methods for specific cases but can also feature a constant negative offset in its time histories, 367 that sometimes can be significant.