Near Stall Unsteady Flow Responses to Morphing Flap Deflections

Abstract

The unsteady flow characteristics and responses of an NACA 0012 airfoil fitted with a bio-inspired morphing trailing edge flap (TEF) at near-stall angles of attack (AoA) undergoing downward deflections are investigated at a Reynolds number of 0.62 × 10⁶ near stall. An unsteady geometric parametrization and a dynamic meshing scheme are used to drive the morphing motion. The objective is to determine the susceptibility of near-stall flow to a morphing actuation and the viability of rapid downward flap deflection as a control mechanism, including its effect on transient forces and flow field unsteadiness. The dynamic flow responses to downward deflections are studied for a range of morphing frequencies (at a fixed large amplitude), using a high-fidelity, hybrid RANS-LES model. The time histories of the lift and drag coefficient responses exhibit a proportional relationship between the morphing frequency and the slope of response at which these quantities evolve. Interestingly, an overshoot in the drag coefficient is captured, even in quasi-static conditions, however this is not seen in the lift coefficient. Qualitative analysis confirms that an airfoil in near stall conditions is receptive to morphing TEF deflections, and that some similarities triggering the stall exist between downward morphing TEFs and rapid ramp-up type pitching motions.

1. Introduction

Most natural flyers (insects, birds, bats, or even extinct Pterosaurs) used shape changing, or morphing, mechanisms in order to harness the full capabilities of their wings [1,2] (Figure 1). Inspired by nature, the aerospace industry has recently seen a surge in the research and development of bio-inspired, intelligent, and adaptive design techniques for the next generation of aircrafts [3]. The promising bio-inspired concept of morphing structures is used on wings not only to optimize their shape during flight, but also to actively interact with the surrounding flow field in order to mitigate possible aerodynamic performance degradation, reduce the profile drag, delay onset of stall, or achieve early recovery from flow separation. However, there are many technical challenges in the integration of morphing configurations with an existing ‘rigid’ wing; one in particular is the difficulty of simulating and resolving the complex flow physics around morphing structures [4,5,6,7]. In addition, one would imagine the difficulty of finding the optimal morphing arrangements (e.g., frequency, amplitude, angle) whilst solving such numerically intensive simulations. William et al. [8] affirmed that a deeper understanding of fundamental flow responses to active morphing actuations is a prerequisite for improved flow control in unsteady aerodynamics.

In order to investigate the effects of morphing on the aerodynamic performance of airfoils and wings, various studies have focused on the numerical modeling of morphing wing flaps using low- to high-fidelity computational fluid dynamics (CFD) methods. Often, however, such studies simplify the problem into a series of steady-state CFD analyses in a ‘frozen’ position during the morphing, neglecting the dynamic effects that exist during such transient morphing processes.

Using a hybrid RANS-LES approach, Abdessemed et al. [9,10] investigated the aerodynamic performance of an airfoil fitted with a downward, statically morphing trailing-edge flap (TEF) at low angles of attack (AoA) and compared it with a traditional hinged TEF configuration. The results showed that on average, about a 6.5% increase in the lift-to-drag ratio could be achieved with the morphing configuration. Similar conclusions were made for a 3D straight wing embedded within a portion of a statically morphing TEF with a seamless side-edge transition to the wing ends, as investigated by the same authors [11].

The use of continuous actuation by means of morphing surfaces (e.g., [12]) or vibrating flaps (e.g., [13]) has been the subject of a few recent studies. For example, the ability of periodic surface morphing to mitigate flow separation and reduce drag at a low Reynolds number of 50,000 was demonstrated experimentally [12,13]. In a different experiment at a high Reynolds number of 1 million, Jodin et al. [14] showed that a substantial reduction in flow instabilities could be obtained using a vibrating TEF, in addition to an increase in the aerodynamic efficiency, if the flap vibrations are set at optimal frequencies and amplitudes. The use of a harmonically morphing TEF was also explored numerically for both aerodynamic and aeroacoustic performance enhancements [15,16] at low AoAs. It was shown that up to a 3% increase in the aerodynamic efficiency could be obtained, together with a 1.5 dB reduction of tonal noise. These works and several others have demonstrated the ability of periodic morphing actuation to maintain flow reattachment and improve aerodynamic efficiency.

Nevertheless, Medina et al. [17] argued that the transient (unsteady) responses to pulse-like actuators could be better utilized to produce higher lift than the one achieved with periodic forcing, basing these conclusions on an earlier work performed by Amitay and Glezer [18]. Based on these assumptions, a series of experiments to investigate the transient flow responses to low-amplitude high deflection rates of conventional hinged flaps were conducted [17,19]. In order to observe the lift coefficient response to various deflection frequencies, an NACA 0006 airfoil equipped with a TEF (at 50% chord) was tested in a water tunnel, first at 0° and 20° AoA [19] and then at a fixed AoA of 20°, where the 0° setting corresponds to attached flow, while the 20° setting represents massively separated flow. Observations showed that a spike, whose amplitude was proportional to the flap pitching rate, was captured in the lift coefficient time history, although there was no clear ‘time lag’ in the lift response to the flap pitching motion. The effect of the direction of the flap deflection (either upwards or downwards) was also investigated in both studies and it was found that the lift and drag forces were both highly dependent on the pitching direction, especially in the separated flow regime. It should be noted that these water tunnel experiments were conducted at a low Reynolds number (~40,000), which is an order of magnitude lower than that of practical wind tunnel tests. Thus, higher Reynolds number investigations need to be carried out to verify the effectiveness of such a mechanism in producing lift transients that could be useful in practical applications, such as the enabling of aggressive manoeuvres, gust loading alleviation, and stall delay at high AoA.

Morphing flaps could potentially offer a better alternative to the hinged flaps used by Medina et al. [17] in their experiments, e.g., the implementation of morphing could result in a lower weight penalty. In addition, morphing flaps would cover a larger deflection range compared with traditional hinged flaps; these are limited by the actuators’ mechanisms, which often constrain them to a smaller range of motion. Furthermore, exploring large-amplitude and reasonably fast flap deflections at a higher Reynolds number and at AoA away from stall would allow a better understanding of the flow response in conditions closer to real life. However, at angles of attack near stall, an airfoil is always susceptible to early stall, e.g., in the form of dynamic stall, because of rapid variations in AoA or freestream turbulence, resulting in a number of dynamic effects, such as increased unsteadiness in aerodynamic loads, bursting of leading edge vortices, and shedding of trailing-edge vortices. Therefore, studying the near-stall flow regime using unsteady morphing is of particular practical interest and could provide a pre-stall control mechanism to achieve increments in aerodynamic performance to extend the operational envelope, or at least to mitigate against abrupt transitions to stall (e.g., due to sudden inflow perturbations) where control becomes more challenging and less energy-efficient.

To the authors’ best knowledge, the study of near-stall dynamic flow response to a rapid morphing flap deflection rather than a hinged flap deflection has not yet been addressed in open literature by means of either experimental or numerical work at high or low Reynolds numbers. Therefore, this study is a timely addition to the recent works on morphing flaps. The paper presents an exploratory CFD study, using a hybrid RANS-LES model, investigating the responses of both lift and drag forces to rapid downward morphing TEF deflection at various frequencies. An NACA 0012 airfoil fitted with a morphing TEF (of a length equal to 25% of the chord) is investigated at a chord-based Reynolds number of 0.62 × 10⁶. CFD analysis is performed using a newly developed computational framework based on dynamic meshing and unsteady parametrization [10]. The hybrid RANS-LES turbulence model, SBES, is used to resolve the unsteady flow field with reasonable accuracy at an affordable computational cost. Again, to the authors’ best knowledge, this is the first time SBES has been used to investigate the near-stall aerodynamic performance of an airfoil in response to rapid morphing deflections. For a fixed large morphing amplitude, flow responses at three angles of attack near stall (8°, 10°, and 12°) and at a range of morphing frequencies are presented and discussed.

2. Problem Definition

All test cases involve an NACA 0012 airfoil with a chord c = 0.2286 m, a Reynolds number of Re = 0.62 × 10⁶ based on the chord, and a Mach number of 0.115, in a free steam flow velocity of $U_{\infty }=40 \mathrm{m}/\mathrm{s}$ with standard sea-level conditions, matching conditions reproducible in most low-speed wind tunnels and used in recent morphing experimental studies [20].

The parametrization approach used for the TEF motion is a modified version of the method introduced by [9,10] and based on the FishBAC concept [20]. It consists of the baseline NACA 0012 thickness distribution added to an unsteady parametrization of the camber line. Interested readers are referred to reference [10] for the mathematical model of the parametrization. Figure 2 illustrates the deforming motion of the TEF, where x_s is the morphing starting location (i.e., 75% c in this study); w_te is the maximum deflection at the trailing edge, which indicates the vertical distance between the initial TE position and the maximum (final) position of the morphed TE$; t_{start}$ is the morphing starting time; and $t_{morph}$ is the time to reach the final (maximum) deflection position.

The cases presented in this paper all deal with a large deflection angle of the TEF, i.e., for w_te = 5% of the chord, which is about 14° of the local deflection compared to the baseline configuration of w_te = 0.

Table 1. Summary of the cases investigated.

wte = 0.05%c, AoA = 8°
Frequencies (Hz)	1	2	4	6
wte = 0.05%c, AoA = 10°
Frequencies (Hz)	2	4	6	8
wte = 0.05%c, AoA = 12°
Frequencies (Hz)	2	4	6	8

summarizes the parameters used for the study.

3. Governing Equations

3.1. Fluid Dynamics

The governing equations are the conservation of mass and the conservation of momentum. Since the flow velocities in the domain are much smaller than the speed of sound, it can be assumed that the density remains constant throughout the flow field. Therefore, the incompressible pressure-based solver included in ANSYS Fluent is used to solve the Navier–Stokes equations [21]. However, the formulation of the conservation equations is different when it comes to moving boundaries problems, such as for a morphing surface. Equation (1) is the integral form of the conservation equation for a general scalar quantity φ on an arbitrary control volume V with moving boundaries:

$$\frac{\partial }{\partial t}{{\displaystyle \int }}_V\rho \phi dV+{{\displaystyle \int }}_{\partial V}\rho \phi \left(\overrightarrow{u}-{\overrightarrow{u}}_g\right)d\overrightarrow{A}={{\displaystyle \int }}_{\partial V}\Gamma \nabla \phi ·d\overrightarrow{A}+{{\displaystyle \int }}_v{S}_{\phi }dV$$

(1)

where ρ is the fluid density, Γ is the diffusion coefficient, $\overrightarrow{u}$ is the time-averaged flow velocity vector, ${\overrightarrow{u}}_g$ is the grid velocity vector, $S_{\phi }$ is a source term, and $\partial V$ represents the boundary of the control volume V.

Turbulence Modeling

Given the high Reynolds number in this study, it would be impractical to use computationally intensive approaches such as DNS or LES. At the same time, the potential for early stall, with possible flow separation, vortex growth, and shedding, necessitates the use of a model beyond simple RANS. Therefore, a hybrid RANS-LES model, which offers a balanced approach between accuracy and computation, is used. The hybrid RANS-LES model named the stress-blended eddy simulation, or SBES, introduced by Menter [22] is used in this work to provide closure to the Reynolds-averaged Navier–Stokes equations. The use of hybrid RANS-LES models in aerodynamic studies has previously been reported in the literature [23,24,25]. More recently, Syawitri el al. [26] and Rezaeiha et al. [27] both used SBES to predict the aerodynamic performance of vertical axis wind turbines (VAWTs). More importantly, Syawitri el al. [26] showed that it is possible to significantly improve the numerical prediction of the flow around a 2D VAWT with three straight blades using SBES, while still accurately predicting the flow field. Rezaeiha et al. [27] demonstrated that 2.5D SBES is a superior approach to both URANS and scale-adaptive simulation (SAS). Finally, Abdessemed et al. [16] conducted a comparative study between 2D SBES and 3D LES for the same configuration and Reynolds number (0.408 × 10⁶), with the results showing that the 2D SBES approach was able to accurately capture the loads and flow features with a fraction of the numerical cost, which is sufficient for the current study to obtain trends and inform future works.

The SBES was used in conjunction with the shear stress transport (SST) k-ω turbulence model for the RANS region, while for the large eddy simulation (LES) region, the wall-adapting local eddy–viscosity (WALE) sub-grid model [28] is considered.

The blending function is the same as that used in the shielded delayed eddy simulation (SDES) [21]. Moreover, a shielding function (Equation (2)) is used to explicitly switch between the RANS and LES models:

$${\tau }_{ij}^{SBES}=f_{SDES}{\tau }_{ij}^{RANS}+\left(1-f_{SDES}\right){\tau }_{ij}^{LES}$$

(2)

where ${\tau }_{ij}^{RANS}$ is the RANS portion and ${\tau }_{ij}^{LES}$ is the LES portion of the modeled stress tensor; $f_{SDES}$ is the shielding function [22]. Unfortunately, the exact formulation of the shielding function is proprietary and to date has not been published by ANSYS for public use.

Amongst the advantages of SBES is that it gives explicit control over which part of the flow the LES is applied to, it provides a rapid transition from the RANS region to LES region, and has less dependency on the mesh compared with the SDES model [22]. The RANS wall boundary layer regions are protected against influences from the LES model when the shielding functions are in use, which protects against early switching to the LES. If the early switching occurs, it can cause a strong decline in the RANS capabilities [29]. The high-performance capabilities of SBES have made it one of the most promising hybrid RANS-LES models available [29].

3.2. Dynamic Mesh Update Methods

Using the parametric method introduced in Section 4, the dynamic mesh update methods included in ANSYS Fluent are utilized in order to deform the mesh and the geometry simultaneously, whilst maintaining a high-quality mesh in the process. Diffusion-based smoothing was used in the present work because it is more robust when it comes to mesh quality preservation (as required for moving boundaries) compared with spring-based smoothing [21]. The exact setup in [16] was replicated using diffusion-based smoothing to deform the computational grid, and for this specific case with large deformation a limited patch of triangular cells was introduced near the TE (Figure 3) in order to use the remeshing feature included in ANSYS Fluent [21]. Readers are directed to [16,30] for in-depth details of the dynamic mesh setup.

In order to drive the dynamic meshing schemes in ANSYS Fluent, a user-defined function (UDF) was developed to incorporate the unsteady parametrization method explained in Section 4. The UDF makes use of the DEFINE_GRID_MOTION macro embedded in Fluent, while the full source code of the UDF is made available in the following repository [31].

4. Mathematical Model of the Trailing-Edge Motion

The parametrization method used for the TEF motion is a modified version of the one previously introduced in [9,10] and repeated here in Equations (3) and (4) below. It consists of the baseline NACA 0012 thickness distribution as defined in Equation (3) [32] added to the unsteady parametrization of the camber line as defined in Equation (4).

$$y_t=\left(\frac{th}{0.2}\right)\left(0.2969\sqrt{\overline{x}}-0.1260\overline{x}-0.3516{\overline{x}}^2+0.2843{\overline{x}}^3-0.1015{\overline{x}}^4\right)$$

(3)

$$y_c=\{\begin{array}{cc}0& , 0\le \overline{x}<x_s and t_{start}\le t\le \frac{T}{4}\hfill \\ \frac{-w_{te}\sin \left(\frac{2\pi \left(t-t_{start}\right)}{T}\right){\left(\overline{x}-x_s\right)}^3}{{\left(1-x_s\right)}^3}& , \overline{x}\ge x_s\hspace{1em} and t_{start}\le t\le \frac{T}{4}\hfill \\ \frac{-w_{te}{\left(\overline{x}-x_s\right)}^3}{{\left(1-x_s\right)}^3}& , \overline{x}\ge x_s\hspace{1em} and t>\frac{T}{4}\hfill \end{array}$$

(4)

where $\overline{x}$ is the nondimensional airfoil coordinate $\overline{x}=x/c$, t is time, and T is the period of the airfoil’s trailing-edge motion; the morphing frequency is obtained by $f=\frac{1}{T}$. In the current study, at t = $t_{start}$ the morphing starts and the flap is deflected downwards until it reaches the maximum deflection value w_te at $t_{morph}=\frac{T}{4}$, effectively simulating one-quarter of a period, as shown in Figure 2.

For convenience, the non-dimensional time $t^{\prime }$ is used in subsequent analysis, where $t^{\prime }=\frac{\left(t-t_{start}\right)U_{\infty }}{c}$, which means the morphing starts from $t^{\prime }=0$.

5. Computational Setup

A fractional step method of the non-iterative time advancement (NITA) scheme was used for time advancement, as it results in significant CPU time savings [33]. A least-squares cell-based spatial discretization approach was used to calculate all gradients. Moreover, the momentum equation discretization employs a bounded central differencing scheme, with a second-order scheme for all remaining quantities and derivatives.

In order to avoid the appearance of highly skewed cells around the TE region during the morphing, a small region around the TE was created to accommodate an unstructured local mesh by making use of the remeshing option in ANSYS Fluent (see Figure 3).

Table 2. Mesh characteristics before and after mesh deformation, where the range of the metrics indicates the minimum and maximum values of each metric and the solver requirements indicate the minimum value that the mesh metrics need to satisfy.

Mesh Metric	Baseline	Deformed at Max Deflection	Difference	Range of the Metrics	Solver Requirements
Average Orthogonal Quality	0.9624	0.9479	−2%	0 (worst) to 1 (best)	the closer to 1 the better
Minimum Orthogonal Quality	0.0609	0.06052	−1%		Larger than 0.01
Maximum Skewness	0.6825	0.7031	3%	0 (best) to 1 (worst)	Lower than 0.98

summarizes the key mesh characteristics before and after the mesh deformation, showing that all of the values are within the bounds of a well-preserved mesh in terms of mesh quality.

The flow domain consists of an NACA 0012 airfoil with a sharp TE, while the pressure far field is placed 30 chord lengths (30c) away from the TE. A structured O-grid-type mesh was generated around the airfoil (Figure 4). The mesh had a total number of 350,000 cells, with 1400 points around the airfoil surface and 200 points in the streamwise wake direction clustered toward the TE. The inflation layers were refined to achieve a first layer height in wall unit Δy⁺ smaller than 1. These mesh characteristics were based on the fine mesh previously used by the authors in [34].

A second-order transient formulation was used for all transient simulations. A time step of Δt = 10⁻⁵ s (which guarantees a Courant number of an order of unity) was used, and all simulations were run for at least 0.8 s (i.e., equal to 4 through-flow times), which was found sufficient to achieve statistically converged results. The simulations were performed using a 16-core Intel Xeon E5-2680 v2 (Ivy Bridge) with 2.8 GHz CPU, and each simulation needed at least 43 h.

6. Results and Discussion

6.1. Verification and Validation

For verification purposes, XFOIL [35] with viscous corrections was used to provide lift and drag predictions for the baseline NACA 0012 airfoil. Figure 5 shows the numerical results obtained for a range of AoA using ANSYS Fluent and XFOIL solvers, compared against the experimental data of Sheldahl et al. [36] for both the lift coefficient (C_L) and drag coefficient (C_D) for a range of AoAs. Overall, it can be seen that a good agreement is achieved between Fluent and XFOIL’s predictions and the experimental data of Sheldahl et al. [36] in the linear region up to AoA = 10°, whereby the average difference between the two numerical predictions is about 5%. The Fluent predictions are in better agreement with the test data in the linear region for both lift and drag predictions, although the experimental data of Sheldahl et al. [36] seem to show an earlier stall for the airfoil at about AoA = 11°, while the Fluent predictions show a higher stall angle at AoA = 13°, resulting in an 8.6% higher C_L,max than the measurement. XFOIL does not capture the stall till AoA = 16° and consistently overpredicts C_L from AoA = 4° onwards.

A few possible reasons could explain those discrepancies observed between the two numerical predictions and the test data. In the experimental work performed by Church et al. [37] investigating an NACA 0021 airfoil at lower Reynolds numbers (less than 200,000), for the same Reynolds number it was suggested that Sheldahl et al.’s data [36] undermeasured both the stall angle and the C_L,max value. It was concluded that Sheldahl et al.’s experiment is not representative of the NACA 0021 airfoil, mainly because the data prior to stall and early stall were calculated using an in-house PROFILE computer code [36]. Furthermore, experimental data obtained by Ladson et al. [38] showed that for a higher Reynolds number of Re = 2,000,000, Sheldahl et al.’s [36] data undermeasured the stall angle by at least 2° AoA, similar to the findings of the present numerical study. The prediction of the C_L,max value by ANSYS Fluent is generally higher than for quasi-3D wind tunnel tests, which can be explained by the fact that using a purely 2D configuration normally results in overprediction of the lift force. Similar trends can be observed for the drag comparisons, where there is good agreement between Fluent, XFOIL, and test data up to the pre-stall region (i.e., AoA = 8°), after which the discrepancies increase monotonically until AoA = 14°. Beyond this point, both Fluent and the test data show a steep increase in drag, indicating that deep stall occurs.

6.2. Morphing TEF Deflection

Once validation was established, the next task was to investigate the dynamic flow characteristics in response to the TEF deflections. Figure 6 shows the transient lift and drag coefficient responses to morphing TEF deflecting at various frequencies. For all simulation cases, unsteady CFD runs are performed to reach a statistically converged status: after a precursor run of the unmorphed baseline airfoil until t = 0.50 s (t′ = 0), morphing commences and transient flow simulations continue for four through-flow times to reach a statistically converged status.

Table 3. Summary of frequencies used, morphing lengths, and times when the morphing stopped.

Frequency (Hz)	Morphing Time Tmoprh (s)	Morphing Stop Time	$\mathbf{Non}\text{-}\mathbf{Dimensional} \mathbf{Time} \left(\frac{\mathit{t}{\mathit{U}}_{\mathit{\infty }}}{\mathit{c}}\right)$
1	0.25	0.75	43.74
2	0.125	0.625	21.87
4	0.0625	0.5625	10.93
6	0.042	0.542	7.29
8	0.03125	0.53125	5.46

summarizes the frequencies used, and the related time quantities associated with each frequency. These frequencies were chosen to reflect the morphing times used in flaps deployments and in other morphing studies [19].

At two moderate AoA values of 8° and 10°, both the lift and drag coefficients exhibit similar transient behaviour. For a fixed AoA, the lift increases linearly from a converged static value up until the flap reaches its final position where the lift asymptotes to a steady-state value, which was the same for all studied frequencies. The slope at which the lift increases during the morphing motion is proportional to the frequency; the higher the frequency, the steeper the slope. However, the transient drag response shows a noticeable difference compared with the lift response. An overshoot is observed in the drag time history before the flap reaches its final position; the amplitude of the spike observed is also proportional to the frequency, with higher morphing frequencies causing larger overshoots of the spike. At AoA = 12°, the flow response is found to be different from the lower AoAs studied. At this AoA, the NACA 0012 airfoil is already in a pre-static stall region, so that a further increase of the local camber due to the TEF deflection coupled with its rapid rate of motion causes the airfoil to enter into a deep stall regime. Prior to the initiation of morphing, both the lift and drag coefficients display similar trends, where a vortex sheet forms in the wake (Figure 8) of the baseline airfoil at a shedding frequency of 50 Hz (obtained with a fast Fourier transform (FFT)), corresponding to a Strouhal number (S_t) of 0.1 based on the maximum airfoil thickness. Similar to the case of moderate AoAs of 8° and 10°, as the morphing flap deflects from its baseline position, a linear increase in the lift and the drag coefficients is observed, with a slight overshoot in the C_D values at t′ = 3.50. However, the airfoil does not recover to the steady state (as it does at 8° and 10°), causing the airfoil to experience dynamic stall, as revealed by the large-scale flow oscillations (S_t = 0.3) and as exhibited in the lift and drag coefficients around higher mean values of both C_L and C_D.

The aerodynamic efficiency (denoted by C_L/C_D) as illustrated in Figure 7 provides an additional understanding of the aerodynamic load response. Here, two representative cases are shown, AoA = 10° and 12°, both at a morphing frequency of 6 Hz; other frequencies exhibited similar behaviour. Starting from the baseline value for both AoAs, the aerodynamic efficiency begins to monotonically decrease until t′ = 1.75, where the rate of decrease is accelerated due to the overshoot in C_D. For AoA = 10°, the efficiency recovers to higher values at t′ = 9.62, where it settles down at an average value of C_L/C_D = 26 after it experiences a 30% instantaneous drop at t′ = 3.50. In contrast, at AoA = 12° after the morphing motion stops at t′ = 7.35, the aerodynamic efficiency does not recover and continues to decrease further until the dynamic stall process has started. Further insight into the flow behaviour during the rapid morphing flap deflection can be gained by qualitative examination of the flow.

Figure 8 shows the vorticity contours for the airfoil during the morphing process at AoA = 10° and a morphing frequency of 8 Hz. At this AoA and before the start of morphing, the airfoil experiences moderate TE flow separation, which causes the shed vortices to travel downstream, forming a vortex shedding pattern in the wake (Figure 9). This change in the local flow field, which induces local separation, is one of the contributors to the decrease of efficiency that can be observed in Figure 7.

Once the morphing motion starts, the deformed TEF geometry accelerates the vortex shedding process at the TE and further induces a large-scale roll-up vortex that travels upstream. This traveling vortex interacts with the thickening shear layer at the start position of the morphing flap (i.e., at 75%c). The interaction results in a large-scale shed vortex which convects downstream along the flap surface; this could be the cause of the spike observed in the drag coefficient. Once this initial large-scale vortex has convected into the wake, the same cycle repeats, however the resulting vortices are smaller than the initial starting vortex.

Medina et al. [19] observed a spike in their lift measurements, for which the flap motion time was of an order of one convective time or less $(~\frac{t{U}_{\infty }}{c}$), whereas in this study the morphing flap motion time was much larger (between 5 to 20$\frac{t{U}_{\infty }}{c}$), which could explain the absence of the overshoot in the lift in the current study. Additionally, Medina et al. [17] only studied the lift response but not the drag response as was the case in this work, so comparing the drag response between the two cases was not possible.

With morphing motion, the lift and drag responses at the higher AoA of 12° are found to be quite different from their behaviours at the slightly lower AoA of 10°. At AoA = 12°, the airfoil is already in a pre-stall condition with local flow separation occurring at the TE, however the flow is still primarily attached to most of the airfoil surfaces, as illustrated by the mean velocity contours (e.g., see Figure 10). The boundary layer is characterized by the presence of a small laminar separation bubble (LSB) near the leading edge (LE). The signature of this LSB is shown in Figure 11 in both the pressure coefficient (C_p) and the skin friction coefficient (C_f), in addition to a close-up view of the mean velocity magnitude. The laminar boundary layer (LBL) instabilities are also captured just upstream of the LSB, as shown with the unsteadiness in the C_f curve.

The LSB triggers the transition from a laminar to a turbulent BL at x/c = 0.0123, where C_f ~ 0. This LSB also causes a separation-induced transition, resulting in a turbulent boundary layer developing downstream of the LSB, starting at x/c = 0.04. This behaviour near the transition was captured through the use of the intermittency transition model; similar behaviour of both C_p and C_f curves has been associated in the literature with LSB triggering the onset transition when captured by transitional models [39].

The transition could also be captured by the intermittency contours (Figure 11), where a value of 0 indicates a laminar BL and a value of 1 indicates a fully developed turbulent BL. On the suction side, the flow is completely laminar. At the TE, the turbulent boundary layer eventually separates due to the larger adverse pressure gradients there, as illustrated in the C_f curve in Figure 11.

Figure 12 and Figure 13 show the vorticity contours at AoA = 12° for two different frequencies: 2 and 8 Hz. Overall, the cases do not show striking differences; perhaps the only noticeable difference is the size and vorticity of the initial vortex forming just before the morphing motion stops. Nevertheless, from the qualitative inspection of the flow responses to the morphing deflection, further insights into the mechanism by which the dynamic stall occurs can be gained.

As soon as the morphing motion starts, the separated turbulent boundary layer flow at the TE begins to propagate upstream, and in contrast to AoA = 10°, the propagation accelerates on the airfoil’s suction surface with increasing camber due to downward flap deflection. The unsteady dynamic stall flow appears to be triggered when the large-scale turbulent flow structures interact with the LSB present near the LE, which causes the latter to burst, forming a large-scale vortex which then convects downstream. In the meantime, it induces a small-scale leading-edge vortex (LEV). The turbulent boundary layer becomes separated at the TE and begins to form a second vortex that also convects downstream and merges with a counter-clockwise TE vortex. Thus, it seems that the dynamic stall is associated with the bursting of the laminar separation bubble (LSB) present in the leading-edge region.

A remarkable similarity could be found in the flow responses observed in this study to those reported by Benton et al. [40] for a ramp-type pitching-up motion at similar AoA. Benton et al. [40] studied an NACA 0012 airfoil experiencing a rapid pitching-up motion at Re = 1,000,000 using a high-fidelity LES. They reached similar conclusions regarding the mechanism that causes the airfoil to enter deep stall. This indicates that the mechanism governing the unsteady behaviour of a morphing flap has similarities with ramp-type pitching-up of a ‘rigid’ airfoil, which would possibly allow future studies of morphing downward flaps to be approximated with a rapid pitching-up motion for simplicity and with low computational overhead costs.

7. Conclusions and Future Work

High Reynolds number transient flow responses to various deflections of a morphing TEF were investigated using a hybrid RANS-LES turbulence model at three angles of attack near stall (8°, 10°, and 12°). Using a pre-defined large-amplitude flap deflection, the effects of various morphing frequencies on the dynamic response of the aerodynamic forces and on the surrounding flow field were explored.

For all AoA, both the lift and the drag coefficient responses showed a proportional relationship between the morphing frequency and the slope of response at which these quantities evolve. Interestingly, an overshoot in the drag coefficient was captured, even at quasi-static conditions, but this was not seen in the lift coefficient. The amplitude of spike in the drag coefficient was also found to be proportional to the morphing frequency, with higher frequencies causing larger overshoots. The results obtained for the range of frequencies investigated at a Reynolds number of 0.62 $\times {10}^6$ show the potential of rapid downward deflections for obtaining instantaneous lift increments. However, the frequencies studied did not seem to produce a recovery from the stall, and instead at the higher AoA of 12° the morphing even drove the airfoil into a dynamic stall.

Qualitative results have given some insights into the mechanism of flow responses to a rapid morphing deflection. At all AoA, the downward morphing motion enhances vortex shedding, although the mechanisms are different at lower and higher AoA. In particular, at an AoA = 12°, the morphing causes the turbulent flow structures that form at the TE to propagate upstream before interacting with the leading-edge laminar separation bubble and causing it to burst, forming a large-scale vortex that convects downstream and initiating a dynamic stall. This behaviour was also observed in the literature in a rapid pitching-up airfoil LES study.

The findings of this study indicate that at the high Reynolds number investigated, a downward morphing TEF can be used to obtain instantaneous lift increments for operational purposes, however care must be exercised to avoid large-scale flow instabilities and deep stall, especially at AoA much closer to stall.

For future work, wall-resolved 3D high-fidelity simulation is needed (e.g., using LES) in order to resolve small-scale turbulent flow structures, which would enrich our understanding of how the flow responds to dynamic morphing flap deflection. Additionally, near-stall AoA should be investigated at higher frequencies (e.g., a deflection time of an order of one convective time) to gain further understanding of unsteady flow physics, to investigate possible stall mitigation, and to explore the capabilities of TE morphing for load alleviation purposes. Lastly, it is important to extend the scope of the current study into different morphing motions, different frequency–magnitude combinations, and a wider range of Reynolds numbers.