Preprint

Article

CFD Simulations and Phenomenological Modelling of Aerodynamic Stall Hysteresis of NACA 0018 Wing

Altmetrics

Downloads

110

Views

Comments

A peer-reviewed article of this preprint also exists.

Mohamed Sereez^*,Nikolay Abramov,

Mikhail Goman^*

Mohamed Sereez^*,Nikolay Abramov,

Mikhail Goman^*

This version is not peer-reviewed

Submitted:

17 February 2024

Posted:

20 February 2024

You are already at the latest version

Alerts

Abstract

Computational simulations of three-dimensional flow around a NACA 0018 wing with an aspect ratio of AR=5 are carried out using the Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations with the Shear-Stress Transport turbulence model closure. Simulations were performed to capture the aerodynamic stall hysteresis using a developed Pseudo-Transient Continuation (PTC) method based on a dual time step approach in the CFD code OpenFOAM. The flow was characterized by an incompressible Mach number of M=0.12 and a moderate Reynolds number Re=0.67×10^6. The results obtained indicate the presence of a noticeable aerodynamic hysteresis in the static dependencies of the force and moment coefficients, as well as the manifestation of bi-stable flow separation patterns, accompanied by the development of asymmetry in the stall zone. The URANS simulation results are in good agreement with the experimental data obtained for the NACA 0018 finite aspect ratio wing from the low-speed wind tunnel under the same test conditions. A new phenomenological bifurcation model of aerodynamic stall hysteresis under static and dynamic conditions is formulated, which is proven to be able to closely match the experimental data.

Keywords:

Subject: Engineering - Aerospace Engineering

1. Introduction

Aerodynamic stall phenomena are important in many engineering applications that influence critical operating limits and system performance. Aerodynamic hysteresis, which is observed both in experiments and in computational simulations under static and dynamic conditions is an important feature of the stall region [1,2,3,4,5,6,7,8,9,10,11]. Stall prediction using wind tunnel tests or computer simulations is challenging due to the high sensitivity of flow separation respectively to experimental testing conditions or computational models and solvers.

Strong turbulence in the wind tunnel flow and vibrations of the test model mounting system can cause premature transitions between the branches of static hysteresis [6]. In computational simulations, there is also an observed sensitivity to the numerical solver, grid resolution, turbulence model, etc. [11]. Albeit the complications, it is important to improve the understanding of the flow physics of this phenomenon to develop more reliable modeling of stall aerodynamics. Applications where static hysteresis exists range from wind turbines and micro air vehicles operating at low Reynolds numbers to the problem of loss of control in-flight (LOC-I) of transport aircraft for which modeling of stall/post-stall aerodynamics at high Reynolds numbers is required [12].

There are many wind tunnel results that show static hysteresis in the separated flow region. For instance, occurrence of static aerodynamic hysteresis in flow past the NACA 0018 airfoil in the range of Reynolds numbers from

R e = 0.3 \times 10^{6}

R e = 1 \times 10^{6}

was shown in [4,5,13]. It was also shown that aerodynamic static hysteresis can exists even at high Reynolds number

R e \approx 6 \times 10^{6}

and may be even enlarged with a circular bump on the leading edge of the pressure side of the airfoil [7]. CFD predictions of such static hysteresis phenomena on 2D airfoils using the URANS equations in OpenFOAM [14,15] with the SA (1-eqn) and SST (2-eqns) turbulence models have also been successful [13,16]. The use of an algebraic turbulence model such as the Baldwin-Lomax also deemed sufficient for capturing static hysteresis on the NACA 0012 airfoil at

R e = 1 \times 10^{6}

, but was associated with a very strong buffeting even after introducing a special stabilizing numerical procedure [10].

In this paper, we consider three-dimensional separated flow with the existence of aerodynamic static hysteresis on NACA 0018 finite aspect wing with

A R = 5

at moderate Reynolds numbers

R e = 0.3 \times 10^{6}

and

R e = 0.67 \times 10^{6}

. For the conducted URANS simulations with the use of SST turbulence model an open source CFD code OpenFOAM was used. The obtained simulation results were compared with wind tunnel data from [8]. To the best knowledge of the authors, this is the first simulation results of aerodynamic static hysteresis loops for a three-dimensional wing configuration matching the experimental data and showing the development of flow asymmetry at high angles of attack. This shows the applicability of URANS simulations for evaluation of separated flows forming with bi-stable flow structures and static hysteresis in aerodynamic loads. The simulations were carried out using a pseudo-transient continuation algorithm based dual time solver developed in OpenFOAM as proposed in [17]. The methodology involved driving the residuals of segregated equations to zero (or at least to a truncation error) in every time step to ensure full convergence of the flow variables [14,15].

This paper also introduces a novel phenomenological method for modelling static hysteresis manifesting bi-stable separated flow structures. A concept of such bifurcation modelling of static hysteresis was proposed in [18] and some of its implementations were discussed early in [13,19,20,21].

The application of the proposed phenomenological model in this paper and its validation against experimental results is shown in Section 4. The CFD simulation results are presented in the paper as follows. The computational framework including the geometry, grid and other numerical setup are discussed in Section 2. Computational results for static hysteresis obtained for Reynolds numbers of

R e = 0.3 \times 10^{6}

and

R e = 0.7 \times 10^{6}

using OpenFOAM are shown in Section 3. This section also includes the skin friction visualization patterns and three dimensional streamlines for bi-stable flow structures. The concluding remarks and future work is presented in Section 5.

2. Computational Framework

Within this section, the discussion revolves around the computational framework employed for simulating static and dynamic hysteresis. This includes the geometrical attributes of the NACA 0018 wing, the governing equations, and the numerical settings applied in OpenFOAM.

The finite aspect NACA 0018 wing with

A R = 5

was made up from the two-dimensional NACA 0018 section profile with rounded tips at its ends. The wing span of the NACA 0018 was measured at

b = 1.21 m

and the chord length is

c = 0.242 m

. A blunt trailing edge was implemented similar to as that in the wind tunnel wing geometry. The incorporation of rounded tips aimed to facilitate a smooth flow attachment to the wing edges and support the proper development of wing tip vortices. The adopted geometry of the considered wing is illustrated in Figure 1. The boundaries of the virtual wind tunnel are positioned at a distance of 50 chord lengths in the upstream, downstream, and sideways directions. Consequently, there is no need for blockage corrections in the Computational Fluid Dynamics (CFD) results obtained in this setup.

The inlet is subjected to a Dirichlet velocity-inlet boundary condition, complemented by a Neumann-type zero-gradient pressure boundary condition. Meanwhile, at the outlet, a zero-gradient velocity outlet boundary condition is applied, and the pressure is prescribed with a fixed value of

p_{outlet} = 0

. Slip boundary conditions are adopted for the wind tunnel walls. To replicate wind tunnel conditions, the turbulent kinetic energy at the inlet is set to a fixed value based on a turbulence intensity of

0.1 %

. For the solid rotating surfaces, specifically the NACA 0018 wing, a "movingWallVelocity" boundary condition is implemented to ensure a zero-flux condition during dynamic or quasi-steady oscillations.

The computational grids are generated using ICEM CFD software, incorporating an O-grid topology for seamless wrapping of the O-type blocking around the wing with a blunt trailing edge. This approach facilitates the creation of a high-quality structured grid with favorable cell determinant values for the hexahedral cells. Utilizing O-type blocking ensures a well-defined boundary layer, maintaining optimal values for cell skewness and orthogonality. The ratios of cell area and volume transitions fall within the range of

0.8 - 1.2

, allowing a maximum change of 20 percent. This ensures that large gradients in flow scalar and vector variables are avoided during the simulation. The boundary layer comprises 30 adjacent layers with a growth rate of

1.15

. The height of the first cell layer is determined by a non-dimensional wall distance of

Y^{+} \leq 1

, enabling the use of a no-wall-function low Reynolds number approach. Following a thorough grid independence study, a mesh size of 3.5 million elements was found to be sufficient for this study as shown in Figure 2.

The flow conditions, characterized by a relatively low Mach number of

M = 0.12

, allows the consideration of incompressible fluid flow. The Navier-Stokes (NS) equations governing incompressible fluid flow are the continuity (1) and the momentum equation (2):

\nabla \cdot u = 0

(1)

\frac{\partial u}{\partial t} + (u \cdot \nabla) u - ν \nabla^{2} u = - \frac{\nabla p}{ρ}

(2)

For flow conditions characterized by high Reynolds numbers, the computational demands associated with the Direct Numerical Simulation (DNS) of Eqns. (1) and (2) typically surpass the current computational capabilities. To address the effects of turbulence in a more computationally feasible manner, the Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations are often employed. These represent a time-averaged approximation of the Navier-Stokes equations, with the averaging process introducing additional terms known as Reynolds stresses. Describing these stresses necessitates the inclusion of empirical equations, either algebraic or differential, to close the computational model. The majority of URANS turbulence models rely on an eddy viscosity concept, analogous to the kinematic viscosity of the fluid, to characterize the turbulent mixing or diffusion of flow momentum. In linear turbulence models, the Reynolds stresses resulting from averaging are modelled using the Boussinesq assumption (3):

τ_{i j} = 2 μ_{t} (S_{i j} - \frac{1}{3} \frac{\partial u_{k}}{\partial x_{k}} δ_{i j}) - \frac{2}{3} ρ k δ_{i j}

(3)

In this study, the SST (Shear-Stress Transport) two-equation turbulence model, as proposed by Menter [22], is utilized. This model is widely applied in aeronautical applications, particularly for external aerodynamics involving adverse pressure gradients and strongly separated flow conditions [11,22]. The authors also have previously employed the SST model for capturing static hysteresis phenomena for the NACA 0018 2D airfoil [13]. The SST model solves two equations, one for the turbulent kinetic energy (k) and the second for evaluating

ω

, which represents the specific dissipation rate of turbulence.

Following a meticulous evaluation and testing of various finite volume schemes and solvers available in OpenFOAM, as outlined in Table 1, the Pre-conditioned Conjugate (PCG) solver coupled with the Geometric Algebraic Multi-Grid method (GAMG) as a pre-conditioner seemed as the most efficient algorithm for driving unsteady residuals to zero at each time level. By employing the GAMG pre-conditioner with 10-30 iterations and applying pre and post-smoothing of residuals for 2-3 levels, it was observed that only 10-30 iterations of the PCG Solver were required to reduce the residual (in the outer iterations) to nearly zero. The gradients of the flow quantities are computed using the second-order accurate Gauss linear scheme with limiters based on cell center values. The divergence of the vector velocity field and the scalar turbulent quantities is estimated with second-order accuracy using the "cellLimited Gauss linear" scheme in OpenFOAM. A linear interpolation scheme is utilized for estimating the contribution of cell center variables to the faces.

The estimation of the time derivatives in OpenFOAM for the case studies involved in this research is accomplished using the dual time stepping method, which has been well established for steady state problems [23,24]. More particularly, the time integration technique proposed in [17] was adopted, which was verified and tested for several dynamic hysteresis and quasi-steady hysteresis cases.

3. Computational Results

In this section, the results of computational simulations for the flow past the NACA 0018 finite aspect wing in two different scenarios are presented and discussed. The first case explores the aerodynamic stall hysteresis obtained during a very slow quasi-steady sinusoidal motion at

{\dot{α}}_{m a x} = 1^{\circ} / s

with a large amplitude of motion,

α_{m} = 30^{\circ} - 35^{\circ}

, at two distinct Reynolds numbers,

R e = 0.3 \times 10^{6}

and

R e = 0.67 \times 10^{6}

. Additionally, insights into the development of flow asymmetry for high angles of attack is also provided. In the second case, the influence of the frequency of the sinusoidal motion on the accuracy of the obtained aerodynamic stall hysteresis is examined. The motion summary for both cases is detailed in Table 2.

3.1. Case1 - Effect of Reynolds Number and Development of Flow Asymmetry

The variation of the normal force coefficient

C_{N}

with physical time is shown in Figure 3. The angle of attack variation during this motion is also mapped on the same Figure. Notably, after a nearly linear increase of

C_{N}

until time

t = 20 s

, there is a substantial drop after reaching

C_{N_{m a x}}

. This is followed by a more moderate increase of

C_{N}

into the higher angle of attack region, i.e.,

25^{\circ} < α < 35^{\circ}

. By comparing the

C_{N}

values at the same angle of attack during the increase and decrease of angle of attack phases, the development of aerodynamic stall hysteresis and bifurcation of aerodynamic loads is distinctly evident in Figure 3.

The acquired computational results for two different Reynolds numbers of

R e = 0.3 \times 10^{6}

and

R e = 0.67 \times 10^{6}

are shown in Figure 4. OpenFOAM results for

C_{N}

R e = 0.67 \times 10^{6}

demonstrate a good agreement with wind tunnel data until

α = 15 . 0^{\circ}

and between

α = 15 . 0^{\circ}

and

α = 20 . 0^{\circ}

the

C_{N}

values (solid black line) slightly exceed the wind tunnel

C_{N}

values (filled black markers). Note that the deviation in

C_{N}

between CFD and wind tunnel data in this range of angle of attack is very small, for instance at

α = 19 . 0^{\circ}

the difference in the normal force coefficient between the wind tunnel data and CFD results is only

Δ C_{N} = 0.035

. The computational simulation results also show that the fully developed stall conditions leading to abrupt drop of

C_{N}

occurs at almost identical angle of attacks i.e.,

23^{\circ} < α < 25^{\circ}

. It is also evident that the level of the abrupt drop in the normal force coefficient at critical angle of attack

α_{c r} \approx 24^{\circ} - 25^{\circ}

is also in agreement with the wind tunnel data from [8]. Some difference between CFD and wind tunnel data occurs at higher angles of attack region

26^{\circ} < α < 40^{\circ}

, where the

C_{N}

coefficient is slightly shifted up for the CFD results having practically the same slope.

During the pitch down motion, the CFD results demonstrate significantly lower

C_{N}

values for the same angle of attack

α

which is identical to the wind tunnel data. The maximum difference in

C_{N}

between the top and the bottom branches occurs at

α = 20 . 0^{\circ}

where

C_{N_{20^{\circ}, t o p}} \approx 1.04

and

C_{N_{20^{\circ}, b o t}} \approx 0.652

leading to a noticeable drop in the normal force coefficient

Δ C_{N} \approx 0.3880

. The static hysteresis loop obtained at this Reynolds number of

R e = 0.67 \times 10^{6}

(see solid black line) is almost as wide as the static hysteresis obtained in wind tunnel (see circle markers). The extent of similarity in the normal force coefficient variation between the CFD and wind tunnel data during the pitch up and the pitch down phases gives a certainty that CFD URANS simulations could be used to capture such non-linear bifurcation phenomena with a good level of accuracy.

OpenFOAM simulations were also conducted at a lower Reynolds number of

R e = 0.3 \times 10^{6}

. The grid used for simulations involving

R e = 0.67 \times 10^{6}

was manipulated in the boundary layer region to adjust the Y+ consistently for the lower Reynolds number flow conditions. The obtained CFD results for

R e = 0.3 \times 10^{6}

shown in Figure 4 (dotted black line) demonstrate an early-onset of fully developed stall conditions with an abrupt drop in the normal force coefficient 2 degrees earlier than that for the case with

R e = 0.67 \times 10^{6}

. An approximately

8 - 10

degrees wide hysteresis loop was obtained at Reynolds number of

R e = 0.3 \times 10^{6}

. Note that there are plenty of 2D experimental and CFD results for the NACA 0018 airfoil hysteresis at

R e = 0.3 \times 10^{6}

[4,5,13], but not for the 3D NACA 0018 wing. The slightly decreased maximum normal force coefficient

C_{N_{m a x}}

and earlier stall for the lower Reynolds number of

R e = 0.3 \times 10^{6}

indicates that the computational results acquired is aligned with the physical intuition.

The pitching moment coefficient

C_{M}

obtained from the CFD URANS simulations using OpenFOAM for two different flow Reynolds numbers are illustrated in Figure 5. Similar to the normal force coefficient

C_{N}

, the obtained CFD results demonstrated noticeable static hysteresis in the pitching moment coefficient

C_{M}

for both cases i.e., Reynolds numbers of

R e = 0.3 \times 10^{6}

and

R e = 0.67 \times 10^{6}

The CFD simulations demonstrate the presence of static stall hysteresis in the aerodynamic coefficients

C_{N}

and

C_{M}

in the range of angle of attack

12^{\circ} < α < 25^{\circ}

. To gain a physical insight in the origin of the static hysteresis phenomenon the flow visualization is extremely important. Figure 6 shows the variation of the skin friction coefficient

C_{f}

contours and the surface streamlines on the NACA 0018 3D finite aspect wing. For easier demonstration, the visualizations have been conducted at specific waypoints in the variation of

C_{N}

against

α

during the pitch up and pitch down phases denoted by A, B, C, D and E as shown in Figure 6. The point A in Figure 6 shows attached flow conditions at

α = 10 . 0^{\circ}

with a minor trailing edge separation as one would expect on this finite aspect ratio wing at the beginning of flow separation. With the increase of angle of attack to

α = 20 . 0^{\circ}

(see point B on Figure 6 ) on the top branch, the surface streamlines seem to reverse from the trailing edge curving towards the central portion of the wing surface, simultaneously the flow from leading edge is still attached and travelling towards the trailing edge of the wing. The reversed flow from the trailing edge and the attached flow from the leading edge meet on a central location on the wing forming a confluence flow towards the middle of the wing where they end in two focuses forming the basement for a 3D arch type vortex.

As angle of attack is further increased to

α = 23 . 0^{\circ}

(see point C in Figure 6) on the top branch, the arch vortex core strengthens, becomes more wider and also deviate towards the opposite ends of the wing in the direction of the wing tips. At this angle of attack, the flow is only attached close to the wing tips on the leading edge and is fully detached from the central part of the leading edge. In the pitch down phase, for the same angle of attack of

α = 23 . 0^{\circ}

(see point D in Figure 6) on the bottom branch, the flow is almost fully detached from the leading edge and the main arch vortex widened and moved close to the wing tips.

The surface streamlines also show a few additional small arch vortices on the top of the wing close to the leading edge. The flow is practically fully separated from the wing at points D and E on the bottom branch of static hysteresis (see Figure 6). Three-dimensional flow topology including the presence of arch type vortices at few selected angle of attacks on the top and bottom branches of static hysteresis is demonstrated in Figure 7.

Flow visualization in wind tunnel using the oil flow images during static hysteresis from [25] shows qualitative agreement with computationally predicted flow patterns. This comparison shown in Figure 8 gives some level of confidence in using URANS simulations for predicting such complex and nonlinear bifurcation phenomena.

An important feature of high angles of attack aerodynamics observed in experiments is the development of aerodynamic asymmetry [8]. The onset of aerodynamic asymmetry occurs in the conducted CFD simulations even with a perfectly symmetric grid. For the simulation carried out at

R e = 0.67 \times 10^{6}

, the onset and development of aerodynamic asymmetry is correlated with the presence of arch vortices on the top branch of static hysteresis and during reattachment flow transition from the bottom branch of static hysteresis. In order to estimate maximum asymmetry from the CFD simulations the rolling moment coefficient

C_{l}

and the yawing moment coefficient

C_{n}

were extracted from the CFD simulations and shown in the Figure 9.

During the pitch up phase, the onset of aerodynamic asymmetry occurs at approximately

α = 10 . 0^{\circ}

. The maximum asymmetry in both

C_{l}

and

C_{n}

during the pitch up phase takes place at

α \approx 23 . 0^{\circ}

which correlates with the maximum normal force coefficient

C_{N_{m a x}}

. The peaks of asymmetry develop during the pitch down phase and is highly agitated in the range of

14^{\circ} < α < 16 . 0^{\circ}

Flow visualizations corresponding to onset of aerodynamic asymmetry are shown on the right side of Figure 10 and the variation of the rolling moment coefficient

C_{l}

vs angle of attack is shown on the left side of the same figure. Point A and point C correspond to the same angle of attack points as shown in Figure 6 i.e.,

α = 10 . 0^{\circ}

and

α = 23 . 0^{\circ}

. Points F and G are new and they correspond to the pitch down phase at positive and negative peaks of aerodynamic asymmetry. The flow visualizations clearly show the process of reattachment of the flow and a highly agitated instantaneous profiles of skin friction contours.

3.2. Case2 - Influence of the Rate of Change of Angle of Attack on Aerodynamic Hysteresis

The effect of different rates of change of the angle of attack

\dot{α}

on the normal force coefficient

C_{N}

around the static hysteresis region is presented in Figure 11. For this analysis the angle of attack was changed periodically with the mean angle of attack

α_{a} = 18 . 0^{\circ}

, the amplitude

10^{\circ}

and three different frequencies

ω

2 ω

and

3 ω

, where

ω = 0.05

rad/s. The maximum rate of angle of attack change for results presented in Figure 11 is respectively

{\dot{α}}_{m a x} = 0.5 d e g / s

(dotted line),

1.0 d e g / s

(dashed line) and

3.0 d e g / s

(solid line).

The obtained computational results demonstrate that the simulations conducted at

{\dot{α}}_{m a x} = 0.5 d e g / s

and

1.0 d e g / s

are suitable for capturing the top branch of static hysteresis, whilst the simulation conducted at

{\dot{α}}_{m a x} = 3.0 d e g / s

slightly over predicts the maximum value of the normal force coefficient. Note that the quasi-static simulations are extremely computationally expensive at very low

\dot{α}

as a transient method must be used which requires the dual time solver to iterate many cycles to reduce residuals to zero. For the bottom branch both

\dot{α} = 0.5 d e g / s

and

1.0 d e g / s

demonstrate identical results and are in good agreement with the experimental data. Figure 11 shows that the impact of the rate of angle of attack change

\dot{α}

is also significant. More information about the dynamic hysteresis effects in aerodynamic loads due to

\dot{α}

can be found in [1] and [2].

4. Phenomenological Bifurcation Model

This section presents a phenomenological model of aerodynamic hysteresis, reflecting the empirical properties of aerodynamic loads in the stall zone. The proposed formulation for the phenomenological model is adapted to the available experimental data [8] and the CFD simulation results presented above.

The formulation of the phenomenological model should reflect the effects observed in experiments and computer simulations of the delay in the onset of flow separation, as well as its reattachment, compared to static conditions. It was shown that the flow separation delay is proportional to the rate of change of the angle of attack

\dot{α}

due to the improvement of the pressure gradient in the region of the leading edge [1]. Flow reattachment when returning to lower angles of attack also exhibits similar behavior. A key feature of the phenomenological model is the need to reflect the existence of bi-stable separated structures that create static hysteresis loops and also dynamic transitions between branches of static hysteresis under unsteady conditions.

The phenomenological model, originally proposed in [18], was based on the use of a first-order nonlinear dynamic system, including a folded equilibrium surface with bifurcation points bounding the static hysteresis zone, representing flows with a non-unique structure. A simplified model [18] for aerodynamic loads taking into account a single separated flow structure in the stall zone, presented in [26], demonstrated good agreement with experimental data obtained in a wind tunnel for a number of delta wings, as well as for aircraft performing Cobra maneuver. The used representation is based on introduction of an internal variable x characterising position of vortex breakdown or flow separation point governed by a first order linear differential equation characterising variation of the internal state variable in transitional motion conditions. The delay of stall was described by using the delay argument

α_{*} = α - τ_{2} \dot{α}

. Recently, this delay-relaxation model attracted attention of researchers and showed its efficiency in many applications ranging from separation flow control and gust alleviation [5,27] to transport aircraft dynamic maneuvers [28].

With introduction of the normalised internal state variable

x \in [0, 1]

the aerodynamic loads depend not only on angle of attack

α

and its rate of change

\dot{α}

, but also on state variable x, for example, the normal force coefficient in the stall region can be represented with nonlinear function

C_{N} (α, x)

, reflecting the influence of internal variable x on the aerodynamic loads. This general function should have the following properties:

C_{N} (α, 1) = C_{N_{a t t}} (α)

and

C_{N} (α, 0) = C_{N_{d e t}} (α)

, where

C_{N_{a t t}} (α)

is the dependence of the normal force coefficient under attached flow conditions which is extended above the stall region, and

C_{N_{d e t}} (α)

represents the dependence of the normal force coefficient under fully separated flow conditions which is extended below the stall region (see Figure 12).

In case of a single flow structure the dynamic behaviour of the internal state variable

x (t)

is described by the following first order linear relaxation differential equation:

τ_{1} \frac{d x}{d t} + x = x_{0} (α - τ_{2} \dot{α})

(4)

where

x \in [0, 1]

is a normalised internal state variable characterising the location of flow separation,

x = 1

corresponds to the conditions of an attached flow and

x = 0

corresponds to the conditions of a completely detached flow,

x_{0} (α)

characterises the location of flow separation under static conditions,

α^{*} = α - τ_{2} \dot{α}

is the delayed angle of attack, where the delay is proportional to the rate of change of the angle of attack

\dot{α}

and

τ_{1}

and

τ_{2}

are the time constants, which characterise the relaxation process and the delay effect, respectively.

The transition of flow from attached conditions to completely detached conditions in the dependence

C_{N} (α, x)

considering that most of the normal force is generated in the area close to the leading edge of the wing, the following approximation is accepted:

C_{N} (α, x) = C_{N_{a t t}} (α) g (x) + C_{N_{d e t}} (α) (1 - g (x))

(5)

where

g (x) \in [0, 1]

. According to the Kirchhoff formula for a fully separated flow with a constant pressure zone behind an airfoil [5,26,28] this function with account of Eqn. 5 is represented as

g (x) = (2 \sqrt{x} + x) / 3

. This formula clearly shows that the normal force coefficient

C_{N}

is mainly generated in the immediate vicinity of the leading edge

x < 0.2

, where about

35 %

of normal force can be lost upon transition to a completely separated flow regime. This is confirmed by the classical distribution of

C_{p} (x)

along the airfoil.

There is a wealth of experimental data showing the existence of static hysteresis in aerodynamic loads within the stall region, such as the experimental data and CFD simulation results for the NACA 0018

A R = 5

presented above in this paper. However, the linear relaxation Eqn. 4 does not allow bi-stable solutions at the same angle of attack, so it must be modified to accommodate the phenomenon of static hysteresis.

To achieve this transformation, we modify Eqn. 4 to the following nonlinear form:

b (α^{*}, x) \frac{d x}{d t} = [x_{0} (α^{*}) - x] F_{h m} (α^{*}) + F_{s d} (α^{*})

(6)

where

F_{h m} (α^{*})

is the hysteresis morphing function,

F_{s d} (α^{*})

is the saddle disturbance function. The choice of these two functions, which we present below, is based on the creation of three different equilibrium solutions in the stall region by shaping the function

F_{h m} (α^{*})

as a closed curve in the form of an ellipse and the function

F_{s d} (α^{*})

, which is non-zero only in the vicinity of two saddle points of a nullcline set of points defined by equation

[x_{0} (α^{*}) - x] F_{h m} (α^{*}) = 0

and function

b (α^{*}, x)

is required for tuning the relaxation time in the phenomenological model.

Figure 13 graphically presents major constructions in designing the nonlinear functions in Eqn. 6. The dashed lines represent the nullcline set of points defined by equation

[x_{0} (α^{*}) - x] F_{h m} (α^{*}) = 0

. The function

x_{0} (α)

is crossing the ellipse defined by the morphing function

F_{h m} (α)

in two saddle points

S P 1

and

S P 2

. These two saddle points are structurally unstable and functions

F_{s d} (α^{*})

applied in the vicinity of saddle points will transform the nullcline set of points into continuous folded curve which shapes a skeletal function for static hysteresis. This curve includes two stable branches of static hysteresis, one on the top with

x > 0.2

and one on the bottom with

x < 0.25

. The intermediate branch (red segment) represents unstable solution of Eqn. 6 which separates regions of attraction of stable branches (two black segments). Two fold bifurcation points

F B 1

and

F B 2

on the skeletal curve indicate jump-like transition from one stable branch of static hysteresis to another stable branch depending on the direction of angle of attack variation. The presented skeletal hysteresis function has been adapted to the experimental positions of

x_{s t} (α)

for the top and bottom branches of static hysteresis

x_{s t - t o p} (α)

and

x_{s t - b o t} (α)

. More details about the constructed functions are given in the Appendix A.

The characteristic time scale

τ_{1}

in the case of Eqn. 6 is now expressed as:

τ_{1} (α, x) = {[b (α, x) \frac{\partial F}{\partial x} |_{F = 0}]}^{- 1}

(7)

where

F = [x_{0} (α) - x] F_{h m} (α) + F_{s d} (α)

. This means that when approaching the bifurcation points of static hysteresis

F B 1

and

F B 2

the relaxation time is approaching infinity

τ_{1} \to \infty

. This explains a significant expansion of the hysteresis loops in unsteady conditions with very slow variation of the angle of attack

\dot{α} = 1 - 3

deg/s (see Figure 11).

Figure 14 shows that the constructed nonlinear functions shown in Figure 13 accurately model the quasi-static hysteresis giving good match with the experimental data. Three different cases of dynamic hysteresis were modelled, (a) the loop around the static hysteresis zone shown in Figure 15, (b) the loop on the top branch of

C_{N}

reaching the bifurcation point

F B 1

shown in Figure 16 and (c) the loop on the bottom branch of

C_{N}

reaching‘ the bifurcation point

F B 2

shown in Figure 17. The images on the left in Figs. Figure 15, Figure 16 and Figure 17 show dynamic variations of the internal state variable

x (α)

with respect to the skeletal function of static hysteresis. The phenomenological modelling results in all three cases (a),(b) and (c) are very close to the experimental wind tunnel data different amplitudes and frequencies of oscillation. The tuning of the nonlinear function

a (α, x)

in the model Eqn. 6 was carried out manually by a trial and error method for an ensemble of the considered processes. A formal mathematical approach for the parameter identification of phenomenological bifurcation models needs further research.

5. Conclusions

The use of the Unsteady Reynolds-Averaged Navier-Stokes equations with the Shear-Stress Transport turbulence model allowed to adequately simulate the aerodynamic characteristics of a NACA 0018 wing with an aspect ratio of

A R = 5

in the stall zone with manifestation of bi-stable separation flow structures and static hysteresis. The flow was characterized by an incompressible Mach number of

M = 0.12

and Reynolds numbers

R e = 0.3 \times 10^{6}

and

R e = 0.67 \times 10^{6}

. The obtained simulation results indicate the presence of a noticeable aerodynamic hysteresis in the static dependencies of the normal force and pitching moment coefficients, which are in good agreement with the wind tunnel data. A qualitative agreement between the CFD obtained surface streamlines on the NACA 0018 wing and the wind tunnel surface oil flow visualization was shown for the top and the bottom branch at various angles of attack. Onset of the aerodynamic asymmetry of separated flow and associated non-zero rolling and yawing moments are also qualitatively similar in CFD simulations and in experiment. The proposed phenomenological model of aerodynamic stall hysteresis under static and dynamic conditions makes it possible to achieve aerodynamic responses very close to the results of wind tunnel tests. The formal procedure for identifying the parameters of the proposed phenomenological bifurcation model requires special effort.

Conflicts of Interest

“The authors declare no conflicts of interest.”

Appendix A

This appendix provides further details for the phenomenological modelling presented in Section 4.

Function

x_{0}

in Eqn. 6 is formulated as:

x_{0} (α^{*}, x) = \frac{1}{1 + {(α^{*} / α_{c})}^{8}}

(A1)

where

α_{c}

is the inflection point of the

x_{0}

curve.

The hysteresis morphing function

F_{h m}

is represented as:

F_{h m} (α^{*}, x) = {(\frac{x - x_{c}}{x_{h}})}^{2} + {(\frac{α^{*} - α_{c} - α_{o f f}}{α_{w}})}^{2} - 1

(A2)

where

x_{c}

and

x_{h}

are the vertical centre position and the vertical semi-diameter of the ellipse,

α_{c}

and

α_{w}

are the horizontal centre position and the horizontal semi-diameter of the ellipse. The ellipse can be shifted from the inflection point vertically and horizontally, for example

α_{o f f}

in Eqn. A2 shifts ellipse horizontally.

To better fit the experimental data it may require a change of the orientation of the ellipse. In the implemented model the hysteresis morphing function was modified as follows:

F_{h m} (α^{*}, x) = {(\frac{x - x_{c}}{x_{h}})}^{2} + γ x (α^{*} - α_{c}) + {(\frac{α^{*} - α_{c} - α_{o f f}}{α_{w}})}^{2} - 1

(A3)

where

γ

is approximately a rotation angle of the ellipse (see Figure 13).

The saddle disturbance function

F_{s d}

is shaped as a spline bell-like function and applied at two saddles

S P 1

and

S P 2

F_{s d} (α^{*}, x) = \pm b e l l ([- 1.0 0 2.0], m x, a r g)

(A4)

where

m x

is maximum magnitude of

F_{s d}

and

a r g

is:

a r g = \sqrt{{(\frac{α^{*} - α_{S P}}{α_{S P}})}^{2} + {(\frac{x - x_{S P}}{x_{S P}})}^{2}}

(A5)

where

α_{S P}

and

x_{S P}

define the location of a saddle point (shown in Figure 13).

References

Ericsson, L.E.; Reding, J.P. Unsteady Airfoil Stall. Technical Report CR-66787, NASA, 1969.
Ekaterinaris, J.A.; Platzer, M.F. Computational Prediction of Airfoil Dynamic Stall. Progress in Aerospace Sciences 1997, 33, 759–846. [Google Scholar] [CrossRef]
Hoak, D.E. The USAF Stability and Control DATCOM. Technical Report TR-83-3048, Air Force Wright Aeronautical Laboratories, 1960.
Timmer, W.A. Two-Dimensional Low-Reynolds Number Wind Tunnel Results for Airfoil NACA 0018. Wind Engineering 2008, 32, 525–537. [Google Scholar] [CrossRef]
Williams, D.R.; Reisner, F., G.; D., Vahl, H.; Strangfeld, C. Modeling Lift Hysteresis with a Modified Goman–Khrabrov Model on Pitching Airfoils. AIAA 2017, 55, 403–410. [CrossRef]
Hoffmann, J.A. Effects of Freestream Turbulence on the Performance Characteristics of an Airfoil. AIAA 1991, 29, 1353–1354. [Google Scholar] [CrossRef]
Svischev, G.P. Investigation of a Low Drag Airfoil with Different Leading Edge Modifications for Increase of Maximum Lift Force. TsAGI proceedings 1946. [Google Scholar]
Zhuk, A.N.; Kolinko, K.A.; Miatov, O.L.; A.N, K. Experimental Investigations of Unsteady Aerodynamic Characteristics of Isolated Wings at Separated Flow Conditions. Preprint TsAGI 1997, pp. 1–56. (In Russian).
Csursia, P.Z.; Siddiqui, M.F.; Runacres, M.C.; De Troyer, T. Unsteady Aerodynamic Lift Force on a Pitching Wing: Experimental Measurement and Data Processing. Vibration 2023, 6, 29–44. [Google Scholar] [CrossRef]
Mittal, S.; Saxena, P. Prediction of Hysteresis Associated with Static Stall of an Airfoil. AIAA 2000, 38, 2179–2189. [Google Scholar] [CrossRef]
Cummings, R.M.; Goreyshir, M. Challenges in the Aerodynamic Modelling of an Oscillating and Translating Airfoil at Large Incidence Angles. Aerospace Science and Technology 2013, 28, 176–190. [Google Scholar] [CrossRef]
Abramov, N.B.; Goman, M.G.; Khrabrov, A.N.; Soemarwoto, B.I. Aerodynamic Modeling for Poststall Flight Simulation of a Transport Airplane. Journal of Aircraft 2019, 56, 1427–1440. [Google Scholar] [CrossRef]
Sereez, M.; Abramov, N.B.; Goman, M.G. Computational Simulation of Airfoils Stall Aerodynamics at Low Reynolds Numbers. Applied Aerodynamics Conference, RAeS, 2016. http://hdl.handle.net/2086/14049.
OpenFOAM. OpenFOAM: The Open Source Computational Fluid Dynamics Toolbox,. http://www.openfoam.com last accessed 8 June 2019.
Chen, G.; Xiong, Q.; Morris.; J., P.; Paterson, E.; Sergeev, A.; Wang, Y.C. OpenFOAM for Computational Fluid Dynamics. Notices of the AMS 2014, 61, 354–363.
Sereez, M.; Abramov, N.M.; Goman, M.G. Prediction of Static Aerodynamic Hesteresis on a Thin Airfoil Using OpenFOAM. Journal of Aircraft 2020. [Google Scholar] [CrossRef]
Sereez, M.; Abramov, N.B.; Goman, M.G. A modified dual time integration technique for aerodynamic quasi-static and dynamic stall hysteresis. Part G: Jnl of Aerospace Eng 2023, 237, 2679–2904. [Google Scholar] [CrossRef]
Goman, M.G. Mathematical description of aerodynamic forces and moments in unsteady regimes of flow with a non-unique structure. TsAGI proc 1983, pp. 29–44. (In Russian).
Abramov, N.B.; Goman, M.G.; Khrabrov, A. N. amd Kolinko, K.A. Simple wings unsteady aerodynamics at high angles of attack: experimental and modelling results. Technical Report paper no: A9936873, AIAA, 1999. [CrossRef]
Abramov, N.B. Modelling of Unsteady Aerodynamic Characteristics for Aircraft Dynamics Applications at High Incidence Flight. Technical Report Thesis, De Montfort University, 2005.
Alieva, D.A.; A., K.K.; N, K.A. Hysteresis of the aerodynamic characteristics of NACA 0018 airfoil at low subsonic speeds. Thermophysics and Aerodynamics 2022, 29, 43–57. [CrossRef]
Menter, F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA 1994, 32. [Google Scholar] [CrossRef]
Jameson, A. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings. Technical Report 10th Computational Fluid Dynamics Conference, AIAA, 1991. [CrossRef]
Nordstrom, J.; Ruggiu, A.A. Dual Time-Stepping Using Second Derivatives. Journal of Scientific Computing 2019. [Google Scholar] [CrossRef]
Kolmakov, Y.A.; Ryzhov, Y.; Stoliarov, G.I.; Tabachnikov, V. Investigation of flow structure of high aspect ratio wing at high angles of attack. TsAGI proceedings 1985, 2290. (In Russian) [Google Scholar]
Goman, M.G.; Khrabrov, A. State-Space Representation of Aerodynamic Characteristics of an Aircraft at High Angles of Attack. Journal of Aircraft 1994, 31. [Google Scholar] [CrossRef]
Sedky, G.; Jones, A.R.; Lagor, F.D. Lift Regulation During Transverse Gust Encounters Using a Modified Goman–Khrabrov Model. AIAA 2020, 58. [Google Scholar] [CrossRef]
Luchtenburg, D.M.; Rowley, C.W.; Lohry, M.; Stengel, R.F. Unsteady High-Angle-of-Attack Aerodynamic Models of a Generic Jet Transport. Journal of Aircraft 2015, 52, 890–895. [Google Scholar] [CrossRef]

Figure 1. NACA 0018 wing geometry.

Figure 2. Normal force coefficient

C_{N}

variation against grid resolution parameter h from the conducted mesh independence study.

Figure 2. Normal force coefficient

C_{N}

variation against grid resolution parameter h from the conducted mesh independence study.

Figure 3. Normal force coefficient

C_{N}

(primary y axis) and angle of attack

α

(secondary y axis) variation with physical time for the static hysteresis loop obtained at

R e = 0.67 \times 10^{6}

Figure 3. Normal force coefficient

C_{N}

(primary y axis) and angle of attack

α

(secondary y axis) variation with physical time for the static hysteresis loop obtained at

R e = 0.67 \times 10^{6}

Figure 4. Normal force coefficient

C_{N}

α

during the pitch up and pitch down phases: CFD simulations and wind tunnel data from [8].

Figure 4. Normal force coefficient

C_{N}

α

during the pitch up and pitch down phases: CFD simulations and wind tunnel data from [8].

Figure 5. Pitching moment coefficient

C_{M}

α

during the pitch up and pitch down phases.

Figure 5. Pitching moment coefficient

C_{M}

α

during the pitch up and pitch down phases.

Figure 6. Contours of skin friction

C_{f}

superimposed on surface streamline patterns on the NACA 0018 wing during the pitch up and pitch down phases of aerodynamic static hysteresis.

Figure 6. Contours of skin friction

C_{f}

superimposed on surface streamline patterns on the NACA 0018 wing during the pitch up and pitch down phases of aerodynamic static hysteresis.

Figure 7. 3D streamlines of velocity field and contours of skin friction

C_{f}

superimposed on surface streamline patterns on the NACA 0018 wing during the pitch up and pitch down phases at

α = 20 . 0^{\circ}

and

23 . 0^{\circ}

Figure 7. 3D streamlines of velocity field and contours of skin friction

C_{f}

superimposed on surface streamline patterns on the NACA 0018 wing during the pitch up and pitch down phases at

α = 20 . 0^{\circ}

and

23 . 0^{\circ}

Figure 8. Qualitative comparison of wind tunnel oil flow visualization [25] and CFD surface streamline patterns on the top and bottom branches of aerodynamic hysteresis.

Figure 9. Rolling and yawing moment coefficients (

C_{l}

C_{n}

) from CFD simulations vs

α

during the pitch up and pitch down phases.

Figure 9. Rolling and yawing moment coefficients (

C_{l}

C_{n}

) from CFD simulations vs

α

during the pitch up and pitch down phases.

Figure 10. Contours of skin friction

C_{f}

superimposed on surface streamline patterns on the NACA 0018 wing during the pitch up and pitch down phases of aerodynamic static hysteresis showing development of flow asymmetry.

Figure 10. Contours of skin friction

C_{f}

superimposed on surface streamline patterns on the NACA 0018 wing during the pitch up and pitch down phases of aerodynamic static hysteresis showing development of flow asymmetry.

Figure 11. Normal force coefficient

C_{N}

from computational simulations at different rotational frequencies with

α_{0} = 18 . 0^{\circ}

and amplitude

α_{m} = 10 . 0^{\circ}

Figure 11. Normal force coefficient

C_{N}

from computational simulations at different rotational frequencies with

α_{0} = 18 . 0^{\circ}

and amplitude

α_{m} = 10 . 0^{\circ}

Figure 12. Normal force coefficient

C_{N}

of NACA 0018 wing at

R e = 0.67 \times 10^{6}

: 1) static test points, 2) slow pitch up and pitch down measurements and 3) approximation lines for fully attached and fully detached flow conditions,

C_{N_{a t t}}

and

C_{N_{d e t}}

Figure 12. Normal force coefficient

C_{N}

of NACA 0018 wing at

R e = 0.67 \times 10^{6}

: 1) static test points, 2) slow pitch up and pitch down measurements and 3) approximation lines for fully attached and fully detached flow conditions,

C_{N_{a t t}}

and

C_{N_{d e t}}

Figure 13. Phenomenological model skeletal function of static hysteresis based on experimental data.

Figure 14. Normal force coefficient

C_{N} (t)

for a very slow change in angle of attack based on phenomenological modeling and experimental results at

R e = 0.67 \times 10^{6}

Figure 14. Normal force coefficient

C_{N} (t)

for a very slow change in angle of attack based on phenomenological modeling and experimental results at

R e = 0.67 \times 10^{6}

Figure 15. Internal state variable

x (α)

and normal force coefficient

C_{N} (α)

: comparison of phenomenological results with experimental data at

α (t) = 25^{\circ} + 24^{\circ} s i n (2 π f t)

where

f = 0.2

Hz.

Figure 15. Internal state variable

x (α)

and normal force coefficient

C_{N} (α)

: comparison of phenomenological results with experimental data at

α (t) = 25^{\circ} + 24^{\circ} s i n (2 π f t)

where

f = 0.2

Hz.

Figure 16. Internal state variable

x (α)

and normal force coefficient

C_{N} (α)

: comparison of phenomenological results with experimental data at

α (t) = 12^{\circ} + 13^{\circ} s i n (2 π f t)

where

f = 0.8

Hz.

Figure 16. Internal state variable

x (α)

and normal force coefficient

C_{N} (α)

: comparison of phenomenological results with experimental data at

α (t) = 12^{\circ} + 13^{\circ} s i n (2 π f t)

where

f = 0.8

Hz.

Figure 17. Internal state variable

x (α)

and normal force coefficient

C_{N} (α)

: comparison of phenomenological results with experimental data at

α (t) = 25^{\circ} + 7^{\circ} s i n (2 π f t)

where

f = 1.6

Hz.

Figure 17. Internal state variable

x (α)

and normal force coefficient

C_{N} (α)

: comparison of phenomenological results with experimental data at

α (t) = 25^{\circ} + 7^{\circ} s i n (2 π f t)

where

f = 1.6

Hz.

Table 1. Finite volume solvers of OpenFOAM used in the study

Method	equations	preconditioner	tolerance	iterations
Pre-conditioned Conjugate Gradient (PCG)	pressure-correction	Multi-Grid GAMG	1e-06	50
smoothSolver	Momentum/ turbulence	none	1e-10	10

Table 2. Summary of the considered case studies.

case	k	M	Re	$α_{0}$	Amplitude, $α_{m}$
1	$0.01$	$0.12$	$0.3 - 0.67 \times 10^{6}$	$0^{\circ}$	$30^{\circ} - 35^{\circ}$
2	varying	$0.12$	$0.67 \times 10^{6}$	$18^{\circ}$	$10^{\circ}$

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.