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

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 $Re=0.3\times {10}^6$ to $Re=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 $Re\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 $Re=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 $AR=5$ at moderate Reynolds numbers $Re=0.3\times {10}^6$ and $Re=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 $Re=0.3\times {10}^6$ and $Re=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 $AR=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.21m$ and the chord length is $c=0.242m$. 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_{\mathrm{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^{+}\le 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 ·\mathbf{u}=0$$

(1)

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\nabla )\mathbf{u}-\nu {\nabla }^2\mathbf{u}=-\frac{\nabla \mathbf{p}}{\rho }$$

(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):

$${\tau }_{ij}=2{\mu }_t\left(S_{ij}-\frac{1}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)-\frac{2}{3}\rho k{\delta }_{ij}$$

(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 $\omega$, 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{\alpha }}_{max}=1^{\circ }/s$ with a large amplitude of motion, ${\alpha }_m={30}^{\circ }-{35}^{\circ }$, at two distinct Reynolds numbers, $Re=0.3\times {10}^6$ and $Re=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=20s$, there is a substantial drop after reaching $C_{N_{max}}$. This is followed by a more moderate increase of $C_N$ into the higher angle of attack region, i.e., ${25}^{\circ }<\alpha <{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 $Re=0.3\times {10}^6$ and $Re=0.67\times {10}^6$ are shown in Figure 4. OpenFOAM results for $C_N$ at $Re=0.67\times {10}^6$ demonstrate a good agreement with wind tunnel data until $\alpha =15.0^{\circ }$ and between $\alpha =15.0^{\circ }$ and $\alpha =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 $\alpha =19.0^{\circ }$ the difference in the normal force coefficient between the wind tunnel data and CFD results is only $\Delta 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 }<\alpha <{25}^{\circ }$. It is also evident that the level of the abrupt drop in the normal force coefficient at critical angle of attack ${\alpha }_{cr}\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 }<\alpha <{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 $\alpha$ which is identical to the wind tunnel data. The maximum difference in $C_N$ between the top and the bottom branches occurs at $\alpha =20.0^{\circ }$ where $C_{N_{{20}^{\circ },top}}\approx 1.04$ and $C_{N_{{20}^{\circ },bot}}\approx 0.652$ leading to a noticeable drop in the normal force coefficient $\Delta C_N\approx 0.3880$. The static hysteresis loop obtained at this Reynolds number of $Re=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 $Re=0.3\times {10}^6$. The grid used for simulations involving $Re=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 $Re=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 $Re=0.67\times {10}^6$. An approximately $8-10$ degrees wide hysteresis loop was obtained at Reynolds number of $Re=0.3\times {10}^6$. Note that there are plenty of 2D experimental and CFD results for the NACA 0018 airfoil hysteresis at $Re=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_{max}}$ and earlier stall for the lower Reynolds number of $Re=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 $Re=0.3\times {10}^6$ and $Re=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 }<\alpha <{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 $\alpha$ 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 $\alpha =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 $\alpha =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 $\alpha =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 $\alpha =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 $Re=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 $\alpha =10.0^{\circ }$. The maximum asymmetry in both $C_l$ and $C_n$ during the pitch up phase takes place at $\alpha \approx 23.0^{\circ }$ which correlates with the maximum normal force coefficient $C_{N_{max}}$. The peaks of asymmetry develop during the pitch down phase and is highly agitated in the range of ${14}^{\circ }<\alpha <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., $\alpha =10.0^{\circ }$ and $\alpha =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{\alpha }$ 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 ${\alpha }_a=18.0^{\circ }$, the amplitude ${10}^{\circ }$ and three different frequencies $\omega$, $2\omega$ and $3\omega$, where $\omega =0.05$ rad/s. The maximum rate of angle of attack change for results presented in Figure 11 is respectively ${\dot{\alpha }}_{max}=0.5deg/s$ (dotted line), $1.0deg/s$ (dashed line) and $3.0deg/s$ (solid line).

The obtained computational results demonstrate that the simulations conducted at ${\dot{\alpha }}_{max}=0.5deg/s$ and $1.0deg/s$ are suitable for capturing the top branch of static hysteresis, whilst the simulation conducted at ${\dot{\alpha }}_{max}=3.0deg/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{\alpha }$ 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{\alpha }=0.5deg/s$ and $1.0deg/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{\alpha }$ is also significant. More information about the dynamic hysteresis effects in aerodynamic loads due to $\dot{\alpha }$ 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{\alpha }$ 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 ${\alpha }_{*}=\alpha -{\tau }_2\dot{\alpha }$. 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 $\alpha$ and its rate of change $\dot{\alpha }$, but also on state variable x, for example, the normal force coefficient in the stall region can be represented with nonlinear function $C_N(\alpha ,x)$, reflecting the influence of internal variable x on the aerodynamic loads. This general function should have the following properties: $C_N(\alpha ,1)=C_{N_{att}}\left(\alpha \right)$ and $C_N(\alpha ,0)=C_{N_{det}}\left(\alpha \right)$, where $C_{N_{att}}\left(\alpha \right)$ is the dependence of the normal force coefficient under attached flow conditions which is extended above the stall region, and $C_{N_{det}}\left(\alpha \right)$ 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\left(t\right)$ is described by the following first order linear relaxation differential equation:

$${\tau }_1\frac{dx}{dt}+x=x_0(\alpha -{\tau }_2\dot{\alpha })$$

(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\left(\alpha \right)$ characterises the location of flow separation under static conditions, ${\alpha }^{*}=\alpha -{\tau }_2\dot{\alpha }$ is the delayed angle of attack, where the delay is proportional to the rate of change of the angle of attack $\dot{\alpha }$ and ${\tau }_1$ and ${\tau }_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(\alpha ,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(\alpha ,x)=C_{N_{att}}\left(\alpha \right)g\left(x\right)+C_{N_{det}}\left(\alpha \right)(1-g\left(x\right))$$

(5)

where $g\left(x\right)\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\left(x\right)=(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\left(x\right)$ 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 $AR=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({\alpha }^{*},x)\frac{dx}{dt}=\left[x_0\left({\alpha }^{*}\right)-x\right]F_{hm}\left({\alpha }^{*}\right)+F_{sd}\left({\alpha }^{*}\right)$$

(6)

where $F_{hm}\left({\alpha }^{*}\right)$ is the hysteresis morphing function, $F_{sd}\left({\alpha }^{*}\right)$ 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_{hm}\left({\alpha }^{*}\right)$ as a closed curve in the form of an ellipse and the function $F_{sd}\left({\alpha }^{*}\right)$, which is non-zero only in the vicinity of two saddle points of a nullcline set of points defined by equation $[x_0\left({\alpha }^{*}\right)-x]F_{hm}\left({\alpha }^{*}\right)=0$ and function $b({\alpha }^{*},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\left({\alpha }^{*}\right)-x]F_{hm}\left({\alpha }^{*}\right)=0$. The function $x_0\left(\alpha \right)$ is crossing the ellipse defined by the morphing function $F_{hm}\left(\alpha \right)$ in two saddle points $SP1$ and $SP2$. These two saddle points are structurally unstable and functions $F_{sd}\left({\alpha }^{*}\right)$ 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 $FB1$ and $FB2$ 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_{st}\left(\alpha \right)$ for the top and bottom branches of static hysteresis $x_{st-top}\left(\alpha \right)$ and $x_{st-bot}\left(\alpha \right)$. More details about the constructed functions are given in the Appendix A.

The characteristic time scale ${\tau }_1$ in the case of Eqn. 6 is now expressed as:

$${\tau }_1(\alpha ,x)={\left[b(\alpha ,x)\frac{\partial F}{\partial x}{|}_{F=0}\right]}^{-1}$$

(7)

where $F=\left[x_0\left(\alpha \right)-x\right]F_{hm}\left(\alpha \right)+F_{sd}\left(\alpha \right)$. This means that when approaching the bifurcation points of static hysteresis $FB1$ and $FB2$ the relaxation time is approaching infinity ${\tau }_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{\alpha }=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 $FB1$ shown in Figure 16 and (c) the loop on the bottom branch of $C_N$ reaching‘ the bifurcation point $FB2$ 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\left(\alpha \right)$ 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(\alpha ,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 $AR=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 $Re=0.3\times {10}^6$ and $Re=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.