2.1 Meshes

The engineer needs to ensure that the computational mesh is sufficient to capture the phenomena that drive and arise from the anticipated separation physics, as well as a timestep sufficient to capture the dynamic motion. The first and most important aspect of meshing is understanding that separation is a three-dimensional process. Two-dimensional Navier-Stokes simulations were performed by necessity in the 1980s due to computational limitations. It was apparent that, while a single plane would provide a general qualitative assessment, it could not accurately capture all of the salient details of separated flow, even without dynamic motion. Unfortunately, even today some users persist in attempting to predict static and dynamic separation with a single plane. They perhaps erroneously assume that physical aerofoil experiments describe an aerofoil undergoing separation in two dimensions, despite the fact that the experimental data were informed by a three-dimensional flow.

Computational fluid dynamic solvers, no matter their formulation, rely on refined meshes in the regions where the flow is rapidly changing, including areas where laminar-to-turbulent transition may occur. Mesh independence studies are routine practices in computational aerodynamics to ensure the mesh is sufficiently accurate to capture the physical phenomena of interest. Since these studies are computationally expensive, most focus on the streamwise – normal planes (e.g., planes at constant wing span or blade radius locations). Exemplar studies for aerofoil separated flows include Szydlowski and Costes [Reference Szydlowski and Costes28] and Smith et al. [Reference Smith, Liggett and Koukol29]. Best practices, in addition to the well-known wall spacing ($y + \le 1$), recommend that the mesh must include sufficient points orthogonal to the surface (20 – 60) to capture the onset of separation within the smooth surface boundary layer, and the normal cell size growth rate should be limited to 10% [Reference Szydlowski and Costes28, Reference Liggett and Smith30, Reference Liggett and Smith31]. The span or radial surface mesh requirements are typically overlooked, but meshes that expand too quickly from the refined edge regions or large aspect ratios can cause problems in capturing separation. Kolpitcke and Smith [Reference Grubb and Smith32] provide further best practices guidance for both translating and rotating lifting surface meshes. In their study, they quantified the role that radial mesh refinement, including chord-to-span mesh aspect ratio, has on the accurate prediction of these separated flows.

Meshes can be constructed of structured hexahedron elements assembled from two-dimensional quadrilateral sides. Unstructured meshes can be composed of multiple element shapes, and a tetrahedon is the most common element applied as it can easily conform to corners and other discontinuities in the geometry. Prisms, which replicate a length-diagonal cut of the hexahedron, provide a more accurate representation of smooth surface boundary layers. Mixed-element or hybrid meshes currently provide the most optimal meshing solution for very complex moving geometries. Some computational frameworks, such as CREATE${{\rm{\;}}^{{\rm{TM}}}}$-AV Helios, enable the use of different computational solvers (structured, unstructured, Cartesian) through its multi-mesh, multi-solver paradigm to provide an optimal mesh solution.

Unstructured meshes can simplify the meshing process of complex geometries, while overset meshes can solve the requirement of motion in different frames of reference, as illustrated by the helicopter in Figure 5. The rotating rotor and hub frame is overset with the translating fuselage and off-body frame. The rotor blade mesh in the upper right illustrates an unstructured mixed-element mesh where smooth body separation may occur, applying prisms near the surface and tetrahedrons beyond that region. The complex fuselage geometry is meshed with tetrahedrons that provide a more cost-efficient mesh where bluff-body separation will be present. The farfield or off-body volume domain is composed of a series of overset Cartesian meshes to capture the wake features without the additional computational overhead of body-fitted mesh mathematics. The Cartesian meshes are layered with increasing mesh volume size to ensure that flow features in the near wake are not overly dissipated by the numerical scheme. Example sizes of volume meshes are included. Adaptive mesh refinement (AMR) can be applied to refine the off-body to more accurately capture these wake features. However, AMR can add significant costs to the solution [Reference Jain33], while fixed refinement regions may be more cost-effective to resolve important flowfield features in some simulations.

Figure 5. Illustration of a modern complex unstructured, overset mixed-element Helios mesh about a UH-60A rotorcraft. Modified from Figure 7 in Roget et al. [Reference Roget, Sitaraman, Wissink, Saberi and Chen34].

While unstructured meshes have the ability to rapidly mesh complex bodies, they also come with increased computational memory and time requirements compared to structured meshes of the same size. Structured meshes, which may require many hours to develop complex meshes and usually involve oversetting or piecing together contiguous meshes, are most computationally cost-efficient. In addition, structured mesh solvers can resolve higher-order spatial schemes (fourth, fifth and sixth order are common), while current unstructured mesh solvers are limited to nominally second-order spatial schemes. Therefore, structured solvers that apply higher-order spatial schemes can achieve the same accuracy as second-order schemes with fewer mesh cells, or these higher-order structured solvers can provide more accurate solutions with the same mesh size.

During validation with experimental model papers, the influence of the wind tunnel walls and model support structures must be assessed. The additional meshing requirements of these features can easily require a 150%–200% increase in the volume mesh cell count, depending on the complexity of the structures. In addition, as many of these validation cases involve subscale models, it may be important to consider the presence of laminar to turbulent transition. The mesh may need to be further refined in the region of transition (Section 2.3), and the timestep adjusted to accommodate these requirements that may not be present in full-scale applications.

While mesh/solver type is highly dependent on the use case, outcome required, and what solvers are available, best practices should be applied – no matter the application. In particular, novice users should perform literature surveys, to include some of the references listed herein, to determine the best practices for their use case and solver. Mesh and timestep independence studies for new applications are also normal practices by diligent engineers.

2.2 Turbulence closures

Accurate turbulence closures hold the key to successful solutions of separated flows, and they are the primary recipient of the blame when computational solutions do not meet expectations of accuracy, in particular when compared to experimental data. Since the 1970s, researchers have sought more accurate turbulence closures, primarily through the development of models that reproduce the macroscopic effects of turbulence without the need to resolve all spatio-temporal scales. A synopsis of the different types of turbulence closures for the resolution of the Navier-Stokes equations, and their ability to capture or model turbulence scales is illustrated in Figure 6.

Figure 6. Types of turbulence closures with their ability to capture (blue) or model (red) physical scales of turbulence.

Reynolds-averaged Navier-Stokes (RANS) or uRANS solvers model all scales of turbulence. Early algebraic models have given way to one-equation partial differential (PDE) models, which in turn have been replaced by two-equation PDE models for modeling flows with separation. These models have a number of tuning variables, which are typically constants, based on validation problems similar to the character of the problem they are trying to resolve. Various corrections for rotation and limiters on dissipation are common.

A recent study by Sridhar et al. evaluated the influence of two popular versions (Menter and Kok) of the two-equation $k - \omega $ turbulence model [Reference Sridhar, Crawford and Smith35]. They found that while the two models predicted comparable results for attached flows for configurations that were both moving and static, once separation was introduced, the methods provided similar but different results with all other numerical considerations kept identical. Differences were due to the numerical growth of the turbulent eddy viscosity in the Kok version which also led to numerical instability in some cases. Limiters in the shear stress transport (SST) algorithm added to the Kok model capped the growth, bringing results in line with other models and experiment.

At the other extreme, direct numerical simulations (DNS) attempt to resolve all scales of turbulence within the flowfield. DNS remains a research tool applied to further extend the understanding of the physics within the complex turbulence process. For example, Yeung and Ravikumar [Reference Yeung and Ravikumar36] further analysed the role of intermittency within turbulence (at low Reynolds numbers) using DNS, which required a mesh of over six trillion mesh points.

Large-eddy simulation (LES) provide a compromise between RANS and DNS. LES captures a range of the larger scales of turbulence, based on the mesh size, and model the smaller, less energetic scales. LES simulations with less refined meshes that resolve a smaller subset of the larger turbulence scales is also known as very large LES (VLES). The use of LES can be prohibitively expensive to resolve the very small Kolmogorov scales of turbulence within the boundary layer at larger Reynolds numbers. These very small scales drive the required mesh refinement and timestep size reduction that limit the use of LES for most practical aerospace engineering applications.

To resolve flows where the primary goal is to capture complex features such as separation, researchers have developed innovative solutions – wall-modeled LES (WMLES) – to address the costs of wall-resolved LES, while still retaining much of the accuracy that LES provides. Figure 7 is descriptive of a general definition of these turbulence closures. The characterisation in this figure focuses on the scales of the physical phenomena rather than the turbulence scales which are inherent within the boundary layer.

Figure 7. Types of turbulence closures for practical application to separated flows at higher Reynolds numbers.

WMLES has taken on many forms in the past two decades of its primary development. As the boundary layers contain the smallest scales of turbulence, alternatives to LES modeling are applied. Two primary focii have emerged in this endeavour: hybridisation of LES with uRANS and development of wall stress models for LES, similar to those applied for uRANS. For both approaches, the smaller energetic scales of the inner layer (nearest the wall) boundary layer are modeled, using the outer layer LES solution as a boundary condition. Larsson et al. provide a detailed overview of WMLES [Reference Larsson, Kawai, Bodart and Bermejo-Moreno37].

In the first approach, the boundary layer – where the smallest, most prohibitive scales of turbulence reside – is modeled with uRANS, while the separated regions are predicted with LES. These hybrid uRANS-LES methods resolve either formal sub-grid scale differential equations for length and turbulent kinetic energy, or the family of detached eddy simulation (DES) approximations [Reference Spalart, Jou, Strelets and Allmaras38]. These hybrid schemes may apply either simplified algebraic relations of the length scale or the length scale and eddy viscosity as the trigger between the two methods. Typically, these hybrid approaches evaluate and weight the uRANS or LES/DES option. This has given rise to a number of different switching options between the inner and outer boundary layer zones. Examples of these options include delayed detached eddy simulation (DDES) [Reference Spalart, Deck, Shur, Squires, Strelets and Travin39], zonal DES [Reference Renard and Deck40] and hybrid RANS-LES (HRLES) [Reference Sánchez-Rocha and Menon41, Reference Sánchez-Rocha and Menon42]. The hybrid uRANS-LES approach has been more widely adopted for use in analyses that involve separation with moving configurations. In particular, the rotorcraft community has adopted this approach as it permits the extension of extant uRANS solutions that include critical capabilities such as multiple reference frames of motion, rotor trim and aeroelastic rotor blades and the CREATE${{\rm{\;}}^{{\rm{TM}}}}$-AV Helios framework.

The wall-stress modeling approach is similar to the uRANS-LES approach, varying in the solution of the inner layer. The uRANS models typically solve the one- or two-equation differential equations (e.g., Spalart-Allmaras or $k - \omega $ models), while the wall-stress approach assumes equilibrium at the edge of the boundary layer and solves an algebraic form of the eddy-viscosity model.

Table 1. Predicted characteristics for various turbulence methods and grids for a semi-infinite circular cylinder at ${\textit{Re}}$ = 3,900. Separation location is given in degrees of azimuth from the leading edge stagnation point. Modified from Lynch and Smith [Reference Lynch and Smith26].

Figure 8. Isocontours of Q-Criterion for a semi-infinite cylinder at Re=3,900 at the same mesh (101 span stations) and timestep, from Lynch and Smith [Reference Lynch and Smith26].

A question or criticism that occasionally arises with the use of these WMLES methods regards the re-use of RANS meshing and/or timestep in the WMLES simulations. The difference between LES that resolves the boundary layer and these wall-modeling LES approaches is that the smallest scales of turbulence that constrain mesh cell size and timestep are avoided with WMLES. Instead, the mesh in the wake, and by extension the timestep, should be appropriate to the larger scales that define the separated flow – usually the size of wake vorticity. Thus the traditional highly refined mesh and timestep sizes of boundary layer LES can be avoided, and best practices associated with the meshes and time-steps for these scales (and uRANS) can be applied. The validity of these assumptions has been demonstrated in a large number of applications, as well as for some canonical LES and uRANS cases.

Lynch and Smith confirmed these assumptions for a semi-infinite cylinder for a wake transitional Reynolds number [Reference Lynch and Smith26], which is a standard LES test case. Although the cylinder is stationary, its complex wake provides some insight into the development of these methods. The two-equation uRANS model captures the primary shed wake or ‘rollers’ but without any wake artifacts due to the instabilities of the spanwise flow, a phenomenon captured by the WMLES, as illustrated in Figure 8. In addition to the lack of wake refinement, the two-equation uRANS model overpredicts the drag by 50% and the separation location by more than 10 degrees, elucidated in Table 1. Their hybrid uRANS-LES approach was then evaluated for different span extents to examine the spanwise mesh requirements, where the results are within experimental error bounds and comparable to the wall-resolved LES computations of Kravchenko and Moin [Reference Kravchenko and Moin43]. The LES wake region acts as a boundary condition for the near-wall uRANS region, which encompassed the cylinder outside the wake and included the attached boundary layer. The LES wake influences and improves the uRANS model, similar to early computational hybrid approaches, such as those that coupled the boundary layer equations with the full potential equations. The centreline pressure distribution around the cylinder confirms this hypothesis in Figure 9. Beyond 60 degrees from the front centreline of the cylinder, the uRANS model underpredicted pressure coefficient (creating a nonphysical region of suction) so that the aft pressure distribution is significantly different from the experiment and the higher-fidelity methods, with a non-physical suction pressure region on the aft portion of the cylinder (Figure 10).

Figure 9. Mean pressure coefficient for the semi-infinite circular cylinder at ${\textit{Re}}$ = 3,900. Figure 3 from Lynch and Smith [Reference Lynch and Smith26], including additional data from Krachenko [Reference Kravchenko and Moin43].

Figure 10. Mean pressure coefficient for a semi-infinite NACA0015 wing at 90${{\rm{\;}}^ \circ }$ angle-of-attack and ${\textit{Re}}$ = 1,000,000. From Smith et al. [Reference Smith, Liggett and Koukol29].

Smith et al. [Reference Smith, Liggett and Koukol29] also confirmed this behaviour for aerofoils at high angles of attack and supercritical Reynolds numbers of one million. The leeward side of the aerofoil, modeled as a semi-infinite wing, exhibits similar underpressure variation with uRANS as the cylinder, while the hybrid RANS-LES approach is able to correctly predict the pressure variation in the separated region, which is key before adding the complexity of unsteady motion.

Hodara et al. extended this evaluation to moving bodies in 2016 [Reference Hodara, Lind, Jones and Smith46], applying a different computational solver to a wing in reverse flow. For rotating systems, the addition of forward flight for rotorcraft or yaw for sustainable energy turbines will yield a situation where the local freestream flow starts at the trailing edge and moves towards the leading edge. The artifacts in the shed flowfield can be captured accurately as they emanate from the sharp trailing edge, making this an excellent case study for computations of static and dynamic bodies with separation. Once again, uRANS was insufficient to capture the separation and the shed wake. DDES, which tend to be more sensitive to the mesh and timestep, accurately captured the integrated loads. The details of the wake phenomena were best captured when a full subgrid-scale differential equation (here, the turbulent kinetic energy or $k$ equation) was resolved. When modeling these unsteady oscillations, some computational vortices experienced a phase-locked lag in release compared to experiment. These lags occurred only during portions of the oscillation cycle when the flow was separated, illustrated in Figure 11. However, all flowfield features were captured and present when the WMLES turbulence closures were applied. The delay in release of the vortices appeared to be related to the underlying uRANS model and were exacerbated with larger timesteps.

Figure 11. Comparison of the experimental and WMLES-predicted flowfield of a wing undergoing oscillations with separation in a reverse flow [Reference Hodara, Lind, Jones and Smith46].

There has also been additional development with traditional LES to make it more tractable for research applications. Implicit scheme LES (ILES) applies larger timesteps than the tradition explicit schemes used by LES. The application of these ILES approaches still have limitations when transition or separation is present, as they may limit the size of the timesteps to maintain stability. The addition of dynamic motion of the vehicle or flowfield can further limit the timestep size. While this approach requires much smaller timesteps than WMLES closures, the introduction of larger CPU or GPU clusters can partially mitigate this time penalty. Using ILES, Visbal [Reference Visbal47] was able to evaluate a semi-infinite wing undergoing plunge oscillations at a series of Reynolds numbers where the transition occurred near the leading edge (10,000; 40,000; 60,000). The analysis, which was correlated with experimental data, indicated the presence of complex flowfield primary and secondary structures (Figure 12(a)) similar to those obtained with WMLES (e.g., Figure 8).

Figure 12. ILES evaluation of a plunging SD7003 aerofoil at a chord-based reduced frequency of ${{\rm{k}}_{\rm{f}}}$ = 3.93 at a ${\textit{Re}}$ = 40,000. ${\rm{\Phi }}$ denotes the location of maximum displacement. Figures modified from Visbal [Reference Visbal47].

A major advantage of the traditional LES (WRLES) approach, in this case ILES, is that these methods are able to accurately capture the dynamic onset of the transition between laminar and turbulent flows, as well as small areas of separation and reattachment, as exemplified in Figure 12(b). However, engineering computations with relevant full-scale Reynolds numbers are still out of reach for most applications.

2.3 Transition models

While traditional LES is able to naturally predict the onset of transition with sufficient mesh resolution in the boundary layer, a reliable transition model for uRANS and WMLES is necessary at higher Reynolds numbers (100,000 and higher) to avoid the high costs associated with smaller turbulence scales as Reynolds numbers increase. The ability to predict transition is arguably more important for computational validations with subscale experiments than for many full-scale aerospace applications. Transition modeling, like turbulence models, are often held accountable for otherwise unidentifiable discrepancies between experimental and computational predictions. There are a number of transition models or techniques that can be coupled with uRANS turbulence models and applied in engineering simulations. Two techniques have been applied to uRANS and uRANS-LES approaches with mixed success: the turbulence intensity/momentum thickness ($\gamma R{e_\theta }$) Langtry-Menter model (LM) and the amplification factor (AFT).

As part of the AIAA Transition Workshop, the LM transition model was applied to several common static test cases [Reference Carnes and Coder48]. The conclusions noted that while overall the model predicted transition reasonably well, it suffered from freestream turbulence decay, which in turn delayed the onset of transition. In addition, mixed results were obtained with varying grid refinement levels, which appeared to be in part related to the type of transition. Flat plate and strongly attached flows had less success than flows with incipient separation. This latter observation suggests that the LM model could be accurate in flows with more separation and when transition is less influenced by the freestream turbulence.

The AFT model, first developed in 2014, also has many features that are well suited for high performance computing. It is rooted in the local turbulence intensity and the local boundary layer shape factor. Coder earlier assessed the AFT approach using the Spalart-Allmaras (SA) model with the same solver, although most test cases were not the same as the LM assessment [Reference Coder49]. The SA-AFT simulations did not appear to be influenced by freestream turbulence decay, as the onset of transition for these cases was noted to be slightly early compared to experimental data, while the LM model was late. Predictions as angles of attack increased followed this trend. The model at the time of the analysis did not include crossflow corrections, so that transition in flowfields that have strong crossflow components, such as rotating blades and separated wing flows, could less reliably be assessed.

In addition to the ability to correctly predict crossflow in separation, the ability to maintain the Galilean invariance (GI), or proper kinematics across multiple frames of motion, is critical. To date, one of the best efforts to demonstrate the importance of cross-flow and GI is a study of rigid rotors by Jain [Reference Jain33]. He computationally evaluated two rotors with both the LM and AFT models and correlated the results with experimental transition locations identified using differential infrared thermography (DIT). For both rotors, good agreement is found between the computed and measured transition locations if the transition model includes the GI terms. The influence of cross-flow is more nuanced. The two-dimensional inboard flows at higher blade pitch angles are dominated by the presence of the laminar separation bubbles (LSB), and their prediction is not influenced by the inclusion of the crossflow terms. However, outboard on the rotor, where the finite tip region is dominated by radial flow rather than LSBs, the crossflow terms are necessary to accurately capture the transition locations during a rotor revolution (${0^ \circ } \le \psi \le {360^ \circ }$). An example of the study results is illustrated in Figure 13 for the pressure-sensitive-paint (PSP) rotor at the $r/R = $ 0.85 radius location. The results from the AFT transition model coupled with the Spalart-Allmaras uRANS model (SA-AFT) are typically not as close with experimental results as the results from the Galilean invariant LM model coupled with the $k - \omega $ SST uRANS model (SSTM-LM-G) when compared to experiment (DIT). The SA-AFT also exhibits a sharp peak that needs further exploration to understand its underlying cause. The SST-LM-G predicts a transition location that is farther aft and delayed (later $\psi $ location) than its counterpart with crossflow (SST-LM-G-CF). The additional of crossflow also resulted in larger regions of turbulent flow in the aft rotor disk (${180^ \circ } \le \psi \le {360^ \circ }$) where the rotor blades can experience separated flows due to interactional aerodynamics.

Figure 13. Comparison of the transition predictions with experiment at a radial location of ${\rm{r}}/{\rm{R}}$ = 0.85 for the PSP rotor at ${\rm{\mu }}$ = 0.3 and ${{\rm{\alpha }}_{\rm{s}}}$ = −3${{\rm{\;}}^ \circ }$, from Jain [Reference Jain33].

Overall, these transition models, when applied with best practices, yield predictions that are close to those observed in physical experiments. Two important numerical considerations must be addressed when transition is present. First, the mesh needs to be sufficiently fine on the surface to capture the transition and LSB, if it is present. This usually requires a finer surface mesh than fully turbulent flow assumptions. In addition, the sensitivity of the solution to the numerical scheme types and limiters is also much greater than other computations.

3.1 Dynamic stall

No discussion on the topic of unsteady separated flows would be complete without discussion of dynamic stall. Classic dynamic stall is characterised by an increasing angle-of-attack over time until stall occurs, as illustrated in Figure 14. The stall persists as angle-of-attack increases and only recovers when the angle-of-attack is low enough so that the flow can reattach, as is illustrated by the different steps. This definition of dynamic stall is typically attributed to the Carr, McAlister and McCroskey [Reference Carr, McAlister and McCroskey50]. The angle-of-attack forcing function is oscillatory and can be characterised by a nondimensionalised reduced frequency, ${k_f} = \frac{{{\omega _f}b}}{V}$. Care must be taken with the ${k_f}$ parameter, as some aerospace sectors define it with the wing or aerofoil chord, $c$, rather than the traditional semi-chord ($b = \frac{c}{2}$).

Figure 14. Pictorial definition and flowfields illustrations of classic dynamic stall. Courtesy N. Liggett.

As the efforts to accurately predict dynamic stall computationally span the past 60–70 years (Figure 15), this phenomenon warrants a discussion here. Dynamic stall occurs in fixed-wing flight conditions, for example, fighter aircraft abrupt stall and aeroelastic phenomena such as control surface buzz. It is also prevalent in rotating systems, such as rotorcraft vehicles and wind energy turbines.

Figure 15. Timeline of dynamic stall computational development, highlighting advances achieved during the past 25 years.

Prior to the year 2000, efforts were primarily restricted to two-dimensional applications due to limits in computational hardware. As denoted in Figure 15, major advances parallel the increase in hardware and software capabilities. Progress generally, but not always, followed the characterisation of dynamic stall in two dimensions from aerofoils to semi-infinite wings, to finite wings, and then to full vehicles. Dynamic stall is a critical phenomenon in rotorcraft, impacting both vibratory loading and performance, so significant resources have been expended in rotorcraft research to understand, predict, and mitigate dynamic stall.

Computational characterisation of dynamic stall moved from simple aerofoils to semi-infinite wings as computational resources improved. Early efforts focused on improving the accuracy of computationally generated C-81 lookup tables. These tables are applied in blade-element methods that comprise the bulk of the rotorcraft engineering analyses [Reference Smith, Wong, Potsdam, Baeder and Phanse51]. The modeling of three-dimensional semi-infinite computational domains was an important advancement since turbulence and separation are three-dimensional phenomena. Computational validations had previously been performed with experimental campaigns that evaluated simple geometries (wall-to-wall rectangular wings of a single aerofoil). Those correlations indicated that the aerodynamic coefficients of the experiments exhibited significant three-dimensional relief effects that were not present in the uRANS-predicted aerodynamic coefficients. When semi-infinite computations were performed, the aerodynamic coefficients more accurately represented the mean cyclic behaviour of the experiments. Both two- and three-dimensional studies identified best practices for the timestep and mesh requirements for dynamic stall [Reference Smith, Liggett and Koukol29, Reference Szydlowski and Costes52].

Concurrent with and following these studies, finite wings undergoing dynamic stall were also being computationally evaluated. Spentzos et al. [Reference Spentzos, Barakos, Badcock, Richards, Wernert, Schreck and Raffel53, Reference Spentzos, Barakos, Badcock, Richards, Coton, Galbraith, Berton and Favier54] demonstrated the ability of uRANS (two-equation $k - \omega $) computations to capture the arch or omega vortex detected in dynamic stall experiments from numerous sources, including Moir and Coton a decade earlier [Reference Moir and Coton55], as illustrated in Figure 16. They also characterised the presence of large radial flows when separation appeared. The latter observation was important as rotorcraft researchers could first study the ability of computational methods to capture radial flows during fixed-wing dynamic stall without adding more complex and costly rotation.

Figure 16. Correlation of computational dynamic stall (right) with experimental smoke visualisation (left) for an aspect ratio 3 NACA0015 wing at a ${\textit{Re}}$ = 13,000, ${\rm{M}}$ = 0.1 and a nondimensional pitch rate (${\rm{\alpha }} + = \frac{{{\dot \alpha \rm {c}}}}{{{{\rm{U}}_\infty }}}$) of 0.16. Planform view across the top; leading edge looking aft across the bottom. From Spentzos et al. [Reference Spentzos, Barakos, Badcock, Richards, Coton, Galbraith, Berton and Favier54] with experimental data originally from Moir and Coton [Reference Moir and Coton55].

These early efforts were followed by further evaluations on oscillating finite wings using uRANS [Reference Kaufmann, Costes, Richez, Gardner and Le Pape56–Reference Kaufmann, Merz and Gardner58], WMLES [Reference Richez, Le Pape and Costes57, Reference Jain, Le Pape, Grubb, Costes, Richez and Smith59] and LES [Reference Visbal and Garmann60]. A series of evaluations by the rotorcraft community were correlated with experimental data for the OA205 finite wing with an aspect ratio of 3 and a Reynolds number of 1 million for static and dynamic stall [Reference Le Pape, Pailhas, David and Deluc61]. Kaufmann et al. confirmed the ability for uRANS to capture the omega vortex [Reference Kaufmann, Costes, Richez, Gardner and Le Pape56]. Richez et al. [Reference Richez, Le Pape and Costes57] found that uRANS was not able to capture many of the important underlying physics, but their WMLES was able to capture them. Jain and his collaborators evaluated the OA205 finite wing with different solvers and uRANS and WMLES (DDES) for both static and dynamic stall [Reference Jain, Le Pape, Grubb, Costes, Richez and Smith59]. They confirmed the earlier observations of radial flows in dynamic separation, and, at the Reynolds number evaluated ($Re = 1 \times {10^6}$), that transition models improved the airload accuracy prediction when correlated to experiment. They evaluated both light (small zones of separated flow) and deep (large areas of separated flow) dynamic stall, indicating that when the best practices adopted from the semi-infinite two-dimensional computational studies were applied, the computational approaches predicted both types of dynamic stall well.

Air Force Research Lab (US) researchers studied finite wing dynamic stall at lower Reynolds numbers ($Re = 200,000$) with their ILES solvers [Reference Visbal and Garmann60]. In this work, they were able to computationally characterise the influence of aspect ratio and wall mounting. The character of the arch or omega vortex was investigated for various scenarios, confirming the earlier findings with uRANS and WMLES. Although computed at lower Reynolds numbers and currently too costly for engineering analyses, the LES analyses are exceptionally valuable in further elucidating the details and importance of the physical phenomena during the dynamic stall. Some phenomena drive accurate airloads, while others may be less important. An example is the development of the omega or arch vortex, depicted in Figure 17. While all approaches have captured the macroscale vortex, the developmental details and strength of the vortex that require much smaller scales can be studied with these LES simulations. The importance of the developmental features then drive the mesh, timestep, transition and turbulence details of the approaches with additional modeling assumptions.

Figure 17. Illustration of the omega vortex interactions during dynamic stall on a pitching NACA0012 aspect ratio 4 wing evaluated with ILES for ${{M}_\infty }$ = 0.1 and ${{R}}{{{e}}_{{c}}}$ = 200,000. Modified from Visbal et al. [Reference Visbal and Garmann60].

Another application that has relevance to dynamic stall and aeroelasticity is the motion of control surfaces with respect to a wing that is either static or moving itself. In this exemplar study, the use of a WMLES turbulence closure ($k - \omega$ uRANS with $k$-equation LES subgrid-scale model), combined with overset meshes in relative motion, illustrates the importance of modeling details such as the gaps between the flaps [Reference Liggett and Smith30, Reference Liggett and Smith31]. Computational evaluations were performed with both a separate flap that included the gap and an integrated flap that did not include the gap. In this study, a NACA0012 wing at ${M_\infty }$ = 0.4, $Re$ = 1.6M and $\alpha $ = 4${{\rm{\;}}^ \circ }$ included a flap oscillating at a reduced frequency, ${k_f}$ = 0.042 while the aerofoil oscillated at ${k_f}$ = 0.021. The flowfield variation during one cycle of the flap oscillation (Figure 18(a)) illustrates some of the complex interactions that occur within the gap. The influence of modeling these gap interactions is clearly observed in integrated aerodynamic variables. Lift coefficient is the exemplar variation in Figure 18(b), though similar trends are observed with drag and pitching moment. While the overall trends are captured with an integrated flap, the predicted magnitudes of the coefficients can be larger than the modeled gap results and experiment [Reference Krzysiak and Narkiewicz62]. The gap provides a relief effect that, with the WMLES turbulence closure accurately follows the experimental results.

Figure 18. Computational simulations of an oscillating NACA0012 aerofoil-flap semi-infinite wing at ${{\rm{M}}_\infty }$ = 0.4 and ${\rm{R}}{{\rm{e}}_\infty }$ = 1.63 million. The main aerofoil oscillates about ${{\rm{\alpha }}_{{\rm{mean}}}}$ = 4${{\rm{\;}}^ \circ } \pm {6^ \circ }$ at ${{\rm{k}}_{\rm{f}}}$ = 0.021. The flap oscillates about ${{\rm{\alpha }}_{{\rm{mean}}}}$ = 0${{\rm{\;}}^ \circ } \pm {6^ \circ }$ at ${{\rm{k}}_{\rm{f}}}$ = 0.042. From Liggett et al. [Reference Liggett and Smith30].

Advances in the accuracy of high-quality computational data over the past two decades in particular have spurred improved and concurrent collaboration with experimentalists. For example, the traditional experimental practice provided a single point per timestep for an analysis, typically time-averaged, or if periodic, phase averaged is no longer sufficient. In some instances, error bars associated with the instrumentation uncertainty have been provided. More recent experimental campaigns [Reference Gardner, Richter, Mai, Altmikus, Klein and Rohardt63–Reference Gardner, Jones, Mulleners, Naughton and Smith65] provide either mean data with standard deviation bounds or unsteady min/max extents of the measured variable over many oscillations to provide an estimate of ‘what is good enough’ to for computational validation.

Recognition of the unsteadiness during experimental campaigns associated with these unsteady flowfields is not a new observation. Carr, McAlister and McCroskey observed the presence of flow variations in their now classic 1977 NASA Technical Note [Reference Carr, McAlister and McCroskey50]. Experimentalists deduced these effects from ‘random turbulent fluctuations’ or other sources by phase averaging over many cycles.

In some instances, this phase averaging can lead to erroneous or biased data. One example of this is associated with dynamic stall. Computational results of integrated aerodynamic coefficients consistently indicated a stronger nonlinear lift peak generated by the leading edge vortex just prior to stall over their experimental counterpart. These differences could not be explained through mesh and time-step refinements, three-dimensional effects, or advanced turbulence models. Ramasamy et al. [Reference Ramasamy, Wilson, McCroskey and Martin66] undertook a more detailed analysis of recent experiment data and established that phase averaging of the data smoothed a slight aperiodicity in lift peak, introducing a bias in the data and reducing the magnitude of the nonlinear lift peak. Recent efforts by Ramasamy et al. [Reference Ramasamy, Sanayri, Wilson, Martin, Harms, Nikoueeyan and Naughton67] identified using cycle-to-cycle experimental data that dynamic stall can include one or more furcations associated with different physics present in a cycle. These different modalities may be caused by a number of physical phenomena, including but not limited to variations in the separation location, leading edge or trailing edge stall onset, presence of dynamic stall vortices, and time of reattachment angle. Test conditions (wall conditions, flow aperiodicity, random turbulence fluctuations) may also induce cycle-to-cycle variations, where phase averaging can introduce bias or erroneous conclusions.

Tran et al. [Reference Tran, Sitaraman and Ramasamy68] further explored this behaviour using both experimental and computational results. Figure 19(a) elucidates the bimodal behaviour with the previous experimental data from Ramasamy et al. [Reference Ramasamy, Sanayri, Wilson, Martin, Harms, Nikoueeyan and Naughton67]. Cluster 1 in blue indicates cycles with a dynamic stall, while cluster 2 cycles in red do not have a stall event. As cluster 2 dominates the experiment for 87% of the time, phase averaging of all of the data cycles (black line) indicates that dynamic stall is not present, leading to an erroneous conclusion. The authors also performed a computational analysis with Helios, over many cycles with refined temporal analysis and WMLES-based turbulence models (DDES). These computations reproduced these bifurcations with excellent accuracy, though the overall combined analysis resulted in some differences as the computations reported different weighting of each cluster compared to its experimental counterpart. This collaboration was key as the computational results provided additional flowfield insights to aid in the explanation of the phenomena driving these variations (Figure 19(b)). Detailed collaborative analyses such as these are providing new insights into these complex separated flows.

Figure 19. Variational cycle-to-cycle cluster comparison in dynamic stall for a modified VR-12 aerofoil at ${\rm{\alpha }} = {18^ \circ } \pm {5^ \circ }$, ${{\rm{M}}_\infty }$ = 0.3, ${{\rm{k}}_{\rm{f}}}$ = 0.10, from Tran et al. [Reference Tran, Sitaraman and Ramasamy68] and experimental data from Ramasamy et al. [Reference Ramasamy, Sanayri, Wilson, Martin, Harms, Nikoueeyan and Naughton67].

The costs of the simulations discussed thus far are very prohibitive, so active research is ongoing to develop methods that provide the same accuracy for aeromechanics at significantly reduced computational time. Hybrid solvers, applied via bespoke coupling algorithms or within frameworks such as Helios, apply the very costly refined viscous simulations near the body to capture the separation, shear layers and vorticity on the body. The mid- and far-field are resolved using potential-based solvers. These second-generation methodologies correct earlier limitations with boundary conditions, simple geometries (i.e., a single rotor blade), and turbulence modeling. The most successful of these hybrid solvers include free-wake methods [Reference Quon, Smith, Whitehouse and Wachspress69–Reference Moushegian and Smith72] and velocity-vorticity wake formulations [Reference Brown and Line73–Reference Whitehouse and Boschitsch75]. Researchers developing the free-wake coupling report comparable uRANS accuracy (within 2% – 5%) for aerodynamic (and structural dynamics) loads for rotors, with 50% and 90% reduction in computational costs for flight conditions with dynamic stall and attached flows, respectively. The reduction in costs is two-fold with a reduction in mesh requirements and a more rapid convergence of the wake. A recent addition to the Helios suite, the reduced order aerodynamics model (ROAM) is another option to address the high cost associated with high-fidelity CFD simulation, in partular for rotorcraft design. ROAM applies an actuator line model for the rotor(s), while non-lifting bodies (fuselage, nacelles, etc.) are resolved via an immersed boundary method [Reference Wissink, Jude and Jayaraman76]. A wall model, detailed in Jude et al., [Reference Jude, Hosseinverdi, Sitaraman, Peron and Bolsard77] was recently added that improved flowfield predictions and better correlation with high-fidelity results for separated flows.

These are but a few highlights of the advances over the past quarter century in dynamic stall. Overviews of progress in the computation and understanding of dynamic, separated rotorcraft flows can be found in Yeo and Ormiston [Reference Yeo and Ormiston78, Reference Yeo and Ormiston79], van der Wall et al. [Reference van der Wall, Lim, Smith, Jung, Bailly, Baeder and Boyd80], Smith et al. [Reference Smith, Lim, van der Wall, Baeder, Biedron, Boyd, Jayaraman, Jung and Min81], Smith et al. [Reference Smith, Jones, Ayancik, Mulleners and Naughton82, Reference Smith, Gardner, Jain, Peters and Richez83], Gardner et al. [Reference Gardner, Jones, Mulleners, Naughton and Smith65] and Smith [Reference Smith84], in addition to individual publications. By applying the software and hardware advances highlighted in Section 2 and the extant literature, the understanding of and ability to computationally predict a complex dynamic phenomenon, dynamic stall, has been significantly advanced over the past quarter century.

3.2 Aerodynamics of transient environments

Another area of emerging critical importance in this area is vehicle interaction with intense and sometimes large transients that are encountered in airwakes near obstacles, such as in urban environments and ship airwakes illustrated in Figure 20. Flight in these unsteady airwakes has implications for military and civilian use cases from very small uncrewed air vehicles (UAV) to advanced air mobility (AAM) to traditional full-size rotorcraft. Flight through these environments affects flight safety [Reference Goericke, Brenner and Lee85], pilot workload, passenger ride quality and operational life cycle and maintenance of vehicles [Reference Smith, Bayoumi and Corman86]. Operations near heliports/vertiports in an obstacle-rich area can pose significant safety hazards to nearby nonparticipants [Reference Branco, Stephens, White and Watson87].

Figure 20. Computational examples of the variation of the winds in the wakes of obstacles.

Traditionally, engineers have applied indicial theory or considered gusts that remain relatively small compared to the forward flight speed of the air vehicles. However, this is not the case when vehicles operate in the roughness sublayer (RL) of the atmospheric boundary layer (ABL), per Figure 21. In the RL, flight vehicles will rapidly encounter a changing environment composed of shear layers, shed vortices, and recirculation zones that will create equivalent transverse and streamwise gusts [Reference Jones, Cetiner and Smith89].

Figure 21. A representative schematic of the atmospheric boundary layer where AAM and UAV flight paths will be, from Salins et al. [Reference Salins, Crawford and Smith90].

Of similar interest are shipboard rotorcraft flight operations, where the atmospheric boundary layer drives the turbulent environment around the ship so that flight deck operations occur in a highly transient environment (Figure 20(b)). The pilot workload and operational safety can be adversely impacted both in forward flight and hover.

During the past decade, experimental and computational groups have been collaboratively investigating these gusts and environments to gain an understanding of the physics and develop the best computational practices for flight through them [Reference Jones, Cetiner and Smith89]. For example, for transverse gusts interacting with lifting surfaces above gust ratios of 0.5–0.75, there is a distinct change in the physics and, correspondingly, the aerodynamic behaviour compared with lower gust ratios (Figure 22), the response can be accurately predicted with Küssner indicial theory as the interaction is linear. As the gust ratio increases, there is an interaction of the leading and trailing edge vortices, along with significant flow separation, so that the overall aerodynamic interaction is nonlinear and linear indicial theory is no longer accurate.

Figure 22. Wing-gust lift response with increasing gust ratio [Reference Bonnet, Kolpitcke and Smith91]. Experimental data [Reference Corkery, Babinsky and Graham92] for ${\rm{GR}} = 1$ is presented as a dashed line.

Current research efforts are ongoing to explore different kinds of unsteady turbulent interactions, how shapes and extents of transients influence the vehicle behaviour, how to predict and model the transients and how to prevent adverse impacts of flight through these transients. Much of the early research focused on the resolution of the physics for fixed lifting surfaces – that is, the gust interacts with the lifting surface, but the surface does not move in response to the gust interaction. For heavy vehicles with significant inertia these assumptions may be good approximations, but as vehicle weight decreases, these are not valid assumptions. A 10 kg UAV and Blackhawk helicopter will respond quite differently to the same gust ratio. The full gust velocity may be absorbed by the larger Blackhawk vehicle, with little uncommanded vehicle displacement. However, the light UAV will respond to the gust velocity primarily through vehicle displacement, so that the effective gust ratio is reduced [Reference Stutz, Hrynuk and Bohl93, Reference Bonnet94]. Therefore, the inclusion of flight controls and/or rigid body modes need to be considering during computational evaluations of these interactions.

The computational modeling of these interactional flows pose some different issues than those previously discussed for dynamic stall. These interactions are not periodic, but pose a single, unique transient. Thus, while for periodic flows the simulation can include many different cycles of the behaviour to ensure accuracy and study the unsteady details, that is not possible here. A practical approach to these computational assessments is the use of overset grids to ‘fly’ the vehicle or lift component through the unsteady environment. For uRANS, WMLES and LES applications, ensuring that Galilean invariance (GI) of the models is maintained is crucial. To illustrate this, a spanwise station of a lifting surface interacting with a transverse gust at a transitional $Re$ is examined in Figure 23. When the simulation involves a traditional single frame of reference with the aerofoil velocity modeled as an inflow variable as in Figure 23(b), if grid motion (or an overset mesh) does not include GI in the LM transition model, an incorrect flowfield results (Figure 23(a)). The correct flowfield is recovered when the GI terms in the LM model are added (Figure 23(c)). This simple example shows the importance of understanding the numerical implication of moving frames of motion.

Figure 23. $\tilde{R}{{{e}}_{{{\theta t}}}}$ field for stationary and translating NACA0012 aerofoil, from Crawford [Reference Crawford95].

With bodies experiencing these transient conditions, where the unsteadiness is not caused by body motion, maintaining convergence and solution accuracy is dependent on the situation and the solution may become unstable using nomimally best convergence practices from periodic simulations [Reference Crawford95]. For these gusts, for example, the simulation accuracy is dependent on achieving temporal convergence prior to the gust interaction. For changing inflows, such as flight through an urban canyon, the numerical convergence should be monitored during the unsteady simulation, and, if required, numerical options modified to ensure an accurate, stable solution. It was also observed that gust solutions are less accurate when dissipative first-order schemes are applied to improve solution stability [Reference Crawford95]. Traditional measures of convergence in separated flows, typically two orders of magnitude reduction, may not be sufficient in these scenarios. This is an area of ongoing research and development.

The simulation of an unsteady environment through the solution of the Navier-Stokes (N-S) equations is problematic due to the dissipative nature of the computational representation of the N-S equations. If a boundary condition is used to impress an unsteady, changing environment, the velocity and vorticity inputs may fully or significantly disappear when they reach the configuration under analysis. There have been various approaches to mitigate this problem using multiple solvers [Reference Wilbur, Moushegian, Smith and Whitehouse70] and jets to model gusts [Reference Bonnet, Kolpitcke and Smith91], for example. Wales et al. [Reference Wales, Jones and Gaitonde96] introduced a more generic approach to model these unsteady conditions without encountering dissipation. Similar to the additional terms that are added for moving grids, additional source terms that propagate the unsteady conditions can be added to the numerical solver. When properly modified, the solver is able to provide a two-way coupling or interaction where the unsteady environment is felt on the vehicle and the vehicle effects modify the environment. In larger transients and gusts, this nonlinear interaction is important in accurately capturing the physics [Reference Jones, Cetiner and Smith89, Reference Bonnet94, Reference Corkery, Babinsky and Harvey97].