1. Introduction
High-lift wing configuration systems are optimized for takeoff and landing, so the effect of proximity to the ground on their aerodynamics is an important part of the design process [1] [2] [3]. Loss of control in flight (LOC-I) due to aerodynamic stall is most critical in close proximity to the ground due to lack of recovery time from aircraft stall upsets [4] [5]. From another perspective, airframe noise reduction during approach and landing is currently an important problem in aviation and therefore computational prediction of the level of noise for a high lift wing with a deflected slat is now an active research area [6]. Additionally, a better understanding of flow separation processes on high-lift wing configurations in close proximity to the ground is essential for the development of adequate flight simulation models for pilot training and improvement of flight safety [7].
Computational analysis of the flow around two-dimensional airfoils and low aspect ratio wing configurations underground effect conditions has been the focus of attention in the literature [8] [9] [10] [11]. For example, the impact of ground on the aerodynamic performance and stability of the CRM-WB-HVT (Common Research Model-Wing-Body-Horizontal-Vertical Tail) aircraft configuration in both the longitudinal and lateral directions, was studied in [12]. Most of the computational and experimental studies have been focused either on predicting the transition to flow separation regimes or on estimating the level of noise [1] [6].
Investigation of aerodynamic stall and dynamic stall hysteresis phenomenon in free air for various flow conditions with a broad range of Reynolds numbers has attracted significant attention in the literature [13] - [23], partly from the perspective of the development of Reduced Order Models (ROM) based on phenomenological principles [17] [24] [25] [26]. However, there is a lack of studies in stall aerodynamics taking into account the ground effect. The motivation of our study lies in the extension of dynamic stall hysteresis simulation in close proximity to the ground, considering a three-element high-lift 30P30N airfoil.
In this paper, the study focuses on the aerodynamics of the 30P30N high-lift airfoil system in close proximity to the ground including the post-stall region in both static and dynamic conditions. The computational prediction of aerodynamic characteristics of the 30P30N high-lift airfoil system at a moderately high Reynolds number of
and Mach number
was performed using the OpenFOAM code [27] and Fluent software [28] in the angle of attack range covering the stall zone. The experimental results for the 30P30N airfoil [2] are used for validation of the obtained simulation results. The Unsteady Reynolds Averaged Navier-Stokes (URANS) equations in both OpenFOAM and Fluent were solved with the use of the Spalart-Allmaras (SA), one-equation eddy viscosity model [29].
The paper is organized as follows. The computational framework including the mesh generation and blocking techniques, time integration methodology, and numerical setup is described in Section 2. The obtained CFD simulation results for fixed angles of attack and oscillatory-pitching motion together with flow visualizations are presented in Section 3. The concluding remarks and plans for the future work are presented in Section 4.
2. Computational Framework
The section describes the 30P30N airfoil geometry, grid generation techniques, URANS equations, turbulence modelling, and the numerical setup which was employed in all simulations in this study. Identical numerical setups were used for the OpenFOAM code and the Fluent software during CFD simulations for consistency of the obtained results.
2.1. Geometry and Grid Generation
The high-lift airfoil 30P30N shown in Figure 1 is comprised of a slat, main airfoil and the flap section. The geometrical features of this airfoil are shown in Table 1, where c is the stowed chord,
and
are the deflection angle of the slat and flap section and
and
are the gap for the slat and flap. The grids for the 30P30N geometry are generated using a structured blocking approach, with two fluid zones connected by an Arbitrary Mesh Interface (AMI). The AMI with two fluid zones helps to further extend the simulations to dynamic loops in both free air and in close proximity to the ground. The structured blocking approach for the high-lift 30P30N airfoil was utilized such that the cell skewness and the cell non-orthogonality are maintained at good quality in order to avoid divergence and improve the convergence process in the URANS solver. The far field extends 50 chord lengths away from the high lift airfoil in all directions. The cell size transition from the far field of the domain to the airfoil is small enough to avoid large variation of gradients of vector and scalar fields. Furthermore, the gaps between the slat, main airfoil and flap are well refined in order to accurately capture any transition or separation, which is often the case for high-lift configuration airfoils. The boundary layers were modelled with an O-type blocking with dimensionless normal wall distance Y+ less than unity i.e.,
, consisting 40 layers growing in the normal direction with a factor of 1.08. A close-up view of the resulting mesh with 267,392 grid points and 131,856 quadrilateral elements is shown in Figure 2.
The boundary conditions used in the OpenFOAM and Fluent are defined as
,
,
to
, and
and
, where U is the velocity,
is the turbulent viscosity,
is the modified turbulent viscosity, and
is the pressure gradient. The boundary condition on the ground is implemented using a “moving WallBoundary” function in both static and pitching motion cases.
Figure 1. Geometrical aspects of the high-lift 30P30N airfoil.
Table 1. Geometrical aspects of 30P30N airfoil.
Figure 2. Structured grid adopted for the 30P30N high-lift airfoil at
.
2.2. Governing Equations
The Navier-Stokes (NS) equations for the case of incompressible flow include the continuity equation:
(1)
and the momentum equation:
(2)
The computational resources required for direct numerical simulations (DNS) of Equations (1) and (2) usually exceed currently available capabilities. The effect of turbulence is normally simplified by solving the unsteady Reynolds-averaged Navier-Stokes (URANS) equations, which are a time-averaged approximation of the NS equations. The averaging of fluctuating velocities generates additional terms, known as the Reynolds stresses. To describe these stresses, additional empirical equations, either algebraic or differential, are required to close the computational system. The majority of URANS turbulence models are based on the concept of eddy viscosity, equivalent to the kinematic viscosity of the fluid and describing the turbulent mixing or diffusion of the flow momentum [27]. The Reynolds stresses appearing in the URANS equations due to averaging, are normally modelled in the linear approximation with the well known Boussinesq assumption as shown in Equation (3):
(3)
2.3. Turbulence Model
The SA [29], one equation turbulence model, which is commonly used in aeronautical applications and specifically developed for external aerodynamics involving adverse pressure gradients and strongly separated flow conditions [17] is employed in this study. Furthermore, the SA model is computationally cheaper than the Shear-Stress Transport (SST) model [30] which is also widely used for flows with strong separation conditions. The SA model solves for the modified turbulent viscosity
and is described by the partial differential equation in Equation (4).
(4)
The turbulent viscosity is then obtained from
(5)
where
(6)
and
(7)
The SA model constants are given in Table 2. In Equation (4), the production of turbulence is described by the term
(8)
the diffusion of turbulence is described by the term
(9)
and the destruction of turbulence is described by the term
(10)
The Dual Time Stepping (DTS) method is used to simulate both the static and the dynamic loops and the DTS can be described by set of equations as follows [27] [31].
(11)
which
(12)
where is n is the time step level, m is the sub-iteration level within each time step,
is the Runge-Kutta coefficient in each stage and
is the integrated volume of each element. The term
is obtained by using either a) the Stabilized Local Time Stepping Scheme (SLTS) or the Courant Number Limited Euler Method (CoEuler) Scheme, which are already implemented in OpenFOAM. The SLTS scheme guarantees a locally stabilized time step size so that the contribution to FVM matrix due to advective/convective terms remain diagonally dominant and the CoEuler scheme maintains the time step size such that the local flow courant number does not exceed the specified values by the user.
The use of a dual time stepping method with 30 - 50 iterations in the inner loop ensures that the continuity equation is fully satisfied, such that the residuals of the flow field at every time step are dropped at least 5 magnitudes of order, which seemed to be sufficient for this study. Furthermore, with the use of the Courant number based local time stepping technique (CoEuler), the local CFL number was kept below 5, enabling accurate capture of the time scale of any eddies that may be formed during static or dynamic hysteresis conditions.
Table 2. SA turbulence model coefficients.
3. Results
The 30P30N airfoil is first simulated in static conditions with no mesh motion and fixed angle of attack,
, using URANS solver and with varying height to chord ratio (
) from the ground. The “free-air” condition is described as far away from the ground (
), such that there is no influence of the ground on the pressure generated either on the suction or pressure side of the 30P30N airfoil. The
parameter is a dimensionless parameter characterising closeness to the ground. A summary of the investigated cases is presented in Table 3.
Once the analysis of the ground effect characteristics of the 30P30N airfoil in static conditions was completed, the airfoil was then simulated in oscillating pitching motion i.e. dynamic loops with specific reduced frequency, k, where
in a) free-air conditions (
) and b) at
. A summary of the setup for dynamic mesh motion cases is shown in Table 4. In Table 4,
is angular rotation rate in rad/s and c is the stowed chord length in metres.
3.1. Validation of CFD Results
This subsection presents the computational results validation using a) comparison of pressure coefficient, CP, against experimental results from [2]; and b) the lift force convergence study for three fixed angle of attacks i.e. (
and 8˚). The validation is carried out at static conditions with no mesh motion and at the Reynolds number,
with respect to airfoil stowed chord length and Mach number,
.
The pressure coefficient, CP distribution obtained in OpenFOAM at angle of attack,
, is shown in Figure 3. There is an excellent agreement of pressure coefficient variation over the slat, main airfoil section and the flaps. The maximum suction pressure is observed at the very start of the leading edge of the main airfoil with
as shown by experimental results from Klausmeyer [2]. This gives us confidence to move forward with other test cases as planned using the currently adopted structured grid generated for the 30P30N airfoil.
Table 3. Summary of case setup in static conditions with fixed angle of attack,
, for the 30P30N airfoil.
Table 4. Summary of OpenFOAM case setup for dynamic loops with
for the 30P30N airfoil.
Figure 3. Validation of pressure coefficient, CP distribution for 30P30N airfoil at
and experimental results from [16].
The obtained results for the force convergence study are shown in Figure 4. The lift coefficient converges to constant values at flow time of
for all the considered angle of attacks.
3.2. Simulation of the Static Stall in Flow past 30P30N Airfoil
The static variation of angle of attack for flow past high-lift 30P30N is carried out at static conditions with no grid motion and covering attached and separated flow regions. The simulations are performed for in and out of ground effect, using two different CFD codes OpenFOAM and Fluent. The obtained results for the lift coefficients for the considered cases (1a, 1b and 1c) are shown in Figure 5. The
in Figure 5 simply means that the airfoil is far away from the ground such that no influence of ground/runway is on the aerodynamic characteristics of the airfoil.
The results for the clean configuration with
from Figure 5, show attached flow conditions with linear variation of the lift force coefficient CL up to angle of attack,
. From 15 degrees of angle of attack, the
slope starts to decline, indicating increasing separation over the airfoil. The
is at
which is the same as obtained in the experimental results shown in [1]. At
a sharp drop in the lift force coefficient is observed. This happens as the flow separation has reached full development with indication of a significantly strong stall for the 30P30N high-lift airfoil.
Figure 4. Force convergence study of the lift coefficient, CL for 30P30N airfoil at
.
Figure 5. Lift force coefficient, CL, for various h/c ratios at
and experimental results from [16].
The results for the 30P30N airfoil in ground effect with close proximity to the ground, shows a considerably large degradation in the aerodynamic characteristics. Both the OpenFOAM simulations at
and Fluent simulations at
show a large drop in the lift force coefficient values in comparison to the “clean” airfoil with
. The results from OpenFOAM code at
in Figure 5 (see red dashed line with circle markers) show that the reduction in the lift coefficient CL values is bigger in comparison to the CL values obtained for fluent simulations at
. For the ground effect case simulated in OpenFOAM with
, maximum lift coefficient,
, is observed at
with a value of
which is almost a 27% reduction in the maximum lift coefficient.
Most of the computational results reported in literature demonstrated that with the decrease of the distance between the airfoil and the ground, typically denoted by the vertical distance of airfoil’s center of moment from ground to the airfoil chord ratio i.e., (
), an increase in the lift coefficient, CL values are shown [8] [12]. However, for the high-lift configurations there have been noticeable deviations in this trend, showing a decrease of the lift coefficient, CL as the airfoil approaches the ground/runway. Such computational case studies with significant reduction in the lift force were reported in [32] [33]. The DATCOM (Section 4.7) [3], presents results showing that a transport aircraft with a slotted flap wing show a noticeable reduction in the lift force in close proximity to the ground. Paper [34] reports similar kind of decrease of the lift force coefficient, CL in close proximity to the ground, where the pressure loss in the suction side was much bigger the increment of pressure caused due to increased closeness to the ground.
A noticeable reduction in the lift force coefficient values is shown in the work by Parth [35], where a modified discrete vortex method-based model along with the Computational Fluid Dynamics (CFD) was used for the prediction of 2D/3D ground effect for various high-lift airfoil and wing configurations. Within their simulation of 30P30N airfoil, for the angle of attack of
, at flow conditions of Reynolds number,
and
, the drop in the lift coefficient CL between
and
was approximately equal to a significant
. This large drop of the lift coefficient by a factor of roughly 1.47 should be kept in consideration, as significant research is now focused on lift and drag force coefficient optimization for high-lift airfoils. To be more specific, the lift and drag optimizations such as found in [36] should be simulated with consideration of how high-lift airfoils behave in ground effect with close proximity to the runway.
A comparison of the pressure coefficient, CP distribution, for the slat section, main airfoil and the flap section at three different angle of attacks i.e.
and 8˚ are shown in Figure 6 As the angle of attack increases the maximum negative pressure generated on the upper surface of the main airfoil increases significantly. However, with the increase of angle of attack, the CP distribution on the flaps do not change by a large extent, while at
the pressure distribution on the slats increase significantly.
Figure 7 shows the difference in pressure coefficient distribution in and out of
Figure 6. Comparison of pressure coefficient, CP, distribution at
for various angle of attack.
Figure 7. Comparison of pressure coefficient, CP distribution at
and
for angle of attack,
.
ground effect for the 30P30N airfoil. It is evident that there is a slight increase in the integral of the pressure developed on the lower side/pressure side of the airfoil. However, in comparison to the free air simulations with
, the main airfoil loses significant pressure at the leading edge when
. Additionally, it can be noticed that the amount of pressure generated in the flap section of the airfoil is reduced for the considered angle of attack of
.
To understand the effect of ground on the aerodynamic characteristics of the airfoil, the skin friction coefficient, Cf is analyzed for the 30P30N airfoil and the results are shown in Figure 8. The higher skin friction coefficient for the clean airfoil in Figure 8 (blue circles) indicates more attached flow conditions in comparison to the ground effect case with
(see red circles), where a lower skin friction coefficient was observed. A lower skin friction coefficient usually indicates separated flow regions leading to lower lift coefficient values as shown in Figure 5.
3.3. Simulation of Dynamic Stall for Flow past 30P30N Airfoil
In this section, the results for the dynamic loops in and out of ground effect for the 30P30N airfoil are presented. The dynamic loops are carried out at three different frequencies of 0.5, 1 and 1.5 Hz (reduced frequencies of
and 0.031). The DTS method along with a local time stepping for the pseudo time integration is used for time advancement. The cell limited least squares evaluation was employed for flow vector/scalar gradients and second order accuracy for all flow quantities such as velocity and turbulence were maintained. The inner fluid zone made of a smaller circle domain inside the viscous fluid outer domain rotates using an angular oscillatory rotation function, while the outer fluid zone is kept in static non-moving conditions. The mean angle of attack,
for the dynamic loops is
and the amplitude of the pitch angle is
. The frequencies
and 1.5 Hz are representative of angular rotation rates often used in wind tunnel tests for estimation of aerodynamic derivatives of aircraft. Large amplitude oscillations with amplitude of
Figure 8. Comparison of skin friction coefficient, Cf distribution in clean conditions vs ground effect at
.
was chosen in order to cover the stall region and evaluate the dynamic stall effects on the 30P30N airfoil in close proximity to the ground. The dynamic loops for variation of the aerodynamic coefficients are repeated for 3 cycles and the data from the final cycle have been extracted in order to eliminate transients in the flow dynamics. The angle of attack variation in time is described by Equation (13)
(13)
The dynamic loops obtained in OpenFOAM for the free air simulations with
are illustrated in Figure 9. The obtained results indicate a significant increase in the
and a large delay in stall angle
for all the considered cases. When reduced frequency of
, the stall is delayed until
. This is followed by a significant drop in the lift coefficient (see Figure 9 magenta line). The re-attachment of CL does not occur until
.
When the reduced frequency is increased to
, the stall angle
is further delayed to
, with a maximum lift coefficient of
(see Figure 9 blue line). The return branch is compromised of high-frequency vortex shedding and shows no re-attachment of flow back to static conditions until
. When the reduced frequency is further increased to
the return branch shows an even large drop in the lift coefficient. For instance at angle of attack,
, the lift coefficient,
, whereas the in the top branch the lift coefficient was a much higher value of
.
The dynamic hysteresis observed for all the three reduced frequencies shows a
Figure 9. Lift force coefficient, CL, for dynamic stall loops at Reynolds number,
, for 30P30N airfoil and
.
large difference in the CL values between the upstroke and downstroke motion with no indication of reattachment for the entire duration of the return loop.
The dynamic hysteresis loops obtained for the cases with close proximity to the ground are shown in Figure 10. Two frequencies of
and 1.5 Hz are used as the angular rotation rate. The results from Figure 10 show similar characteristics as shown in Figure 9. The dynamic loops lead to a large delay in the stall angle with an increase in
coefficient for the top branch. For the bottom branch the lift coefficients are significantly low and do not show return to static conditions until the angle of attack is equal to
.
A comparison of the obtained dynamic hysteresis loops for free air and in close proximity to the ground is presented in Figure 11. In comparison to the free air simulations (see solid lines in Figure 11), the dynamic hysteresis loops obtained in ground proximity (see dashed lines in Figure 11) indicates lower lift coefficient CL values for the entire range of angle of attack variation.
The flow visualization for
and frequency,
is shown in Figure 12. The top branch shows almost fully attached flow conditions while the bottom branch shows large separation regions above the main airfoil section which corresponds well to the CL values for the dynamic hysteresis loops shown in Figure 9. The contours for turbulent viscosity at two different angle of attacks
and 18˚ in Figure 13, show no eddies or large separation regions during the upstroke, while during the downstroke motion for the same angle of attack, turbulent eddies are shed from the airfoil.
Figure 10. Lift force coefficient, CL, for dynamic stall loops at Reynolds number,
, for 30P30N airfoil and
.
Figure 11. Comparison of lift force coefficient, CL, for dynamic stall hysteresis loops at Reynolds number,
, for 30P30N airfoil with
and
.
Figure 12. Streamlines of velocity for dynamic stall conditions at reduced frequency of
and angle of attack,
, top branch (left) and bottom branch (right).
4. Conclusion
Stall aerodynamics of the three-element airfoil 30P30N has been investigated through CFD simulations using the OpenFOAM software both in free air and in close proximity to the ground at high Reynolds number
and low Mach number
. The simulation results in static conditions show that the ground effect causes the significant reduction in the lift coefficient CL (
). This unusual effect on the lift force due to the ground proximity is likely associated with the influence of a slotted flap on the pressure distribution of the airfoil [3] (Section 4.7). The validation of the CFD simulation results was done for static points only in free air conditions because of the absence of experimental data in close proximity to the ground. Dynamic stall simulations show an additional decrease in the lift coefficient
of the 30P30N airfoil during the downstroke both in free air conditions and in close proximity to the ground due to the development of a fully separated flow, delayed to lower angles of attack. This loss of lift requires special attention from a flight safety point of view, as dynamic stalls can occur during takeoff and landing due to pilot error or gusts of wind. An investigation of the finite aspect ratio wing with the high-lift 30P30N airfoil in free air and in the ground proximity is considered by the authors in the future.
Nomenclature
= normal force coefficient
= lift coefficient
= pitching moment coefficient
= reference chord length and area
= maximum lift coefficient
= reference flow velocity
k = reduced frequency,
Re = Reynolds number
t = physical time
= aspect ratio
= angle of attack
= mean angle of attack
= amplitude of angle of attack variation
= height to chord ratio