1 Introduction

Flow over oscillating airfoil is encountered in many systems such as sail surfaces, compressors, wind turbines, micro air vehicles, and aircrafts. Therefore, many researchers have taken into account different aspects to investigate the flow over oscillating airfoils considering inviscid and viscous flows. This also includes the unsteady flow over airfoils with plunge, pitch, combined pitching, and plunging motion. Furthermore, previous studies have considered the unsteady flow from low to high Reynolds numbers and from subsonic flows to transonic flows using numerical and experimental methods. Gnesin and Rzadkowski^[1] solved the Euler equations using the Godunov-Kolgan volume scheme to investigate the three-dimensional inviscid flutter. The dynamic stall was investigated in the study by Ekaterinaris and Platzer^[2]. The unsteady, compressible, and viscous fluid flow over the NACA 0012 airfoil was studied by Murthy et al.^[3]. They used the Roe scheme to explore the laminar separation. Lee and Gerontakos^[4] used hot-film sensors to study boundary layer transition and separation over an oscillating airfoil experimentally. Furthermore, they discussed mechanisms responsible for the stall events. The two-dimensional heaving airfoil in a low Reynolds number was studied by Lewin and Haj-Hariri^[5], in which the effects of oscillation frequency on the flow characteristics and power coefficient were investigated. Gopalakrishnan and Tafti^[6] used a collocated finite volume multi block approach to study the incompressible viscous flow over the heaving airfoil. The effects of compressible flow on the heaving airfoil were investigated by Raveh and Dowell^[7] considering self-sustained shock oscillation in a transonic flow over the NACA 0012 airfoil at reduced frequency and angle of attack.

Akbari and Price^[8] studied the laminar incompressible viscous flow over an elliptic airfoil with pitching oscillation at a large angle of attack. They used the Joukowski transformation to solve the Navier-Stokes equations in the computational domain to investigate the influence of the oscillation frequency, the location of pitch axis, and the thickness ratio of ellipse. In their next study^[9], they investigated the unsteady flow field over the pitching NACA 0012 airfoil using the vortex method to solve the Navier-Stokes equations. They concluded that the reduced frequency has a significant effect on the flow field and force coefficient of the airfoil. The influence of the pitch angle on the vertex structure and aerodynamic forces was explored in the study of Sarkar and Venkatraman^[10], where they studied the dynamic stall in a wide range of the reduced frequencies. In their next study^[11], they discussed the thrust generated by the pitching airfoil at high frequencies of oscillation and investigated the effects of various parameters such as the amplitude of oscillation and the location of pitching on thrust generation. Fluid flow with low Reynolds numbers over the pitching airfoil was studied by Amiralaei et al.^[12]. Their study indicated the effects of aerodynamic parameters on the maximum lift coefficient and hysteresis loops of lift and drag coefficients. Also, fluid flow with high Reynolds numbers over a symmetric rotor blade with pitching oscillation was investigated by Sarkar and Bijl^[13]. They studied the effects of initial conditions and structural nonlinearity. Transonic flow over the pitching airfoil was studied experimentally by Hillenherms et al.^[14].

Guglielmini and Blondeaux^[15] investigated the flow over a foil with a combined pitching and plunging motion, where they solved the vorticity equation for a wide range of parameters. The study by Xiao and Liao^[16] investigated the effects of angle of attack profile on pitching/plunging motion. They achieved a modified cosine profile for the angle of attack. Furthermore, the particle image velocimetry (PIV) was used by Prangemeier et al.^[17] to study the flow over an airfoil with plunge and quick-pitch motion. They modified the pitching/plunging oscillation using various quick-pitch motion. Ashraf et al.^[18] numerically studied the effects of airfoil thickness and camber on the combined pitching and plunging airfoil considering laminar and turbulent flow regimes. Moreover, the transonic flow over the thin airfoil with pitching and plunging motion was studied by Guruswamy and Yang^[19], and the transonic flow over heave/pitching flutter was experimentally investigated by Dietz et al.^[20]. In addition, the Navier-Stokes equations for the transonic flow over the NACA 0012 airfoil with pitch and flap motion were solved by Iovnovich and Raveh^[21]. They discussed that the aerodynamic coefficients can be regulated by planned structural motion.

The incompressible viscous flow around a flapping airfoil with active chordwise flexibility was investigated by Tay and Lim^[22], where they studied the effects of flexing on three airfoil configurations and showed that flexibility does not necessarily improve the airfoil performance. Flow with low Reynolds numbers over a pitching/plunging airfoil with nonlinear deformation was studied by Jaworski and Gordnier^[23]. They explored the effects of the Reynolds number variation on the vortex dynamics and aerodynamic coefficients and discussed the differences between a rigid airfoil and an airfoil with low to moderate oscillation amplitudes. Unger et al.^[24] investigated the laminar-turbulent transition of incompressible fluid flow over a flapping flexible airfoil, where the effects of trailing edge deformation on the lift coefficient and the location of transition were investigated. They compared the numerical results obtained for the jig-shape of airfoil with the corresponding results obtained by the PIV measurements. The turbulent incompressible flow over an airfoil with a vertical displacement, rotation, and deformed position for aileron was studied by Svacek^[25]. In this study, a reliable and robust numerical model for the simulation of fluid-structure interaction problems was presented using the finite element method and the k-ω turbulence model.

The incompressible flow over the pitching/plunging airfoil has been completely investigated in the above studies. However, the effects of airfoil deformation in the compressible fluid flow particularly at the transonic regime remain mostly unknown. Here, the unsteady transonic inviscid compressible flow over the two-dimensional deformable NACA 0012 airfoil with plunging motion is investigated numerically to explore various aspects of flow patterns. Our study covers a wide range of various parameters involved in the unsteady transonic flow. To do so, the Roe scheme in a generalized coordinates system is developed to solve the governing arbitrary Lagrangian-Eulerian equations. The effects of different governing parameters such as the phase angle, the deformation amplitude, the initial angle of attack, the flapping frequency, and the Mach number on the unsteady flow field and aerodynamic coefficients are studied in detail.

2 Governing equations

Taking into account the unsteady inviscid compressible flow over a deformable airfoil with plunging motion, the two-dimensional Euler equations in the generalized coordinates (ξ, η, t) are used. To write the equations in the non-dimensional form, the following non-dimensional parameters are used:

where x is the dimensionless coordinate in the x-direction, y is the dimensionless coordinate in the y-direction, c is the airfoil chord, t is the dimensionless time, u is the dimensionless velocity in the x-direction, v is the dimensionless velocity in the y-direction, x_t is the dimensionless grid speed in the x-direction, y_t is the dimensionless grid speed in the y-direction, u_∞ is the free stream velocity, ρ is the dimensionless density, T is the dimensionless temperature, p is the dimensionless pressure, ρ_∞ is the free stream density, and T_∞ is the free stream temperature.

Therefore, the non-dimensional Euler equations in the conservative form can be written in the generalized coordinates as

where where ξ_x, ξ_y, η_x, and η_y are the matrices of transformation, and J is the Jacobian of transformation. The vectors of conserved variables Q and the convective fluxes F and G are given by

The contravariant velocities, U and V , are defined as

Σ is the sum of kinetic and internal energies per unit volume, which can be expressed in the non-dimensional form as where γ is the heat capacity ratio, and Ma_∞ is the Mach number of free stream, The equation of state for a perfect gas in the non-dimensional form is defined as In the present study, the flow over a deformable NACA 0012 airfoil with an initial angle of attack (α₀) with the plunging motion is considered. The plunging motion, shown in Fig. 1, is expressed as follows: where h is the dimensionless instantaneous position of airfoil, h₀ is the dimensionless stroke amplitude, and ω is the dimensionless flapping frequency, ω = ωc/u_∞. To facilitate the analysis of deformation amplitude, the deformation of airfoil is simplified using a time variant profile assumption, i.e., where $x'$ and $y'$ indicate the coordinates attached to the airfoil body, a₀ is the dimensionless deformation amplitude, and ψ is the phase angle.

Aerodynamic characteristics consisting of the lift, drag, moment, and pressure coefficients are, respectively, as

where f_x and f_y are the non-dimensional forces applied on the airfoil surface in the x- and y-directions, respectively. Moreover, m_c/4 is the non-dimensional torque around a point at a quarter of chord from the leading edge, which is given by the following equation: Finally, the averaged lift, drag, and moment coefficients are expressed as where T_p is the non-dimensional flapping period, which is given by 2π/ω. The aerodynamic performance of the deformable airfoil can be defined by its lift-to-drag ratio ${\bar C_1}$/${\bar C_d}$ ^[26] .

3 Numericalmethod

In the present study, the Roe scheme in the generalized coordinates^[27] is developed to solve the arbitrary Lagrangian-Eulerian fluid flow equations. Figure 2 shows the grid in both physical and computational domains, where the quadrilateral grids are used to maintain the non-diffuseness of the Roe scheme.

Here, the cell-centre method is used to store the physical quantities as shown in Fig. 3, in which the cell faces are described by dashed lines represented by E (east), W (west), N (north), and S (south). The cell faces are located at the midpoint between the nodes, and the unit vector $\hat n$_ξE normal to the cell face E can be expressed as

where θ^E is the angle between the normal vector to the face E and the x-axis. Similarly, $\hat n$_ηN is the normal vector to the face N, which is determined as where βN is the angle between the normal vector to the face N and the x-axis. The angles θE and βN are computed by the following equations: where ξx, ξy, ηx, and ηy are the matrices of transformation. The geometric properties of the grid at the cell face are specified as

Using the first-order upwind method, we can write

where q represents a primitive variable. Therefore, Roe’s averaged conditions are where

Considering Fig. 3, the Roe’s numerical flux in the generalized coordinates crossing the face E is obtained by the following equation:

where $\hat \lambda $_E^(k) are the eigenvalues, $\hat T$_E^(k) are the eigenvectors, and δw_E^(k) are the wave amplitudes. However, in a moving mesh, the eigenvalues are changed, while the eigenvectors and wave amplitudes remain unchanged with respect to the static mesh. $\hat \lambda $_E^(k) , $\hat T$_E^(k) , and δw_E^(k) are specified by Roe’s averaged conditions as follows: where

As shown in Fig. 4, u^⊥_E and u_||^E are the components of fluid velocity vectors vertical and parallel to the face E, respectively. Considering $\hat c$_E as the speed of sound, we can write

Also, x_{t⊥ E} , the mesh velocity component vertical to the face E, is determined as

To avoid the expansion sonic waves, the following entropy correction^[28] is used:

The Roe’s numerical flux in the generalized coordinates at the north face G1N is determined by the following equation:

The eigenvalues, eigenvectors, and wave amplitudes at the cell face N are computed by the similar equations used for the cell face E using the β angle in the relations. By computing numerical fluxes, the Euler equation in the generalized coordinates can be discretized explicitly as follows:

4 Boundary conditions

In the subsonic flow, among the four eigenvalues in the two-dimensional flow, three of them are positive, and one is negative. Hence, at the inflow boundary, there are three incoming and one outgoing characteristics. As a result, three variables consisting of the velocities in the xand y-directions and the density (ρ) are specified, but the pressure (p) is extrapolated from the interior. At the outflow boundary, there is one incoming characteristic and three outgoing characteristics. Therefore, only one variable, the pressure, is specified, and the other variables are extrapolated from the interior using an explicit procedure. In the inviscid flow, the slip condition at the solid boundary is imposed. Therefore, the value of relative velocity (V_r) at the solid boundary, which is tangent to the surface as shown in Fig. 5, is assumed to be equal to its value at the point near the boundary. Consequently,

where ϕ_i is the local surface slope at the node i.

The pressure at the surface is specified by projecting the momentum equation normal to the surface, which results in the value of pressure gradient at the boundary.

where

To determine the variables at the wall, the following iterative method^[29] is used:

(i) Assume ρi, 1 = ρi, 2.

(ii) Calculate u and v using Eqs. (29) and (5).

(iii) Determine the pressure at the surface by solving Eq. (30) implicitly.

(iv) Extrapolate the temperature from the interior. (v) Use the equation of state for a perfect gas, Eq. (8), for density computing.

(vi) Check the convergence as follows:

where ε is the small number, which is assumed to be equal to 0.001 in this study.

5 Dynamic mesh algorithm

One of the important aspects of flow over a deformable airfoil with plunging motion is to handle the airfoil deformation in the computational domain. Here, we use a dynamic mesh approach to resolve the plunging and deformable motion of the airfoil. To do so, an o-type orthogonal structured mesh^[29] is used, where the orthogonality of grid lines at the surface is imposed by employing a forcing function. In addition, the grid is clustered near the surface and in the vicinity of leading and trailing edges using appropriate functions. The dynamic mesh algorithm presented by Batina^[30] is applied, where the spring stiffness for each edge i-j is inversely related to the length of the edge as follows:

New locations of the interior grid points are specified by solving static equilibrium equations at each time step iteratively. This is done using a predictor-corrector method, where the displacements of grid points are predicted by linear extrapolation, specified by

Then, the displacements are corrected by solving the static equilibrium equations in the xand y-directions which are given by Eq. (35). These equations can be solved by the Gauss-Seidel iteration method.

The summations in Eq. (35) are taken over 4 edges surrounding the node i. Finally, the new locations of grid points are specified by the following equations:

The grid points at the outer boundary are held fixed, while the grid points near the solid boundary up to a specified η (here 7 15ηmax) are moved by the oscillating function defined for the airfoil in Eqs. (9) and (10). Therefore, the spring method is not applied in these nodes, but it is applied for the nodes staying between 7 15 ηmax and ηmax (the outer boundary), which enables large movements to be taken into account at the higher h₀ and a₀ without any numerical problems. A sample of dynamic mesh used in this study is shown in Fig. 6 in eight different time intervals within a cycle of motion.

To avoid errors due to mesh movement, the geometric conservation law (GCL) is satisfied during the solution. Use the following expression^[30]:

It can be discretized as The summation of Eq. (38) is taken over 4 edges of the control volume. In each time step, the calculated cell areas are compared to satisfy the GCL condition.

6 Grid independency and code validation

To show the grid independency of the numerical code developed here, three different mesh sizes at h₀ = 0.1, a₀ = 0.1, α0 = 2°, and ψ = π/2 are considered. The lift coefficients with respect to time for eight cycles of motion at Ma_∞ = 0.755 and ω = 0.162 8 for three different meshes are presented in Fig. 7. The comparison shows good agreement of the results obtained in three different mesh sizes.

To verify the accuracy of the developed numerical method and the validation of the numerical code, the unsteady inviscid compressible flow over the NACA 0012 airfoil pitching at Ma_∞ = 0.755, α0 = 0.016°, αm = 2.51°, and ω = 0.162 8 (α = α0+αmsin(ωt)) is solved. The numerical results obtained are compared with the experimental data and the numerical results of Batina^[30] in Fig. 8, where the curves of the lift and moment coefficients versus the instantaneous angle of attack are shown. As it is evident from the figure, the numerical results obtained in the current study are in satisfactory agreement with the previous experimental and numerical results.

7 Results and discussion

Here, the effects of different parameters such as the phase angle (ψ), the dimensionless deformation amplitude (a₀), the Mach number (Ma_∞), the dimensionless frequency (ω), and the initial angle of attack (α0) are investigated in the unsteady transonic compressible flow over a deformable airfoil with plunging motion in a wide range. The dimensionless plunging amplitude (h₀) is fixed at 0.1, and in each case, eight cycles of motion are calculated to obtain a steady state periodic solution. Table 1 shows the effects of phase angle on the aerodynamic coefficients. Since the phase angle has slight influence on the lift and drag coefficients, it is fixed at π/2 for the remaining cases.

In order to investigate the effects of deformation amplitude on the flow over the airfoil, six different values for a₀ are considered from 0 to 0.25 at the intervals of 0.05. Figure 9 indicates the lift, drag, and moment coefficients versus the instantaneous angle of attack at various deformation amplitudes and flow conditions when α0 = 2°, ω = 0.162 8, and Ma_∞ = 0.755. As shown in this figure, increasing the deformation amplitude causes the maximum amounts of lift, drag, and moment coefficients at the same angle of attack to increase. Consequently, the loops are enlarged for the greater a₀.

Time variations of C_l, C_d, and C_m at various deformation amplitudes are shown in Fig. 10. As it is evident from the figure, the instantaneous lift coefficients of the rigid airfoil and the deformable one vary inversely in comparison with each other, i.e., the lift coefficient for the deformable airfoil is decreased while increased for the rigid one. Furthermore, Fig. 10 indicates that the amplitudes of lift, drag, and moment coefficients are increased by the growth of the deformation amplitude. In Fig. 11, the aerodynamic coefficients of the deformable airfoil with plunging motion at different deformation amplitudes are compared with those obtained for the rigid plunging airfoil a₀ = 0. It is shown that the averaged lift coefficient of the deformable airfoil is decreased by the growth of deformation amplitude. Therefore, the deformation has a negative effect on the airfoil aerodynamic performance. This effect intensifies when the deformation amplitude increases.

The effects of deformation motion on the pressure distribution on the airfoil are displayed in Fig. 12, where the pressure coefficients on the upper and lower surfaces of the airfoil with the deformation amplitudes of 0, 0.1, and 0.25 at α0 = 2°, ω = 0.162 8, and Ma_∞ = 0.755 are shown. As it can be seen from the figure, when the deformable airfoil with plunging motion is bent upward through the down stroke, the shock weakens on the upper surface due to the decreasing angle of attack. However, shock is relocated on the lower surface by increasing the negative angle of attack. As the airfoil moves towards its lowest position from (7+1/4)T_p to (7+1/2)T_p, the airfoil returns to its initial shape with the positive angle of attack without deformation. However, the shock remains on the lower surface and transfers from the lower surface to the upper surface with a time lag. Consequently, the lift coefficient on the large portion of down stroke is negative due to the relatively high pressure on the upper surface. As the airfoil deforms downward when bending in up stroke, shock reappears on the upper surface and strengthens by increasing the deformation. Therefore, the lift coefficient becomes positive as the shape of the airfoil reverses. Moreover, Fig. 12 shows that shocks generally strengthen as the deformation amplitude increases. Furthermore, the shock displaces remarkably over the airfoil. It is observed that the pressure difference between the lower and upper surfaces is increased as the deformation amplitude becomes larger.

The computed Mach contours around the deformable airfoil with plunging motion at the largest deformation amplitude, a₀ = 0.25, are shown in Fig. 13 at eight different non-dimensional time within one cycle of motion. It can be seen that at t = (7+1/2)T_p, with the positive angle of attack, shock appears on the lower surface, and at t = (7+1/8)T_p, with the negative angle of attack, shock appears on the upper surface. As discussed previously, this is related to the time lag of shock displacement between upper and lower surfaces as the sign of angle of attack is changed. Moreover, at the large deformation of the airfoil, the shock appears on both lower and upper surfaces at the same time within a complete cycle of motion. This event is also shown in Fig. 12.

To investigate the effects of the initial angle of attack, the time histories of lift, drag, and moment coefficients are displayed in Fig. 14 at various initial angles of attack at a₀ = 0.1, ω = 0.162 8, and Ma_∞ = 0.755. As it can be seen from the figure, increasing the initial angle of attack causes the lift coefficient amplitude to decrease but causes the drag coefficient amplitude to increase. Moreover, the peak value of lift coefficient lags as the initial angle of attack decreases.

In general, at the small initial angles of attack, the airfoil deformation controls the instantaneous angle of attack, and the frequency of the drag coefficient is almost twice as much as the one in the large initial angle of attack at which the instantaneous angle of attack is almost not affected by the airfoil deformation. The computed aerodynamic coefficients of the airfoil with the deformation amplitude of 0.1 at various initial angles of attack are compared in Fig. 15 to those obtained for the rigid airfoil at the same conditions. It can be seen that the averaged lift coefficient of the airfoil with deformation motion is smaller than the rigid airfoil. Furthermore, it is observed that the aerodynamic performance of the airfoil with deformation at small angles of attack is much less than the one of the rigid airfoil. However, the aerodynamic performances of the deformable airfoil and the rigid one get closer to each other as the initial angle of attack increases.

The effects of the flapping frequency on the aerodynamic coefficients are shown in Fig. 16 at a₀ = 0.1, α0 = 2°, and Ma_∞ = 0.755. As it is evident from this figure, at ω = 0.2, there is a transition point in the lift and drag coefficients trend, which appears as an inflection point in the aerodynamic performance coefficient in Fig. 17. Similarly, the computed averaged aerodynamic coefficients at various frequencies are shown in Fig. 17. It can be seen that increasing the frequency causes the averaged lift coefficient to increase and the averaged drag coefficient to decrease. Consequently, the aerodynamic performance of the airfoil improves considerably as the frequency increases. Furthermore, the slope of aerodynamic performance diagram for ω < 0.15 is slightly larger than that for ω > 0.2.

Finally, the effects of the Mach number on the aerodynamic coefficients are investigated. The time histories of lift, drag, and moment coefficients at various Mach numbers for a₀ = 0.1, α0 = 2°, and ω = 0.162 8 are displayed in Fig. 18. It is observed that at the Mach numbers equal or greater than 0.9, the amplitudes of lift, drag, and moment coefficients are decreased. In Fig. 19, the computed averaged aerodynamic coefficients of the deformable airfoil at various Mach numbers are compared with the corresponding values of the rigid airfoil at the same flow conditions. As it is evident from the figure, the trend in the lift coefficient of the deformable airfoil in comparison with the rigid airfoil is not monotonic. For the Mach number less than about 0.88, the lift coefficient of the rigid airfoil becomes larger than the one of the deformable airfoil. However, beyond that, the reverse effect occurs. In addition, the aerodynamic performance of the deformable airfoil degrades compared with the rigid airfoil at Ma_∞ < 0.9. However, it gets closer for the two airfoils at Ma_∞ ≥ 0.9.

As a final point, compared the deformable airfoil and the rigid one in Fig. 15 and Fig. 19, it can be concluded that the deformable airfoil with plunging motion has more smooth and stable aerodynamic performance. Consequently, the deformable airfoil improves maneuverability.

8 Summary and conclusions

In this article, we study the unsteady transonic inviscid compressible flow over a deformable NACA 0012 airfoil with the plunging motion. The Roe scheme in the generalized coordinates system is developed to solve the governing arbitrary Lagrangian-Eulerian fluid flow equations. A structured dynamic mesh is used to resolve the instantaneous motion of the airfoil. In a relatively wide range, the effects of parameters such as the phase angle, the deformation amplitude, the initial angle of attack, the flapping frequency, and the Mach number are investigated.

The numerical results obtained indicate that the phase angle has a negligible effect on the amplitudes and the average lift and drag coefficients. Furthermore, the deformation motion degrades the aerodynamic performance of the airfoil, and the aerodynamic performance weakens by increasing the deformation amplitude. However, the aerodynamic performance of the deformable airfoil with plunging motion stays monotone, indicating its better maneuverability compared with the rigid airfoil. Our results indicate that the aerodynamic performance of the deformable airfoil with plunging motion significantly improves by increasing the flapping frequency. In addition, at small initial angles of attack, the aerodynamic coefficient of the airfoil is remarkably affected by the airfoil deformation due to its dominant effect on the instantaneous angles of attack. However, at large initial angles of attack, the airfoil deformation has a very small effect on the instantaneous angles of attack.