Numerical Study of Reduced Frequency Effect on Longitudinal Stability Derivatives of Airfoil under Pitching and Plunging Oscillations

J. Aerosp. Technol. Manag., São José dos Campos, Vol.8, No 3, pp.272-280, Jul.-Sep., 2016 ABSTRACT: In this study, incompressible, unsteady and turbulent flow over an airfoil with pitching and plunging oscillations is numerically studied in order to investigate the effect of reduced frequency on stability derivatives of oscillating airfoil. Linear k – ε model called Launder-Sharma and Rhie and Chow model are used for turbulence modeling and overcoming pressure checkerboard problem. This means that a co-located approach is proposed in this paper to study a moving grid problem, and the results demonstrate the high accuracy of the method. Control volume and Crank-Nicholson discretization method are also used for the numerical solution. It is shown that the longitudinal stability derivatives of plunging and pitching motions trend change intensively beyond the stall angle of attack while pitching rate has a completely opposite behavior. The results also show that increasing reduced frequency leads to stability reduction in plunging oscillation but it does not have significant effect on pitching oscillation case in pre-stall, stall and post-stall conditions. Grid convergence is examined to assess the accuracy of the numerical method that shows the high accuracy of it and this is a prominent achievement of the present study. The results of the proposed method in forces and moment show a good agreement with the experimental data.


INTRODUCTION
One of the most important problems in the aerodynamics theory is analyzing the stability of aircraft in order to be assured of its controllability and safety.The accurate prediction of moment and forces applied on aircraft plays a prominent role in determining its stability.Unsteady aerodynamics of oscillating airfoils has received attention in order to analyze the aerodynamic forces on aircrafts.Glauert (1935), in the 1930s, commenced studying oscillating airfoils.Following him, in 1940, Theodorsen and Garrick (1948) started much more accurate investigations on flutter problem of aircraft wings and paved the way for future studies on inviscid and incompressible flows over low-amplitude oscillating airfoils.Although plenty of numerical and experimental studies have been carried out to determine aerodynamic coefficients of static airfoils, more investigations on unsteady aerodynamics of oscillating airfoils are demanded.Among experimental studies, there are the papers by Panda and Zaman (1994), who investigated the flow field and estimation of lift from wake survey on a pitching airfoil.The Reynolds numbers were 44,000 and 22,000 based on chord length and the airfoil was NACA 0012 with sinusoidal oscillation.Tolouei et al. (2004) studied the flow around a pitching airfoil in a frequency range of 0.022 -0.066.They particularly investigated the pressure distribution of an oscillating airfoil.Sadeghi et al. (2010) also experimentally studied unsteady wake of a pitching airfoil and investigated pitching amplitude effect on the wake thickness.
Tuncer (1986), in the 1980s, started a numerical investigation on rapidly pitched airfoils.He studied unsteady fields over Numerical Study of Reduced Frequency Effect on Longitudinal Stability Derivatives of Airfoil under Pitching and Plunging Oscillations pitching airfoil and dynamic stall in his paper.In the 1990s, Tuncer (1997) used overset grid method to numerically study a plunging airfoil and compared the results with single grid solution.In that paper, he used Baldwin-Brath as a turbulence modeling and the Mach 0.3 flow was considered.Guilmineau and Queutey (1999) numerically studied dynamic stall on different airfoil sections.They used Baldwin-Brath and k -ω SST turbulence models.These studies are devoted to find a numerical approach to deal with oscillating airfoils.However, there is still a compelling desire to achieve highly accurate numerical method on this subject.
Additionally, plenty of experimental and numerical investigations have been carried out on stability derivatives and stability of aircrafts.Bak et al. (2000) experimentally investigated force coefficients and stability derivatives of NACA 63-415 airfoil in wind tunnel tests.Park et al. (2003) predicted stability coefficients using Euler equations in pitching motion.They used multi-grid method to calculate pitch and roll moments.Altun and İyigün (2004) measured the dynamic stability derivatives of a generic combat aircraft model in wind tunnel applying force oscillation on their model.Gopinath and Jameson (2005) used time-spectral method to study force and moment coefficients on 2-D and 3-D bodies.They used multi-grid method to solve Navier-Stokes equations for pitching airfoils and pitching wings.Schmidt and Newman (2010) also numerically investigated stability derivatives of a generic combat aircraft model.In recent studies Bhagwandin and Sahu (2014) numerically calculated stability derivatives for finned projectiles in subsonic and transonic regimes.Considering these papers, more investigations are required to study parameters that have significant effects on stability derivatives and oscillating motions such as oscillation amplitude or reduced frequency.
In the present research, a co-located numerical method is proposed to a moving grid problem.The purpose is to accurately estimate the stability derivatives of an oscillating airfoil and investigate the reduced frequency effect on the aerodynamics of oscillating airfoils.The proposed method is highly accurate; therefore, it can be a solution for the most common issue of numerical methods in unsteady aerodynamics, which is low accuracy.Moreover, the proposed method is applied to a 2-D problem and is able to be applied on more complex 3-D problems either.Numerical finite volume method with linear k -ε, Launder-Sharma turbulence model is used to study unsteady, incompressible and turbulent flow over a pitching and a plunging NACA 63-415 airfoil.After computing force and moment coefficient diagrams against time, longitudinal stability derivatives are obtained using curve fitting in both oscillations.This curve fitting is based on standard approaches that help the method to be able to model non-linear responses on flying vehicles.In another remarkable point of this numerical study, Rhie and Chow (1983) method is devoted to overcome pressure-velocity fluctuations (pressure checkerboard problem).After obtaining stability derivatives, the reduced frequency effects in pre-stall, stall and post-stall conditions on these derivatives are studied.Grid convergence index (GCI) and order of accuracy are also obtained to verify the numerical method and approve its accuracy.

GOVERNING EQUATIONS
Continuity and momentum equations in Reynolds-Averaged Navier-Stokes (RANS) equations for unsteady incompressible flow are as follows: where: t represents time; x is the space direction; U is the mean velocity; i and j are indexes, with i, j = 1, 2, 3; P is the pressure; U i is the mean velocity of the flow; ρ is density; ν is viscosity; u i u j is called Reynolds stress terms that are unknowns of the problem and have to be modeled.
Due to the oscillation motion, the source term S i is zero.Boussinesq (1987) proposed a method in which Reynolds stress could be modeled using an equivalent turbulent viscosity and is expressed as: (1) where: δ ij is the Kronecker delta; v t is the turbulent kinematic viscosity and is calculated as: where: C μ and f μ values are given in Table 1.

Coefficients
Launder-Sharma model C ε1 1.44 In order to model Reynolds stress terms, linear k -ε turbulence model called Launder-Sharma is used.Turbulence kinetic energy (k) and dissipation rate (ε) equations are: where: V and s are integral volume and surface, respectively; dn j is the differential component of surface normal vector.
Note that, for a moving mesh problem, these equations are changed in some ways to model new control volume shapes.These corrections are made using Leibnitz rule.Let W j be the boundary velocity of a control volume.Using Leibnitz rule, integrated continuity and momentum equations are considered just like the previous section; the only difference is that it is written for a moving mesh and a control volume. (5) where: σ k , σ ε , C ε1 , C ε2 , f 1 , f 2 and E are given in Table 1.In these relations, P* is the production term in turbulence flow and is calculated as: where: S ij is the mean rate of stress tensor and is calculated this way:

NUMERICAL METHOD
In order to solve flow equations in this paper, control volume method is used.For stationary control volumes, integrated continuity and momentum equations for incompressible flow are: In the proposed numerical method, Crank-Nicholson method is used for discretization.To prevent pressure and velocity field fluctuation, Rhie and Chow method for incompressible flow is utilized.In this co-located method, velocity components on moving control volume boundaries are given as: where: u shows the velocity on volume boundaries while U indicates the velocity at the center of the control volume; e indexes are related to the east of the control volume and E indexes are related to the neighboring control volume centers; P indexes indicate the particular control volume that is being solved.
Figure 1 shows the scheme of a cell.
Transporting velocity is also modeled to prevent pressure checkerboard problem: Comparing different solutions, the most appropriate value for k΄ is obtained; k΄ = 1 is selected for Eq. 14 for an accurate and smooth solution.For boundary conditions in this numerical method, at inlet, velocity, kinetic turbulent energy, dissipation rate, and molecular viscosity are set, and, at outlet, pressure is applied.
where: α • is the variation rate of angle of attack and q is the pitch rate; C ma is the pitch damping that can be neglected because of its small value.

SIMULATION OF PITCHING AND PLUNGING OSCILLATIONS
In order to calculate longitudinal stability derivatives, unsteady aerodynamics are considered.Moreover, force and moment coefficients against time in different angles of attack have to be computed.Two types of oscillating motions are used in this paper: (1) pitching and (2) plunging.Thus, a harmonic oscillating motion with constant amplitude is applied on airfoil.Pitching motion with sinusoidal variation in angle of attack is described as: where: α m is the mean angle of attack; α A is the oscillation amplitude; ω is the circular frequency and is determined using the reduced frequency relation k c = ωl/2V , where l is the chord length and V is the freestream velocity.
In addition, for plunging oscillation, the equation of motion is: After applying pitching and plunging oscillations to the airfoil, force and drag coefficients are computed in time; then, using curve fitting, the following equation is assumed for the pitching moment: where: C mA is the amplitude; C m0 is the initial pitching moment; δ is the phase difference that is assumed by moment coefficient variations and IEEE-std-1050 standard for curve fitting.
Comparing Eq. 19 with Eq. 18, the following equation is obtained: Additionally, in order to use the previous equation for plunging airfoil, an equivalent angle of attack is defined:

GRID AND INITIAL CONDITIONS
Airflow Reynolds number is .The NACA 63-415 airfoil that is used in this study has chord length of 0.6 m and the inlet air velocity is 40 m•s -1 .Density is considered constant and is set at 1.225 kg•m -3 ; the pressure at the outlet is set at 1 bar.Since the reduced frequency effects on the stability longitudinal derivatives are investigated, constant mean angle of attack and oscillation amplitude are determined, and the values are 1.5° and 1.3°, respectively.The reduced frequency interval for this study is 0.09 to 0.16, C-type grid is used for the airfoil and the domain is shown in Fig. 2. Three different where: a is the mean oscillation amplitude.Using Taylor series for pitching moment and neglecting higher-than-2nd-order derivative terms, the pitching moment can be described based on reduced frequency for pitching airfoil (Ronch et al. 2011): grids with 11,200, 21,600, and 25,200 cells are used for grid refinement.The 21,600-cells grid is used for numerical solution.
Note that the y plus value for the grid is about 1.1 and the wall distance is about 10 −5 m.The grid convergence study has also been done and it is proposed in the next section.

GRID CONVERGENCE STUDY
For grid study convergence, Richardson extrapolation is used to calculate a higher-order estimation of the flow fields from a series of lower-order discrete values (f 1 , f 2 , …, f n ).For the case of grid refinement study, the value estimated from the Richardson extrapolation is the one that would result when the cell grid size is tended to zero.Roache (1994) also suggests a GCI to provide a consistent manner in reporting the results of grid convergence.To compute the GCI, 3 levels are recommended in order to accurately estimate the order of convergence and to check that the solutions are within the asymptotic range of convergence:

GEOMETRIC CONSERVATION LAW
Since a moving grid is considered in the present paper, in this section, a test case is introduced to demonstrate the correct performance of the algorithm on moving grids.The case is defined to prove that grid motions would not affect the solution.For this purpose, the fluid is initially stationary with the following parameters: P = 1 bar and ρ = 1.225 kg•m -3 .With these initial conditions, fluid flow equations are solved in the domain while the elements oscillate with the motion of displacement.After 3 s, with the time step of 0.001 s, the contours of dV (V − 0), variation of velocity magnitude after 3 s, that is physically expected to remain zero, are plotted in Fig. 3.As can be seen, they are almost zero.This demonstrates the correct implementation of Geometric Conservation Law in the solution.
Table 2 shows grid convergence study results for the proposed numerical method.Since moment coefficient is a significant parameter in this paper, it is used to study grid convergence and obtain the order of accuracy.The convergence is studied for moment coefficient in 2.1 o and 0.5 o angles of attack.Subscripts 3, 2 and 1 represent, respectively, 11,200-, 21,600-and 25,200-cells grids and ε ij = f i -f j .
As Table 2 shows, the order of accuracy for the numerical method is about 1.92 while the theoretical order of accuracy is 2. The difference is most likely due to grid stretching or turbulence modeling.There is a reduction in GCI value for the successive grid refinements (GCI 12 < GCI 23 ).The GCI for finer grid (GCI 12 ) is relatively low if compared to the coarser grid (GCI 23 ), indicating that the dependency of the numerical method on the cell size has been reduced.Additionally, as the GCIreduction from the coarser grid to the finer grid is relatively high, the grid independent solution can be said to have been nearly achieved.Further refinement of the grid will not give much change in the results.So, it can be said that the solution on the proper grid resolution is nearly grid-independent.
Figure 4 indicates the independency of solution from grids for lift and moment coefficients of a plunging airfoil.In Fig. 5, lift and drag coefficient diagrams in different angles of attack for 3 different time steps are shown to demonstrate the independency of results from time step.As Fig. 5 shows, 0.001 s is assumed for the time step of the solution.Note that the results are shown for a pitching airfoil in this figure.

NUMERICAL RESULTS
As Figs. 4 and 5 show, lift, drag and moment coefficient diagrams against angle of attack have hysteresis loop in a 0.25 to 2.75 interval.Due to sinusoidal oscillations in airfoil motion, force and moment coefficients have hysteresis loops in particular intervals of angles of attack for both pitching and plunging airfoils.Figure 6 shows forces and moment coefficients for mean angle of attack of 1.5°, reduced frequency of 0.09 and pitching amplitude of 1.3° for a pitching NACA 63-415 airfoil.Hysteresis loops are clearly observed in these diagrams.This figure also compares the results of this research with wind tunnel test results found in Bak et al. (2000) to validate the proposed numerical method.It is shown that the results of the present study have good accuracy in comparison with the existing experimental results.Numerical method errors are 3.2, 1.4 and 3.9 percent for lift, moment and drag coefficient,   In order to predict the time variation of longitudinal stability derivatives, force and moment diagrams against time have to be determined.Curve fitting these diagrams using Eq. 19, longitudinal stability derivatives can be calculated by Eq. 20.After this procedure, longitudinal stability derivative variations with different angles of attack can be obtained.Figure 6 shows the variation for longitudinal stability derivatives in different angles of attack.For this case, the mean angle of attack is 1.5° and the reduced frequency is 0.09 for both pitching and plunging airfoils.Longitudinal stability derivative variations for pitching (C mα •+ C mq ) and plunging airfoil (C mα •) with angle of attack are shown in Fig. 7. Numerical results show that the longitudinal stability derivatives of pitching and plunging airfoils have the same trend in different angles of attack, whilst the longitudinal stability of pitch rate (C mq ) has a completely different attitude in pre-and post-stall conditions.In other words, longitudinal stability derivatives for pitching and plunging airfoils values are increased in pre-stall condition and decreased in post-stall condition.Longitudinal stability derivatives of pitch rate have an opposite trend.
Figure 7 also shows that longitudinal stability derivatives of pitching and plunging airfoils have negative signs in all angles of attack, which describes the longitudinal stability in pitching and plunging motions.Longitudinal stability derivatives of pitch rate have negative sign after stall angle of attack and become more stable.In addition to that, in high angles of attack, longitudinal stability derivative values of plunging motion are positive, which means that this motion becomes unstable in post-stall condition.Longitudinal derivative values of pitching motion are still negative in this area, which indicates its stability.
Since reduced frequency is one of the main elements in oscillating motions, the effect of this parameter on stability derivatives of a plunging and pitching airfoil is investigated in this paper.For these 3 conditions, 3 particular mean angles of attack are determined: 1.5°, 11° and 18° for pre-stall, stall and poststall conditions, respectively.As shown in Fig. 8, by increasing the reduced frequency from 0.09 to 0.16, longitudinal stability derivative values for plunging airfoil move toward positive values, which leads to unstable oscillations in higher reduced frequencies.This can be justified by fluid displacement effect on airfoil.On the other hand, the increasing in the reduced frequency does not have significant influence on longitudinal stability derivatives of pitching airfoil.Longitudinal stability of pitch rate is also stable in the studied reduced frequency interval.In plunging motion, by increasing the reduced frequency to 0.06, negative damping effects of wake on moment coefficient become evident and, by increasing further, these effects are more severe.In pitching airfoil, longitudinal stability derivatives in low reduced frequency have a different trend in comparison to plunging airfoil due to the pitch rate.
In Fig. 9, which is stall condition, with mean angle of attack equal to 11°, the situation is as the same as the pre-stall condition.The only notable difference is that the variation of longitudinal stability derivatives of plunging and pitching airfoils is not significant as in the pre-stall condition.
The mean angle of attack equal to 18° shows the post-stall condition in Fig. 10.In this condition, as in the previous ones, longitudinal stability derivatives of plunging and pitching airfoils have negative signs, which means they are stable.However, by increasing the reduced frequency, longitudinal stability derivatives of plunging airfoil move toward the positive values and become more unstable.It is important to note that the variation of longitudinal stability derivatives of plunging airfoil is more than that of the pitching airfoil.In other words, plunging airfoil is more sensitive to reduced frequency than pitching airfoil.
Studying pitch rate in all of the conditions shows the reduced frequency effect on the stability of the oscillation.Figure 8 shows that the pitch rate effect is elevated by increasing the reduced frequency and the negative value of C mq is also increased, which means it is always in stable condition.

CONCLUSION
A numerical method was proposed in this paper to predict the stability derivatives of a NACA 63-415 airfoil under pitching and plunging oscillations and to study the effect of one of the most important oscillation parameters, which is the reduced frequency.This method was a co-located approach and was utilized to solve a moving grid problem.Grid convergence was also examined for numerical method and results indicated that the grid convergence index was below 5%, which demonstrates the high accuracy of the numerical approach.The Reynolds number was 1.6 × 10 6 , linear k -ε model was used as turbulence model and Rhie and Chow method was also devoted to overcome pressure checkerboard problem.
The paper investigated the behavior of stability derivatives in different conditions.Furthermore, it is shown that the stability of a plunging airfoil was decreased by increasing the reduced frequency, primarily in pre-and post-stall conditions.Finally, the research shows that the pitching airfoil stability is relatively independent of reduced frequency especially in low angles of attack due to the existence of the pitch rate in this motion.

Figure 1 .
Figure 1.Control volume and its direction indexes.W and E indicate the center of control volumes on the western and eastern sides of the cell, respectively; w and e indicate the boundaries on the western and eastern sides of the cell, respectively.

Figure 2 .
Figure 2. C-type mesh for numerical study.
h 1 , h 2 2 grid spacing, with h 1 being the finer grid, r is the grid refinement ratio, r = h 2 /h 1 ; p΄ is the order of accuracy of the numerical method and 3 different solutions are needed for this term, which is calculated as(Stern et al. 2001):

Figure 4 .
Figure 4. C m and C l against angle of attack for different grids.Figure 5. C d and C l against angle of attack for different time steps.

Figure 5 .Figure 7 .
Figure 4. C m and C l against angle of attack for different grids.Figure 5. C d and C l against angle of attack for different time steps.Angle of attack [deg]

Figure 9 .
Figure 9. Reduced frequency effect on longitudinal stability derivatives in stall condition.

Figure 10 .
Figure 10.Reduced frequency effect on longitudinal stability derivatives in post-stall condition.

Figure 8 .
Figure 8. Reduced frequency effect on longitudinal stability derivatives in pre-stall condition.

Table 2 .
Grid convergence index and order of accuracy for C m .
Table 3 also shows curve fitting data for stability derivatives in an angle of attack equal to 1.5°.

Table 3 .
Curve fitting data for angle of attack of 1.5°.