Analysis of the vortex-dominated flow field over a delta wing at transonic speed

Abstract

The present work provides an advancement in the prediction of delta wing flow and an improved understanding of various flow physical phenomena which occur over the wing in transonic conditions. Scale-resolving simulations of the vortex-dominated flow around a sharp leading-edge VFE-2 wing have been performed using the SA-based IDDES model. The complex leading-edge vortex pattern with embedded shocks and subsequent shock-vortex interaction is investigated. A promising accuracy has been achieved using the high-fidelity flow field data provided by the scale-resolving simulation results. Besides the assessment of sensitivity to spatial and temporal resolution, physical aspects are presented, which are not accessible in experimental data in such detail and require scale-resolving simulation approaches. This includes the observation of the vortex system and the shocks in the fully three-dimensional flow field data. Finally, turbulence-related quantities such as eddy viscosity and resolved Reynolds-stresses and their behaviour during the vortex formation and sustaining process are analysed.

Nomenclature

${C_L}$

Lift coefficient

${C_{{m_y}}}$

Pitching moment coefficient

${c_p}$

Pressure coefficient

d

RANS length scale [m]

${d_w}$

Distance to the wall [m]

$\tilde d$

Hybrid length scale [m]

${H_n}$

Normalised helicity

K

Normalised turbulence kinetic energy

L

Characteristic length [m]

$LESI{Q_\nu }$

Index of Resolution Quality

Ma

Mach number

${p_0}$

Total pressure [ $N/{m^2}$ ]

Q

Q-criterion [ $1/{s^2}$ ]

${R_{ij}}$

Normalised specific Reynolds-stress tensor

${R_t}$

Turbulent eddy viscosity over dynamic molecular viscosity

Re

Reynolds number

U

Velocity magnitude [ $m/s$ ]

${u_t}$

In-plane tangential velocity [ $m/s$ ]

$u,v,w$

Velocity components [ $m/s$ ]

${u^{\prime}},{v^{\prime}},{w^{\prime}}$

Velocity fluctuations [ $m/s$ ]

$x,y,z$

Cartesian coordinates [m]

Greek symbol

$\alpha $

Angle-of-attack [deg]

$\Delta $

Subgrid length scale [m]

$\mu $

Dynamic molecular viscosity [ $kg/(ms)$ ]

${\mu _t}$

Turbulent eddy viscosity [ $kg/(ms)$ ]

$\xi ,\eta ,\psi $

Dimensionless Cartesian coordinates

${\omega _x}$

x-vorticity [ $1/s$ ]

$\nabla \rho $

Density gradient [ $kg/{m^4}$ ]

1.0 Introduction

Complex flow fields dominated by vortex systems are generated at extreme flight conditions over low-aspect-ratio wings with highly swept leading edges. Delta wings are used in a variety of aerospace vehicles, such as highly swept wings of a fighter aircraft and for supersonic civil transport. The accurate prediction of the flow physics strongly depends on the treatment of turbulence, which is one major challenge in this case.

A large amount of research has been carried out on aerodynamics of delta wings. A wide variety of conditions can be computed applying classical Reynolds Averaged Navier-Stokes (RANS) models that are very efficient in terms of computational time. However, they are not capable of predicting the flow sufficiently accurate at the flow conditions considered in this work, like transonic speed and high angles of attack. To assess the performance of different RANS turbulence models, vortical flows around a generic triple-delta wing geometry have been investigated at transonic speeds by Werner et al. [1]. At Mach 0.85 large discrepancies connected with vortex-shock interaction have been observed between the results of the various turbulence models. This leads to the conclusion that a common one- or two-equation RANS model, even a complex one, will not provide the accuracy needed in all varieties of highly unsteady vortical flows. Different approaches can be found in literature to overcome the general deficiencies of the eddy-viscosity approaches, such as the Reynolds-stress models and the Quadratic Constitutive Relation [2]. Both approaches try to remedy some of the shortcomings of the linear eddy-viscosity models due to the Boussinesq hypothesis, since the latter assumes isotropy and locality of the Reynolds-stress field. Besides, Moioli et al. [3] aimed to adapt the one-equation SA model specifically to delta wing applications, by using an experimental database and a gradient-based optimisation process. Such development leads to very specialised turbulence models suitable only for specific applications.

On the other hand, it may be possible to resolve turbulence in a more general way like employing Large Eddy Simulation (LES). However, this method is too expensive in terms of computational effort to apply it on a routine basis. At this point, the hybrid RANS/LES approach is considered an alternative to accurately capture the unsteady characteristics of the vortices at reasonable effort by resolving portions of the turbulence spectrum in the main flow detached from the wall, instead of modeling the entire turbulence spectrum.

Hybrid RANS/LES is an emerging field in the analysis of vortical, turbulent structures, and many researchers have been trying to assess the method in order to correctly simulate the aforementioned flows. For example, Peng and Jirásek [4] have performed RANS and hybrid RANS-LES computations taking the flow around the VFE-2 delta wing at a low-speed subsonic Mach number $Ma = 0.14$ and an angle-of-attack (AoA) $\alpha { = 23^ \circ }$ into account. The hybrid RANS-LES computation has reproduced the mean flow more reasonable than the RANS computation, conducted with the Spalart-Allmaras (SA) model and the Explicit Algebraic Reynolds Stress Model (EARSM), in view of the resolved secondary vortex and the predicted surface pressure. Moreover, numerical simulations of the flow at subsonic conditions for the VFE-2 delta wing configuration with rounded leading-edges have been carried out by Cummings and Schütte [5] using RANS and several hybrid turbulence models. The simulated flow field resulting from SA-DDES simulations has shown significant improvements compared to the other hybrid turbulence model simulations, and the results showed promise for gaining an understanding of the flow field.

The present work provides an advancement in the prediction of the vortex-dominated flows for the understanding of the various flow phenomena that occur over the delta-wing particularly at transonic conditions. Simulations of the sharp leading-edge VFE-2 wing [6] have been performed with $M{a_\infty } = 0.8$ , $R{e_\infty } = 2{\rm{e}}6$ and $\alpha { = 20.5^ \circ }$ . Since the flow separation, which forms the initial stage of vortex formation, is fixed by the sharp leading-edge, the main challenge within the simulation is to correctly produce formation and further development of the vortical flow system along the wing surface, which is primarily affected by the treatment of the turbulence in terms of modeling or resolving turbulent eddies. In case of low-aspect-ratio delta wings, the shear layer emanating from the leading edge is rolling up to form a leading-edge vortex, thereby inducing additional velocities on the wing surface. The generated vortex sheet is highly influenced by the pressure gradients in its vicinity, and its separation at the swept leading edge causes a local low-pressure region on the suction side, which contributes to the overall lift [7]. The so-called vortex lift has a limiting AoA at which the vortex breaks down. This consists of an abrupt change in the flow topology where the flow decelerates and diverges. The flow physics over a delta wing further complicates in the presence of single or multiple shock waves. The interaction between shocks and vortices, which can trigger the appearance of the vortex breakdown [8], then becomes a relevant phenomenon to analyse. Both the complex vortex system and the shock-vortex interaction have been investigated in detail.

The vortex-dominated flow has been investigated comparing Improved Delayed Detached Eddy Simulation (IDDES) and unsteady RANS (URANS) results with experimental data [6]. The IDDES method based on the Spalart-Allmaras One-Equation Model with corrections for negative turbulent viscosity (SA-neg) is applied in the scale-resolving computations, whereas the One-Equation SA-negRC (with corrections for Rotation/Curvature as well) is employed to close the RANS equations [9]. The IDDES approach has been selected for the scale-resolving simulations in the present work because on a sharp leading edge delta wing the generation of turbulence usually starts soon after the leading edge separation in the separated shear layer yet in close distance to the wall. The IDDES model essentially switches to URANS mode in the wall layer while running in LES mode in the off-wall region. Mesh refinement in the onset of the turbulent shear layer emanating from the leading edge helps to drive the IDDES model into wall-modeled LES mode and thereby benefits from its capability of controlling the transition between RANS and LES in the region immediately after separation by trying to mitigate the grey-area effects. The sensitivity of the results to temporal and spatial resolution will be addressed and discussed as well. The simulations have been performed employing the DLR TAU-Code [10].

2.0 Test case analysis

A detailed numerical investigation of the VFE-2 delta wing configuration at transonic conditions with $M{a_\infty } = 0.8$ , $R{e_\infty } = 2{\rm{e}}6$ and $\alpha { = 20.5^ \circ }$ has been performed. The selected free-stream conditions are representative for highly agile delta wing aircraft and therefore relevant to aerodynamic design topics such as manoeuverability, stability and control.

2.1 Geometry and mesh

The VFE-2 wing with a leading-edge sweep angle of $\phi { = 65^ \circ }$ features a sharp leading-edge. The experimental data consist of PSP and PIV measurements provided by DLR [6]. The operating conditions used in the CFD analysis are set according to the experimental conditions. Dimensionless Cartesian coordinates are introduced as follows, $\xi = x/L$ , $\eta = y/(x \cdot tan(\phi ))$ and $\psi = z/(x \cdot tan(\phi ))$ , where L is the characteristic length, the chord length in the symmetry plane.

The computational mesh, denoted “extra-fine” in Table 1, is employed to investigate the VFE-2 geometry. It is shown in Fig. 1. The outer boundary is formed by a spherical farfield boundary located 50 L from the wing. The unstructured mesh consists of about 30 million cells and features up to 30 prism layers along the walls with the wall-normal growth factor equal to 1.1 and the first cell thickness such that ${y^ + } \lt 1$ . The mesh is symmetric to the plane $y = 0$ . The cells sizes vary within the computational domain and the finest cells, whose size is around 0.002 L, are located within the vortex region as well as close to the leading-edge to resolve the shear layer onset. In order to capture the development of turbulent scales, the mesh refinement follows the vortex region over the wing. The number of grid points inside the vortex diameter at chordwise location $\xi = 0.4$ provides a justification that the grid resolution can be assumed to be adequate for the given flow. As suggested by Landa et al. [11], the vortex diameter ${d_\omega }$ has been computed from the vorticity distribution ${\omega _x}$ and has been found to be at the order of ${d_\omega } \approx 0.3L$ . ${N_\omega } \approx 150$ is therefore the number of grid cells inside the vortex diameter ${d_\omega }$ . Although the cell size slightly increases along the wing, the ratio of the vortex diameter to the cell size rises due to the vortex expansion.

Table 1. Details of the meshes

Figure 1. VFE-2 geometry and mesh.

2.2 Mesh convergence study – URANS and IDDES runs

A grid-resolution study has been performed in order to analyse the grid effects and keep them as low as possible. Table 1 summarises the main mesh characteristics. The finest cell size has been halved step by step. The URANS simulations have been performed with all the four meshes, whereas the IDDES runs have been performed employing only the two finest meshes to demonstrate the suitability of the grid resolution.

Figure 2 shows the relative deviation of lift and pitching moment coefficients with respect to the finest mesh for which results will be presented. Taking both approaches (IDDES and URANS) into account, the plot of the relative deviation demonstrates the monotonous and effective reduction of the grid effects. In particular, the mesh convergence is clearly visible comparing the RANS results. Besides, the “fine-I” column in Fig. 2 shows that there is no relevant difference in the prediction of the aerodynamic coefficients between the IDDES results achieved with the “fine” and “extra-fine” mesh, confirming mesh convergence. Further considerations on the grid-sensitivity will be made in Section 4.1 with the analysis of mean surface pressure coefficient.

Figure 2. Relative deviation bar plot of the aerodynamic coefficients. “I” denotes the comparison between the IDDES results on the “fine” and “extra-fine” mesh [12].

3.0 Modeling and solution

The DLR TAU-Code has been used to perform the CFD simulations. A brief code introduction with an algorithmic overview, which shortly describes the code functionality, has been provided by Gerhold [13]. In the present setup the flow solver TAU solves the compressible, three-dimensional, time-accurate Reynolds-Averaged or filtered Navier-Stokes equations using a finite volume formulation. Governing equations are presented within the work of Langer et al. [10] and will not fully be shown here for sake of conciseness.

3.1 Model equations and approaches

The SA-neg model is based on the following single equation by Spalart and Allmaras [14] for the eddy viscosity

(1) \begin{align}\frac{\partial }{{\partial t}}(\rho \tilde \nu ) + \vec u \cdot \vec \nabla (\rho \tilde \nu ) = \vec \nabla \cdot \left( {\frac{{\mu + \rho \tilde \nu }}{\sigma }\vec \nabla \tilde \nu } \right) + \rho \frac{{{c_{b2}}}}{\sigma }{(\vec \nabla \tilde \nu )^2} + {c_{b1}}\rho \tilde S\tilde \nu - {c_{w1}}{f_w}\rho {\left( {\frac{{\tilde \nu }}{d}} \right)^2},\end{align}

where the RANS length scale $d = {L_{RANS}}$ in the destruction term is the distance to the nearest wall [15]. In the SA-neg version [15], turbulent eddy viscosity is set to zero in case it becomes negative in Equation (1). This version is used to improve stability and robustness without changing the (converged) results of the SA model.

The SA turbulence model often provides excessive eddy viscosity production in the vortex core with implications on the unburst vortex size, type and velocities. Shur et al. [16] proposed a streamline curvature correction (SA-RC), applied in the current work within the URANS computations, which alters the source term with a rotation function. In the RANS context a comparison with and without RC has been performed resulting in the selection of the SA-negRC model for the presented data as it produces superior results.

In the IDDES model, ${L_{RANS}}$ is replaced with $\tilde d = {L_{IDDES}}$ , which is defined as

(2) \begin{align}{L_{IDDES}} = {\tilde f_d}(1 + {f_e}){L_{RANS}} + (1 - {\tilde f_d}){L_{LES}}\end{align}

where ${L_{RANS}}$ and ${L_{LES}}$ are the RANS ( ${L_{RANS}} = {d_w}$ for the SA model) and LES ( ${L_{LES}} = {C_{DES}}\Psi \Delta $ with the function $\Psi $ defined in literature [9]) length scale respectively. In order to define a proper IDDES blending function ${\tilde f_d}$ the authors of IDDES propose a set of semi-empirical functions designed to provide for both a correct performance of the WMLES and DDES branch [9]. The objective is to prevent an excessive reduction of the RANS Reynolds-stresses as could be caused by the interaction of RANS and LES regions in the vicinity of their interface [9]. Figure 3 shows the instantaneous hybrid length scale over RANS length scale ( $\tilde d/d$ ). It illustrates where the IDDES approach switches from RANS to LES. The thin regions close to the wall are fully modeled by the RANS mode and a ratio around unity is expected. Thereby, the shear layer transition takes place in the LES region. Since the vortex core region is covered by the LES mode, the RC correction has not been applied within the hybrid RANS/LES model. Spatial and temporal resolution in this zone have been investigated and will be discussed in the following.

Figure 3. Instantaneous hybrid length scale over RANS length scale $\tilde d/d$ and instantaneous LES Index of Resolution Quality $LESI{Q_\nu }$ at chordwise locations $\xi = 0.2,\ 0.6,\ 0.95$ .

To demonstrate that the current mesh resolution allows to resolve a major part of the turbulence spectrum, the LES Index of Resolution Quality has been included in Fig. 3. The $LESI{Q_\nu }$ relates turbulent viscosity to laminar viscosity using the formulation proposed by Celik et al. [17]. The $LESI{Q_\nu }$ is a dimensionless number between zero and one. It is calibrated in such a way that the index behaves similar to the ratio of resolved to total turbulent kinetic energy. An index of quality greater than 0.8 is considered a good LES, 0.95 and higher is considered DNS [18]. Overall, the plots indicate a good spatial resolution of the vortex region. Slight weaknesses appear in the downstream parts of the shear layer which, however, are too far downstream to have a significant impact on the vortex core formation and also on the shock-induced phenomenon.

3.2 Numerical approach

An implicit dual-time stepping approach, employing a Backward-Euler/LUSGS implicit smoother, has been used in the unsteady simulations. To ensure convergence of the inner iterations in the IDDES runs, Cauchy convergence criteria of several quantities, namely volume averaged turbulence kinetic energy, maximum Mach number and certain aerodynamic coefficients, with tolerance values of $1{\rm{e}} - 5$ have been applied. The matrix dissipation model has been selected while performing the computation of the fluxes with a central scheme. However, in hybrid RANS/LES the artificial dissipation has been reduced to prevent excessive damping of the resolved turbulent structures and a (hybrid) low-dissipation low-dispersion discretisation scheme (LD2) has been used, as suggested by Probst et al. [19].

Figure 4 shows the temporal evolution of the lift coefficient ${C_L}$ and provides the duration of simulated time for the different cases. Time has been characterised by the Convective Time Unit (CTU) which is computed as follows

(3) \begin{align}T = L/{U_\infty } = 1\quad CTU,\end{align}

where L is the characteristic length and ${U_\infty }$ the free stream velocity.

Figure 4. Time history of the lift coefficient, instantaneous (on the left) and running time averaged (on the right).

For URANS, the time step size equals $100\mu s \approx 2.64{\rm{e}} - 2\,CTU$ . Ten CTU have been computed before starting the time-averaging in order to overcome the initial transient and five flow-through times have been taken into account in the computation of the mean flow field values.

To fully resolve the convective transport and consequently to capture the flow characteristics accurately, the maximum allowed time step size for IDDES runs has been computed. On the “extra-fine” mesh, the time step size $\Delta t = 1\mu s = 2.64{\rm{e}} - 4\,CTU$ adequately resolves the time scales of the energy containing eddies in the region of interest and keeps the convective Courant-Friedrichs-Lewy (CFL) number lower than unity [20], as shown in Fig. 5. To accelerate the simulation and save computational time, a simulation has been performed increasing the time step to $\Delta t = 10\mu s = 2.64{\rm{e}} - 3\,CTU$ resulting in a convective CFL number approximately one order higher in magnitude. The IDDES run using the “extra-fine” mesh has been initialised from scratch with a steady run and afterwards five CTU have been performed with the time step size equal to $\Delta t = 10\mu s$ . Ten overflows have been considered to compute the mean values of the flow properties. Besides, after five CTU performed with the time step size equal to $\Delta t = 10\mu s$ , five CTU have been computed with $\Delta t = 1\mu s$ before starting the time-averaging, to overcome the initial transient. Finally, in this case ten overflows have been considered as well to compute the mean flow field values. The effect of time step size can be observed comparing the black and the green lines in Fig. 4. The lift coefficient remains almost constant in time with the coarser temporal resolution whereas it starts to drop suddenly when reducing the time step. The increased temporal accuracy provides an overall improvement of the achieved numerical results.

Figure 5. Mean convective CFL number at chordwise locations $\xi = 0.2$ , comparison between the IDDES results with the time step size $\Delta t = 10\mu s$ and $\Delta t = 1\mu s$ [12].

Figure 4 shows the evolution of the IDDES run with the “fine” mesh as well. These simulations have been initialised with URANS results and no difference has been found in the duration of the initial transient. It is also noteworthy that the time step sizes have been doubled (see the legend of the Fig. 4) as the size of the finest cell has also been doubled compared to the “extra-fine” reference mesh, as summarised in Table 1.

4.0 Results

The flow physics is described and illustrated structured in the following sub-topics. The analysis of mean flow features is presented in Section 4.1 providing a comparison of several aspects, the instantaneous flow features are discussed in Section 4.2. In Section 4.1.1 an emphasis is given to the validation of the numerical results by comparison with experimental data regarding the mean surface pressure. Besides, the location and intensity of the shock waves are discussed. In Section 4.1.2 the vortex flow pattern and the vortex-shock interaction phenomenon are then investigated with the help of the scale-resolving results. Finally, the turbulence-related quantities, such as instantaneous eddy viscosity, resolved turbulence kinetic energy and components of the Reynolds-stress tensor are presented in Section 4.3.

4.1 Mean flow features

Over the delta wing, the flow separates at the leading-edge and subsequently rolls up to form a stable, separation-induced primary vortex. The flow reattaches on the surface as the primary vortex is rolling up and the spanwise flow underneath subsequently separates a second time to form a counter-rotating secondary vortex outboard of the primary one. The direct effect of the vortices is an additional suction footprint and, consequently, increased lift which eventually results in a non-linear dependence of the lifting force on the angle-of-attack. The understanding and prediction of vortex and shock waves, their generation and evolution, are of essential importance. The interaction between vortices and shock waves is a crucial feature of the flow physics at transonic conditions and therefore is assessed in detail.

Figure 6. Mean surface pressure coefficient $\overline {{c_p}} $ , comparison between experimental and numerical data. The black contour lines indicate the sonic pressure coefficient $\overline {c_p^*} = - 0.43$ [12].

4.1.1 Result validation and shock-wave locations

The suction footprint on the wing surface is mainly caused by high tangential velocity around the inner vortex core. Figure 6 shows the mean surface pressure coefficient. To investigate the prediction of $\overline {{c_p}} $ , several slice planes have been extracted, as indicated in Fig. 6(a). The $\overline {{c_p}} $ distribution along the spanwise direction at chordwise positions $\xi = 0.2,\ 0.4,\ 0.6,\ 0.8,\ 0.95$ is plotted in Fig. 7, comparing experimental data and simulation results (URANS, IDDES with $\Delta t = 10\mu s$ and $\Delta t = 1\mu s$ using the “extra-fine” mesh and IDDES with $\Delta t = 2\mu s$ using the “fine” mesh). The comparison is presented to demonstrate that the best results have been achieved using the IDDES approach with extra-fine (reference) mesh and $\Delta t = 1\mu s$ .

Figure 7. Mean surface $\overline {{c_p}} $ distribution at chordwise locations $\xi = 0.2,\ 0.4,\ 0.6,\ 0.8,\ 0.95$ and at symmetry plane $\eta = 0$ (lower right) comparison between experimental data, URANS and IDDES results [12].

The plots highlight the sensitivity of the results to the temporal and spatial resolution. The $\overline {{c_p}} $ distribution shows that some differences appear between the IDDES results achieved with the two finest meshes (i.e. “fine” and “extra-fine”), even though Fig. 4 had shown almost identical lift coefficients. These discrepancies indicate that additional criteria need to be considered for the mesh convergence since the integral aerodynamic coefficients do not seem to be sufficient. However, the uncertainties of the IDDES results due to mesh resolution are considered significantly reduced and IDDES on the “extra-fine” mesh hence is taken as reference. Figure 7 shows a satisfactory and suitable agreement between the IDDES results ( $\Delta t = 1\mu s$ ) and experimental data. The IDDES results ( $\Delta t = 1\mu s$ ) provide the closest match with the experiment of all numerical results, especially at section $\xi = 0.4$ . Hybrid RANS/LES considerably improves the predicted suction footprint by the vortices, although the experimental results are only roughly reproduced in the front part of the wing. The fine resolution is particularly beneficial in the apex region to reduce the grey area, which will be further discussed in Section 4.3.

The secondary and tertiary vortices are indicated with Roman numbers (II and III) at the smoother peaks of ${c_p}$ at $\eta \approx 0.7$ and 0.85 from chordwise location $\xi = 0.2$ . The tertiary vortex then disappears in the IDDES results with $\Delta t = 1\mu s$ at chordwise location $\xi = 0.4$ and matches the experimental data perfectly, whereas it is still erroneously present in the URANS results. The experimental data indicates the overall absence of a tertiary vortex but, taking the data resolution into account, this phenomenon also cannot be fully excluded. Besides, the secondary vortex peak is weaker at $\xi = 0.8$ in the experimental data and it is not present any more at $\xi = 0.95$ . This means that it breaks down between $\xi = 0.8$ and $\xi = 0.95$ . Figure 7 indicates that in the IDDES with $\Delta t = 1\mu s$ it bursts further upstream ( $\xi \lt 0.8$ ), as it will be also discussed in Section 4.2.

Moreover, it is very interesting to note that increasing the time step size does not considerably affect the surface $\overline {{c_p}} $ caused by the primary vortex. On the other hand, the higher time step $\Delta t = 10\mu s$ leads to a wrong prediction of the secondary vortex due to large unsteady fluctuations.

As it can be seen in Fig. 7 at $\xi = 0.95$ , the trend of the $\overline {{c_p}} $ curve has been captured well by the numerical results except for the middle part (where the primary vortex core is located). The experiment shows a decay in this region whereas the simulation still predicts a strong drop of $\overline {{c_p}} $ caused by the vortex. The vortex in the experimental data is weaker which could indicate the tendency to break down near the trailing edge. Besides, as it will be discussed in Section 4.1.2, close to the trailing edge, the shear layer is not rolling up to form a stable leading-edge vortex over the wing and the suction footprint behind the transported shear layer leading-edge separation abruptly drops. The primary vortex is no longer fed by the generated turbulence and the coherent vortex becomes more vulnerable. A possible explanation of this phenomenon will be also given in Section 4.2 by means of unsteady flow characteristics.

The investigation of the supersonic area over the wing and consequently the presence of the shock waves can be observed in Fig. 6 where the sonic pressure coefficient $\overline {c_p^*} = - 0.43$ is indicated by the black contour lines. The IDDES approach replicates the experimental results best in terms of predicted supersonic area. In front of the secondary vortex a cross-flow shock wave appears, where Fig. 6 shows that the outboard directed flow underneath the primary vortex close to the surface is supersonic. The so-called “lambda” shock is also illustrated in Fig. 9. Besides, a shock wave typical for transonic free-stream conditions is captured close to the wing apex, where critical flow conditions are reached. Finally, the contour lines of the sonic pressure coefficient in the centre of the delta wing give the indication of shock waves in proximity of the sting fairing: a first one located at about $\xi = 0.55$ , a second one at about $\xi = 0.75$ and a third one closer to the trailing edge downstream of $\xi = 0.9$ . The same shock waves are shown and enumerated in Fig. 11(a). Besides, Fig. 7 shows the mean surface $\overline {{c_p}} $ distribution at symmetry plane $\eta = 0$ to quantify the intensity of the aforementioned shock waves. Although the simulation results predict a higher pressure peak close to the wing apex ( $\xi = 0$ ) and to the sting tip ( $\xi \approx 0.62$ ), the pressure trend and the position of the shock waves above the wing are captured well by the numerical data.

4.1.2 Vortex pattern and vortex-shock interaction

Considering the simulation results valid based on the discussion above, the detailed analysis of the IDDES results with $\Delta t = 1\mu s$ can be performed. The vortex pattern itself and how it is modified by the shock-vortex interaction needs to be understood in detail. For this purpose, Fig. 8 shows the field of mean x-direction vorticity $\overline {{\omega _x}} $ at chordwise locations $\xi = 0.5,\ 0.6,\ 0.7,\ 0.75$ from both experiment and IDDES. PIV data [6] have been collected only at these locations illustrated even in Fig. 6(b). Figure 9 then shows the mean density gradient magnitude $||\overline {\nabla \rho } ||$ , the normalised mean x-direction velocity $\overline u /{U_\infty }$ and the normalised mean in-plane tangential velocity $\overline {{u_t}} /{U_\infty }$ distribution, where ${u_t} = \sqrt {{v^2} + {w^2}} $ . Further, Fig. 10 shows the primary vortex core line extracted considering the local maximum mean x-velocity.

Figure 8. Mean x-direction vorticity $\overline {{\omega _x}} [1/{s^2}]$ distribution at chordwise locations $\xi = 0.5,\ 0.6,$ $0.7,\ 0.75$ comparison between experimental data (top) and IDDES results (bottom). The black contour lines in the experimental data are related to the divergence of the in-plane velocity vectors, $\nabla \cdot \vec u[1/{s^2}]$ [12].

Figure 9. Mean density gradient magnitude $||\nabla \rho ||$ , mean normalised x-direction velocity $\overline u /{U_\infty }$ and mean normalised in-plane tangential velocity $\overline {{u_t}} /{U_\infty }$ distribution, IDDES results at chordwise locations $\xi = 0.5,\ 0.6,\ 0.7,\ 0.8,\ 0.9$ [12].

Figure 10. Primary vortex core line extracted from the IDDES ( $\Delta t = 1\mu s$ ) results.

Figure 11. IDDES ( $\Delta t = 1\mu s$ ) results for the VFE-2 geometry with $M{a_\infty } = 0.8$ , $R{e_\infty } = 2{\rm{e}}6$ , $\alpha { = 20.5^ \circ }$ .

Since a vortical structure comes from the concentration of the low energy parts of the flow, initially contained in the boundary layer, the vortex can also be seen as a region of rotational flow whose stagnation conditions, especially the pressure, are inferior to those of the outer flow field. Thus, the vortex is like an energy “hole” or entropy “hill” which will be more sensible to disturbing entities, such as shock waves [21]. The unburst vortex is characterised by a coherent structure with high values of axial and rotational velocity. As it can be seen in Fig. 8, the IDDES results capture the primary and the secondary vortex in a very accurate matter. Between $\xi = 0.6$ and $\xi = 0.7$ , the x-direction vorticity starts to decrease, as the vortex does not expand uniformly in the streamwise increase of the $\eta $ -coordinate. As the self-induction breakdown theory explains [22], the axial vorticity of the vortex induces an azimuthal velocity, which in turn tilts the vorticity vector towards the azimuthal direction. Due to the gradient of azimuthal vorticity caused by the increased circulation, the leading-edge vortex expands radially. This is followed by a change of sign of the azimuthal vorticity caused by the region downstream of the vorticity gradient which rotates slower than the upstream region. These actions proceed together up to the turning point where the vortex filaments turn inward onto themselves causing a change of sign of the axial vorticity. This phenomenon is visible in a slice plane close to the trailing edge in Fig. 11(b), where the instantaneous x-vorticity is shown. It means that, after the interaction with the shock upstream of the sting fairing ( $\xi = 0.62$ ), the vortex appears more vulnerable and prone to breakdown. The detached shock caused by the sting tip interacts with the vortex core and causes a decrease of the x-velocity (see Fig. 10(b) at chordwise $\xi = 0.62$ ) with a subsequent reduction of the suction footprint also visible in Fig. 7 at $\xi = 0.8$ . The same behaviour of ${\overline u _{max}}$ can be seen in proximity of the other shock-vortex interactions (at about $\xi = 0.2$ and $\xi = 0.35$ ) as already discussed in Section 4.1.1. The vortex core loses velocity and kinetic energy across the shock. The width of the shock wave is of the same order of the magnitude as the local mean free path of fluid particles. Hence, the shock wave can be regarded as a sharp discontinuity in common aerodynamic flow and, across the shock, the Mach number and the normal velocity suddenly drop. A more gradual drop can be seen in the figures because the mean variables are plotted. This means that the shocks location fluctuates slightly on the wing over time. Besides, the in-plane (tangential) velocity ideally is parallel to the shock surface, which means that it should not be altered significantly by the shock wave, as Fig. 9 demonstrates.

At the present flow conditions the shock wave is not strong enough to abruptly trigger the vortex breakdown. The intensity of the shock close to the sting fairing onset increases as the angle-of-attack rises since it is directly depended on the incoming Mach number. At transonic conditions, the increase of the angle-of-attack may thus induce the vortex breakdown. If, like in this case, the shock wave is not strong enough to induce a vortex breakdown after the interaction with the shock, the leading-edge vortex tends to recover its stability and overcomes the perturbation. Since the main vortex is fed continuously along the wing by the shear layer emanating from the leading edge (illustrated in Fig. 11(b)), the x-velocity of the vortex consequently starts to increase again in the vortex core, as it can be seen in Figs 9 and 10(b).

In Fig. 8 also regions of negative divergence of the in-plane velocity vectors ( $\nabla \cdot \vec u$ ) are marked, which gives an indication of the location and shape of a cross-flow shock wave in the experimental results. These shocks are well captured by the IDDES results and the distribution of $\overline {{u_t}} $ is altered, as it can be seen from the iso-contour lines of $\overline {{u_t}} $ in Fig. 9. The magnitude of the mean density gradient $||\nabla \rho ||$ in Fig. 9 shows this phenomenon as well. A complex shock system is located beneath the leading-edge vortex, as it has been introduced in Section 4.1.1. The “lambda” shock, predicted from the experimental data as well, is mainly caused by the acceleration of the cross-flow in this region and its presence deeply alters the local flow topology. Driven by the interaction of this shock with the boundary layer on the upper side of the wing, the flow separates and feeds the secondary vortex [23]. Moreover, the leading-edge vortex takes a kidney shape for $Ma = 0.8$ and cross-flow shocks appear around the leading-edge vortex and on the top of the shear layer. They are assumed to be caused by the curvature of the flow trajectory leading to an acceleration up to thermophysical conditions which cannot be sustained [23]. These shock waves affect the behaviour of the velocity curve with reversed flow upstream of the sting fairing shock shown in Fig. 10.

Figure 9 shows further phenomena which occur in the rear part of the wing. As mentioned in Section 4.1.1, the secondary vortex breaks down within the range $0.7 \lt \xi \lt 0.8$ , and is not visible at $\xi \gt 0.8$ anymore. Besides, here the primary vortex is not further fed by the shear layer and becomes more vulnerable. The streamwise boundary layer separation highlighted at $\xi \gt 0.8$ might be considered as starting point for the vortex breakdown. Indeed, the breakdown of the secondary vortex generates a chaotic motion of the turbulent shear-layer and prevents rolling up to sustain the primary vortex, as it can be seen at $\xi \gt 0.9$ .

4.2 Instantaneous flow features

Figure 11(a) shows an illustration of the vortices by plotting the instantaneous Q-criterion iso-surface coloured by the normalised helicity ${H_n}$ . The sign of helicity (positive in blue and negative in red) indicates the sense of rotation to separate primary from secondary and tertiary vortices [24]. The IDDES results in Fig. 11(a) allow for a qualitative assessment of the turbulence resolution in the LES regions. Turbulent fluctuations are clearly visible in the vortices. The level of resolution seems to be adequate to thoroughly investigate the flow physics, a large spectrum of turbulent structures has been resolved by the grid.

Moreover, as already introduced, at the transonic speed of $Ma = 0.8$ the flow is much more complex compared to subsonic conditions, because the flow above the delta wing reaches supersonic speeds and shock waves appear. Figure 11(a) with the iso-surface of the density gradient in x-direction $\nabla {\rho _x}$ plot confirms that three main shocks are present over the delta wing in proximity of the sting fairing. Other shock waves are also present close to the wing apex, as already pointed out.

The formation onset of the primary vortex caused by the separated shear-layer is visualised in Fig. 11(b) where the instantaneous total pressure and the primary vortex stream-traces are shown. The total pressure is shown since it represents an energetic quantity and gives an indication of the vortex strength and position. Streamlines of different colours are used to visualise the individual effects. The primary vortex formation starts immediately downstream of the wing apex and the high velocity core of the primary vortex is exclusively built from flow coming off the shear-layer that separates at the wing apex as indicated by the black streamlines. The primary vortex then grows in diameter because it is fed by the shear-layer all along the wing, as it can be seen from the green, brown and purple stream traces in Fig. 11(b). This flow reinforces the main core by rotating around it and feeding it with kinetic energy. In this way the vortex is sustained, remains coherent and the axial velocity increases.

Figure 11(b) also shows the instantaneous x-vorticity, which distinguishes between the primary and the secondary vortex by the rotational direction and the red stream traces of the secondary vortex (in the right half). This plot can be used to observe formation and breaking of the secondary vortex. The secondary vortex formation is not directly caused by the shear-layer separation. The stream traces of the particles that form the secondary vortex do not come from the leading edge in proximity. The secondary vortex is more a secondary effect caused by the primary vortex. Both are closely connected and the secondary vortex only appears well structured if the primary vortex is strong enough. Then, as it can be seen in Fig. 11(b), the secondary vortex breaks down in the second half of the wing downstream of the sting tip and this phenomenon makes it impossible for the shear layer to roll up fully feeding into the primary vortex.

These considerations lead to the hypothesis explaining the process of primary vortex breakdown to be seen at higher angles of attack. Over the wing, the flow separates at the leading edge and subsequently rolls up forming a stable, separation-induced primary vortex. Typically, the flow reattaches when the primary vortex reaches the wing surface. Under certain conditions, the spanwise flow below the primary vortex subsequently separates a second time to form a counter-rotating secondary vortex outboard of the primary one. The secondary vortex is well structured and clearly visible only at mid-wing, though its suction footprint starts about at the same location as the primary vortex (see Fig. 6(b)). Feeding the shear layer into the primary vortex is supported by the presence of the secondary vortex. When the bow shock generated by the sting tip interacts with the primary vortex, its flow conditions change, and the secondary vortex becomes not sustainable anymore. At this point, the boundary layer separates in streamwise direction in proximity of the lambda shock forming a recirculation zone. The fluid previously forming the secondary vortex then breaks into turbulent motion of smaller scales. The shear layer thereby is prevented from rolling up and feeds the incoherent turbulent flow instead of the primary vortex. Consequently, the primary vortex, which loses its source of kinetic energy, becomes more vulnerable. This process can be regarded as the onset of primary vortex breakdown, which only is fully established at slightly higher angles of attack as known from the experimental data.

4.3 Turbulence-related variables

Figure 12 shows the contours of ${R_t} = {\mu _t}/\mu $ , where $\mu $ is the molecular dynamic viscosity. The modeled turbulent eddy viscosity, ${\mu _t}$ , is compared for the chordwise locations $\xi = 0.2,\ 0.4,\ 0.5,\ 0.6,\ 0.7,\ 0.8,$ and $\quad 0.95$ between URANS and IDDES. Overall, the URANS approach produces very high levels of turbulent eddy viscosity in the vortex core. The region with large ${\mu _t}$ values in RANS corresponds to vortex motions with high production of turbulence energy due to strong flow rotation and deformation [4]. As it can be seen in Fig. 12, the RC correction avoids the excessive eddy viscosity production in the front part of the wing. Unfortunately, it is not sufficient in the rear part of the vortex region, where the turbulent viscosity production is by far too high. Regarding the hybrid RANS-LES, it is worth to recall Fig. 3 that shows the instantaneous hybrid length scale over RANS length scale ( $\tilde d/d$ ). The IDDES approach employs the SGS eddy viscosity in the off-wall region and the IDDES ${R_t}$ -contours show that the SGS eddy viscosity is much smaller than its RANS counterpart. The relatively low level of modelled ${\mu _t}$ in LES mode is associated with the local fine grid resolution required for LES. On the other hand, the region with higher ${\mu _t}$ in LES mode indicates strong local flow rotation/deformation and/or coarse grid resolution, usually inducing intensive energy dissipation of resolved large-scale turbulence [4]. A slight grid refinement could be indicated in the rear part of the wing close to the leading edge where the shear layer separates and rolls up.

Figure 12. Instantaneous turbulent eddy viscosity, represented by ${R_t}$ , comparison between URANS and IDDES results.

Figure 13. Resolved turbulence-related quantities along the wing from the IDDES results (note the different scale in the legend boxes).

Figure 13(a) shows the resolved turbulence kinetic energy normalised by the squared free-stream velocity ${U_\infty }$ . As it can be seen from the first slice plane and assuming the grey-area issue, K is probably under-predicted in the initial development of the vortical flow. The formation of resolved turbulence is slightly delayed. This phenomenon is the so-called grey-area problem rooted in the IDDES modeling and it affects the downstream development of the turbulent process including the vortex breakdown onset position. The grey-area problem happens because the formation of the vortex initially takes place in near-wall layers that are modeled by the RANS mode. The transition of the shear layer between the RANS and the LES mode is then supposed to be the main reason for the slight discrepancies highlighted in the $\overline {{c_p}} $ results in the front part of the wing and it generates a rather stiff resolved vortex motion leading to a delayed vortex breakdown. In the rear part of the wing, the resolved turbulence kinetic energy shows the separated shear layer and turbulent downstream transport close to the leading edge, whereas a cross-flow shock alters the K distribution below the primary vortex.

Figure 13 also shows the components of the normalised specific Reynolds-stress tensor ${R_{ij}} = \overline {u_i^{\prime}u_j^{\prime}} /U_\infty ^2$ along the wing. ${R_{ij}}$ represents the intensity of the turbulent fluctuations along the three directions and their covariance. The components of the Reynolds-stress tensor have been normalised by the free-stream velocity and the local mean density to obtain the same order of magnitude of K and to show the turbulent fluctuations without considering the contribution to the energy transportation, respectively. Besides, knowing the high values of the axial velocity in the vortex core shown in Section 4.1.2, it can be affirmed that the density considerably drops in the vortex region. For this reason, the contribution of the turbulence kinetic energy to the total kinetic energy becomes negligible in the vortex core. It acts mainly around the primary vortex core and within the shear ayer where the density value is at least one order higher in magnitude. ${R_{11}}$ and ${R_{22}}$ are illustrated in Fig. 13(b). ${R_{11}}$ shows the turbulent behaviour of the transported turbulent shear layer. The turbulent motion becomes more intense once the secondary vortex breaks and the turbulence kinetic energy is then transported downstream without feeding the primary vortex, as explained in Section 4.2. ${R_{22}}$ indicates that the fluctuations are generated from the leading edge. It is also the main origin of high turbulence kinetic energy in the vortex core. Figure 13(c) shows the covariance ${R_{12}}$ and the normal component ${R_{33}}$ . ${R_{12}}$ can be mainly used to visualise the core of the primary and secondary vortices, where negative values are located (positive on the other wing side), and to identify the boundary of the primary vortex, where the highest values are present. The z-direction of these turbulence fluctuations is visible from ${R_{33}}$ . The elements transported downstream over the wing do not roll up to form the primary vortex but tend to grow moving away from the wing and influencing the primary vortex. Finally, the covariances ${R_{13}}$ and ${R_{23}}$ are shown in Fig. 13(d). The same phenomenon identified with ${R_{12}}$ can be also discussed with the ${R_{13}}$ component. ${R_{13}}$ further indicates the location of the shear layer as well as its thickness and the coherence of the vortex. In the rear part it becomes less coherent and tends to breakdown. Finally, as it can be seen from the legend scale, ${R_{23}}$ is the strongest covariance component and it mainly appears in the shear layer where the complex process of separation and rolling up appears. It is also negative in the core centre (positive on the other wing side), but it is not as localised and compact as ${R_{12}}$ .

5.0 Summary and conclusions

The scale-resolving simulation has provided detailed insight to the transonic flow field around the delta wing. The flow physics has been described, explained and illustrated in detail divided into the analysis of the mean and instantaneous flow features.

Provided adequate cell and time step sizes, the simulation reveals important flow features and represents a valuable method to improve the understanding of the physics of delta wings at transonic conditions. Therefore, the sensitivity of the results to temporal and spatial resolution has been addressed and discussed to ensure that the presented data yields a high prediction accuracy. Numerical results have further been validated by comparing them with the experimental data regarding the mean surface pressure coefficient. Besides, the behaviour of shock waves has been investigated. The surface pressure and the position of the shock waves have been captured well by the numerical data.

The vortex flow pattern and vortex-shock interaction have then been assessed in detail. The scale-resolving method is capable of producing the interaction between the leading-edge vortex and the shock waves. The magnitude of the mean density gradient has been shown to provide insight to the vortex and shock structures and to analyse their characteristics.

The formation of the primary vortex caused by the separated shear layer emanating from the leading edge has been visualised by streamlines. Formation and breakdown of the secondary vortex have been examined as well, providing a hypothesis to explain the process and mechanism for vortex breakdown.

The turbulence-related variables, such as instantaneous eddy viscosity, resolved turbulence kinetic energy and the components of the Reynolds-stress tensor have been discussed. It has been shown that the initial development of resolved turbulence is slightly delayed. Resolved K is predicted insufficiently in the very initial stage of the vortical flow but thereby affects the downstream development of the turbulent process, including the vortex breakdown onset position, which has been predicted slightly too far downstream on the wing by the IDDES results. Being situated within the well-known range of grey area issues, this observation provides potential to introduce slight adaptations into the turbulence model within future work. The components of the Reynolds-stress tensor have been shown to discuss the turbulent behaviour of the flow and to visualise several phenomena, such as the thickness and location of the shear layer and shape and size of the vortex.

IDDES can be then considered as a promising approach to simulate the flow around a delta wing even at transonic conditions and this approach might serve as a reference model for such test cases, even at more challenging conditions and configurations. In continuation of this analysis, further steps will focus on analysis of the vortex structure and the unsteady nature of the flow.