Investigation of transonic aerodynamics on a triple-delta wing in side slip conditions

The vortex-dominated flow around the triple-delta wing ADS-NA2-W1 aircraft is investigated in order to achieve a better understanding of the flow physics phenomena that occur over the aircraft particularly at the transonic speed condition. Both URANS and scale-resolving DDES have been employed in order to explore the range of suitability of current CFD methods. The Spalart–Allmaras One-Equation Model with corrections for negative turbulent viscosity and Rotation/Curvature (SA-negRC) is employed to close the RANS equations, whereas the SAneg-based DDES model is applied in the scale-resolving computations. The DLR TAU-Code is used to perform the numerical simulations. The deficiencies of the URANS results are illustrated and promising improvements are reached employing the SAneg-DDES numerical method. The hybrid method results show great advancement in the prediction of the multiple-delta wing flow by revealing physical aspects which have not been seen from URANS with sufficient accuracy like vortex–vortex interaction and shock-vortex interaction. These phenomena furthermore explain in a clear way the improved prediction of the surface pressure coefficient over the aircraft and consequently of the aerodynamic force and moment coefficients.


Introduction
Combat aircraft configurations typically feature low aspect ratio wings with highly swept leading edges in order to provide the required agility. At extreme flight conditions, complex flow fields dominated by vortex systems, which are challenging for numerical flow simulations, are generated. The key challenge for producing the flow correctly in the numerical simulation is the treatment of turbulence.
The investigation of leading-edge vortices of swept wings with low aspect ratio has been subject to several research projects within the past decades [1]. Also, unsteady phenomena like the vortex breakdown at high angles of attack have been investigated in detail [2,3]. In many configurations the flow separation, which forms the initial stage of vortex formation, is fixed by the sharp leading edge, therefore the main challenge of turbulence models is to correctly produce formation and further development of the vortical flow system along the wing surface. With increased complexity of the configuration such as multiple leading edge angles, variations of edge contours and other devices of flow control, it is impossible to predict the flow behaviour without detailed simulation or very expensive wind tunnel testing.
The complex turbulence fluctuations in the flow field are captured by the underlying turbulence models. Generally, these turbulent fluctuations are represented by the Reynoldsstress tensor in the RANS momentum equation. Different assumptions are used for modeling the Reynolds-stress tensor, which categorizes the type of the turbulence model used in the solver. The widely used Boussinesq assumption relates the stress tensor linearly to the velocity gradients by means of the turbulent viscosity identifying the so-called EVM. In case of the one-equation eddy viscosity turbulence model, one transport equation is used to describe the transport of one scalar (the turbulent viscosity) [4]. However, although the classical RANS models are very efficient in terms of computational time and can be applied for large ranges of computations, they are not capable of predicting the flow in these configurations sufficiently accurate. At higher AoA, the vortex flow pattern further complicates, for example due to the presence of the vortex breakdown, and the numerical simulations often deviate from experimental data, especially in the vortex regions. Different methods are then present in literature in order to overcome the deficiency of the RANS simulations.
Moioli et al. [5] for example aimed to adapt a turbulence model to a specific application of vortex dominated flows. The model terms have been modified in order to achieve better agreement with measured data from experiments. Since the Boussinesq assumption limits the potential accuracy of RANS numerical simulation of vortex flow types, different approaches, such as the Spalart-Allmaras One-Equation Model with Quadratic Constitutive Relation (QCR) and the Reynolds-Stress-Transport (RST) models [4], have been proposed to remedy some of the shortcomings of the linear eddy-viscosity models. However, it is unlikely that a RANS model, even a complex and costly one, will provide the accuracy needed in all varieties of vortical, highly unsteady, flows.
On the other hand, resolving turbulence employing the Direct numerical simulation (DNS) or Large eddy simulation (LES) numerical methods is far too expensive in terms of computational time to apply it on a routine basis. For this reason, instead of modeling the entire turbulent spectrum, it is possible to resolve parts of the spectrum by means of a scale-resolving simulation, using hybrid RANS/LES which becomes an alternative to accurately capture the unsteady characteristics of various scale vortices at slightly higher cost. A promising research contribution in the field of hybrid RANS/LES is given by Zhou et al. [6]. The computations have been performed for the turbulent flow around a delta wing at a low subsonic Mach number and the delayed detached eddy simulation with shear-layer adapted (SLA) subgrid scale model has been applied to predict the vortex breakdown phenomenon. RANS and hybrid RANS-LES computations were carried out for the flow about the VFE-2 delta wing by Peng and Jirásek [7]. The hybrid RANS-LES computation has reproduced the mean flow in a more reasonable pattern than the RANS computation, conducted with the Spalart-Allmaras (SA) model and an Explicit Algebraic Reynolds Stress Model (EARSM), in view of the resolved secondary vortex and the predicted surface pressure. Besides, Cummings and Schütte [8] have presented numerical simulations, performed using RANS, DES, and several DDES turbulence models, of the flow for the VFE-2 delta wing configuration with rounded leading edges. The simulated flow field using SA-DDES has been significantly improved over any of the other hybrid turbulence model simulations and the results showed promise for gaining a fuller understanding of the flow field.
The present work aims to provide a contribution in the field of hybrid RANS/LES numerical methods in order to understand if they are applicable to simulations of a complex delta wing model at transonic flow conditions. The main goal is to provide a significant advancement in the prediction of multiple-delta wing flow for the understanding of the several flow physics phenomena that occur over the aircraft. For this reason, the vortex dominated flow around the triple-delta wing ADS-NA2-W1 aircraft has been investigated comparing URANS with scale-resolving DDES results and experimental data, which consist of integral force and moment coefficients and surface coefficient of pressure over the wing [9]. The Spalart-Allmaras One-Equation Model with corrections for negative turbulent viscosity and Rotation/Curvature (SA-negRC) [10] has been employed to close the RANS equations, whereas in the scale-resolving computations the SAneg-based DDES model has been applied [11]. The transonic regime of Ma ∞ = 0.85 and Re ∞ = 12.53 × 10 6 with = 20 • , 24 • and = 5 • has been selected, which represents realistic conditions for a delta wing and is challenging for the numerical simulation. The DLR TAU-Code release.2019.1.0 flow solver has been employed to perform the simulations [12].

ADS-NA2-W1 test case: geometry and mesh
The ADS-NA2-W1 geometry is a 1:30-scaled version of a generic combat aircraft and it is characterized by a tripledelta wing with three different leading-edge sections, as it can be seen in Fig. 1. The wing thickness is equal to around 0.014L and the wing configuration is equipped with different flat-plate wing planforms including sharp leading edges and sets of corresponding control surfaces. The leading-edge sweep angle of the outer main wing section is 3 = 52.5 • , while the strake section exhibits two different leading-edge sweep angles of 1 = 52.5 • and 2 = 75 • [9]. The nautical labeling is used: leeward starboard ( y > 0 ) and windward portside ( y < 0 ). The dimensionless Cartesian coordinates are furthermore introduced as follows = x∕L , = y∕(b∕2) , = z∕(b∕2) . Figures 1 and 2 show the computational domain employed to investigate the ADS-NA2-W1 geometry. The unstructured mesh consists of about 40 million cells and it is formed by 35 prism layers close to the aircraft and tetrahedral volumes everywhere else. The mesh is symmetric to the plane y = 0 but it is not homogeneous, the cells size varies within the computational domain. The finest cells, whose size is around 0.001 times the characteristic length, = V 1∕3 ≈ 0.001 L , are located close to the leading edge, where the two main vortices are generated, and the mesh refinement roughly follows the vortices in order to capture the turbulent fluctuations along the energy cascade. In order to provide a justification that the grid resolution can assumed to be adequate for the given flow, the number of grid points inside the vortex diameter is analysed at chordwise location = 0.35 . The vortex diameter has been computed from the vorticity distribution x of the vortex, denoted by d or from the distance d of the two extreme of the tangential velocity V , as it has been suggested by Landa et al [13]. Figure 2 shows the qualitative vortex measures of the diameter, indicated by an arrow, and the computational domain at chordwise location = 0.35 . Table 1 summarizes the quantitative data, N = d ∕ and N = d ∕ are the number of grid-points inside the vortex diameter d and d , respectively. Although the cell size slightly increases inside the vortex core along the wing, the ratio of the vortex diameter to the cell size rises due to the vortex expansion.
Since it is challenging and computational expensive to perform meaningful grid-resolution studies for LES-type simulations, three more grid levels are considered for the URANS runs in order to analyze the grid effects.

Mesh convergence study
The mesh convergence study has been performed with four different meshes. The lift and pitching moment coefficient results have been compared. Table 2 summarizes the main mesh characteristics and the results of the aerodynamic coefficients. Figure 3 shows the absolute deviation of the aerodynamic coefficients with respect to the 40 M mesh for which results will be presented.
Even if the 35 M mesh is already adequate to perform URANS, it is further refined to build the 40 M mesh that is used to perform the simulations for both approaches, in order to eliminate any influence of the mesh between URANS and DDES.

Hybrid RANS/LES Method
In the present work the Spalart-Allmaras model that represents a standard RANS closure for aerodynamic applications is employed for the hybrid RANS/LES simulations. The SAneg-DDES model [14] is based on the following oneequation by Spalart-Allmaras [15] for the eddy viscosity where the production term P and the destruction term are This is exactly the original SA-model, except that the length scale d or d in the destruction term is modified. In the SAmodel, d is the distance to the nearest wall (RANS length scale) [16]. In the DDES model, d is replaced with d (hybrid length scale), which is defined as with max = max( x, y, z) , where x, y, z denote the grid spacing in x-, y-, and z-direction, respectively, and f d is a shielding function designed to be unity in the LES region and zero elsewhere [17].
The SAneg is used only in order to improve stability and robustness without changing the (converged) results of the SA model. Equation 1 would be modified and the turbulent eddy viscosity in the momentum and energy equations would be set to zero just in case the kinematic eddy viscosity becomes negative [16].

URANS turbulence model
The Negative Spalart-Allmaras One-Equation Model (SAneg) with corrections for Rotation/Curvature (SA-negRC) is employed to close the RANS equations. The streamline curvature correction was proposed by Shur et al. [18] and it alters the source term with a rotation function, written as follows: with the constants c r1 , c r2 and c r3 calibrated as 1, 12 and 1, respectively, by being multiplied with the production term of the eddy viscosity transport Eq. 1.

Numerical approach
Unsteady simulations have been performed with an implicit dual-time stepping approach, employing a Backward-Euler/ LUSGS implicit smoother. To ensure convergence of the inner iterations in the DDES runs, Cauchy convergence criteria of the variables volume averaged turbulent kinetic energy, maximum eddy viscosity, total vorticity, maximum Mach number and some aerodynamic coefficients with tolerance values of 1e-05 has been used. The CFL number is reduced starting from a large value in order to find the best compromise between speed and stability. The computation of the fluxes have been performed with a central scheme and the matrix dissipation model has been selected. However, in hybrid RANS/LES the artificial dissipation should be reduced in order to prevent excessive damping of the resolved turbulent structures. A (hybrid) low-dissipation low-dispersion discretization scheme (LD2) has been used. It is based on a 2nd-order energy-conserving skew-symmetric convection operator that is combined with a minimal level of 4th-order artificial matrix dissipation for stabilization. A local switch to 1st order scheme by according artificial dissipation is used to stabilize the simulations at shocks locations. The central flux terms employ an additional gradient extrapolation that increases the discretization stencil and is used to reduce the dispersion error of the scheme [19].
The time that a fluid element takes to pass the aircraft has to be taken into account during an unsteady simulations in order to understand how much physical or computational time is necessary to obtain a reliable solution. The so-called convective time unit (CTU) can be computed as follows: where L is the characteristic length (see Fig. 1) and U ∞ the free-stream velocity. (2) Regarding the URANS runs, the selected time step size is equal to 5 × 10 −2 CTU . 10 CTU have been computed before starting the time-averaging in order to overcome the initial transient and 5 flow trough times have been taken into account in order to compute the mean values of the flow proprieties. Approximately 13,000 cpu hours (2-3 days with 8 nodes and 32 cores per node) are required to complete one test case. Figure 4 shows the time history of the pitching moment coefficient C m y and all the described information regarding the time length series.
The DDES runs have been initialized with the URANS results in order to reduce the initial transient, afterwards 3 CTU have been run with the time step size equal to 5 × 10 −3 CTU and then 7 CTU have been computed with t = 2.5 × 10 −4 CTU before starting the time-averaging in order to overcome the initial transient. In the end, 10 overflows have been taken into account in order to compute the mean values of the flow proprieties. Approximately 2 million cpu hours (8-9 days with 200 nodes and 48 cores per node) are needed in the LRZ SuperMuc-NG environment to complete one test case.
Regarding the DDES, in order to fully resolve the convective transport and consequently to capture the flow characteristics accurately, the maximum allowed time step size has been computed. As it has been explained in [20], the chosen time step size, t = 2.5 × 10 −4 CTU , resolves adequately the time scales of the energy containing eddies in the flow of interest. Indeed, the convective CFL number, CFL conv = U t∕ , is lower than unity in each cell of the computation domain, as Fig. 5 shows.

Results and discussion
The transonic regime of Ma ∞ = 0.85 and Re ∞ = 12.53 × 10 6 has been selected. Different URANS simulations have been performed with constant side slip angle = 5 • , focussing on the asymmetry of the turbulent flow, and varying the angle of attack between 12 • < < 28 • . DDES only have been performed for = 12 • , = 20 • and = 24 • due to the high computational costs, but results will be only discussed in detail for the two more challenging, higher angles of attack. The results section contains two main parts. A brief overview of its structure is given as follows: where the instantaneous Q-criterion, the instantaneous x-density gradient iso-surface, the mean surface coefficient of pressure and the mean x-velocity contour are plotted to provide a comparison between the different available data. This section is the core of the present work and aims to provide a significant advancement in the prediction of multiple-delta wing flow for the understanding of the several flow physics phenomena that occur over the aircraft. The flow physics is described, explained and illustrated in detail by dividing the analysis of the unsteady (instantaneous) and the mean flow features. -In Section 4.2 the numerical and the experimental data are compared taking the integral force and moment coefficients into account. In particular, the lift, the rolling and the pitching moment coefficient curves are presented by doing a comparison between URANS, DDES results and experimental data. Moreover, based on the interesting behaviour of force and moment coefficients, several conclusions are drawn and four different flow regimes are identified.

Flow physics analysis
In case of low aspect ratio delta wings, 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 [21]. The suction footprint on the wing surface is mainly caused by the high tangential velocity around the inner vortex core. The so-called vortex lift has a limiting AoA at which the vortex bursts or breaks down. This consists of an abrupt change in the flow topology where the flow decelerates and diverges. The location and mode of breakdown depends on various parameters such as adverse pressure gradients, type of delta wing planforms, angle of attack, sweep angle, interaction with shock waves. The understanding and prediction of the vortex and shock waves, generation and evolution, are of essential importance and are described in the present section. The flow pattern of the ADS-NA2-W1 test case is further complicated by the presence of vortex merging caused by the presence of multiple sweep angles. The side slip angle of = 5 • introduces an asymmetry of the flow and generates two different flow conditions on the two wings (leeward and windward). Moreover, the transonic condition generates a supersonic area over the wing and consequently different shock waves which interact with the vortices and enhance the vortex breakdown. In order to visualize all these phenomena in more detail, the simulation results at = 20 • and = 24 • are visualized graphically. The results are presented by travelling along the wing, from the front to the rear part. All the captured phenomena are analyzed and discussed by focusing on the several physics aspects and doing a comparison between URANS SA-negRC, SAneg-DDES results and experimental data (if available). Figure 6 shows an illustration of the vortices plotting the Q-criterion iso-surface. The iso-surface is colored by the normalized helicity H n , where the positive and negative values are in red and blue, respectively. The rotation sense of a vortex is determined by the sign of the helicity density, so it is possible to differentiate between counter-rotating vortices. This can be used to separate primary from secondary vortices [22], identified with the numbers 1 and 2 in Fig. 6a, b.

Unsteady (instantaneous) Flow Features
Over the wing, the flow undergoes a primary separation at the wing leading edge and subsequently rolls up to form a stable, separation-induced leading-edge vortex. As it can be seen in Fig. 6, two well-distinguished vortices are present on the (leeward) starboard wing and two less-distinguished (more merged) vortices are captured on the (windward) portside. The two vortices (first and second, I and II 1 ) are generated in correspondence with the two increasing sweep angles  The DDES results in Fig. 6b allow for a first, qualitative assessment of the resolution of turbulence in the LES areas. Turbulent fluctuations are clearly visible in the vortices, but the level of resolution is not very high, i.e. only the larger turbulent structures appear to be resolved by the grid (especially on the starboard side). This may be connected with the necessity for further grid-resolution studies, as it will be discussed in Sect. 4.1.3.
As it can be seen in Fig. 6c, d, where the instantaneous iso-surface of the x-direction density gradient ∇ x is shown with flood contour by instantaneous Mach number, several shock waves are present over the aircraft. The interaction phenomenon between leading edge vortices and shock waves that is crucial for the understanding of the flow physics at transonic conditions needs to be assessed in detail, since it could affect the vortex breakdown formation.
The main difference between the two models that is worthy to note is the effect of the highlighted shock wave in Fig. 6d, not shown by the URANS results. Across the shock wave, The other way around, the Mach number, velocity and normal velocity component and total pressure decrease. It means that the vortex core loses velocity and kinetic energy across the shock wave. The interaction between the aforementioned shock wave and the vortex core triggers the breakdown of the first vortex (I.1) on the portside and only the DDES run is able to capture this fundamental phenomenon. Figure 6d shows it also with the reduction of the Mach number behind the shock in combination with the onset of chaotic structures of the vortex breakdown. The breakdown in the transonic regime could be the consequence of a shock/leading-edge vortex interaction. For this reason, the correct prediction of the shock wave location is important. The numerical dissipation in CFD may smear some shock waves making the discontinuity less "sharp" affecting the shock/vortex interaction. Moreover, it is worthy to note that the ability to predict the vortex breakdown influences consequently the prediction of the rolling and the pitching moment, as it will be discussed in Sect. 4.2.
Mean flow features As it can be seen in Fig. 7, where the mean surface coefficient of pressure on the aircraft is shown, four different slice planes have been extracted. The mean surface coefficient of pressure c p distribution along the spanwise direction and the mean x-velocity ( u ) contour at chordwise positions = 0.35, 0.55, 0.75, 0.85 are then plotted in Fig. 8.
In the front part of the aircraft, hybrid RANS/LES considerably improves the results on both wings, as it can be seen in Fig. 8 from the surface coefficient of pressure at the chordwise location = 0.35 . The x-velocity contour plot at the same location shows that only the DDES is able to capture the separation and the reversed flow over the front part of the starboard wing and this phenomenon could explain the great prediction of the experimental data in this area of the wing. However, the simulations are not able to correctly predict the flow physics close to the fuselage (around −0.1 < < 0.1 at this station). The wrong prediction of c p close the fuselage is present more or less all over the aircraft. This might be caused by the presence of two vortices over the fuselage, as the Q-criterion in Fig. 6b shows. Since these vortices are very close to the fuselage surface, they fall into the URANS region that cannot correctly resolve the turbulent flow. A mesh refinement may improve this inaccurate prediction by trying to get them within the LES mode. For this reason it will be of interest to analyze in more detail how the DDES actually behaves, where the hybrid model (DDES) switches from RANS to LES, and this investigation has been conducted in Sect. 4.1.3 for AoA= 24 • .
The DDES captures the secondary vortex formation in particular on the starboard side at the chordwise station = 0.55 , which, however, is not as accurate as desired, and the negative coefficient of pressure is overestimated, as it can be seen in Fig. 8 ( 0.45 < < 0.35 ). Moreover, the same figure shows a better agreement between DDES results and experimental data on the portside, where the secondary vortex is well captured instead. The counter-rotating secondary vortex affects the velocity field. The opposite sign of the vorticity field induces a negative x-velocity. This generates a reduction of the total mean x-velocity, as it can be seen in Fig. 8, close to the leading edge, where areas of low speed flow are visible.
As said before, the separation onset occurs in correspondence with the two sweep angles ( 1 and 3 illustrated in Fig. 1) on the first and third leading edge section. The two generated fully developed vortices (I.1 and II.1) interact with each other in the rear part of the aircraft. In the DDES results the two primary vortices (I.1 and II.1) are still distinguished in Fig. 8 at the chordwise location = 0.75 where the two peaks of axial velocity are located. Taking a look at the surface coefficient of pressure at the same location, the two still separate vortices are confirmed by the presence of the two peaks of negative c p in the experimental results for 0.4 < < 0.8 , even though the two suction footprints are overestimated by the DDES. Instead, the URANS results do not show a wellformed second vortex (II.1) on the starboard side, the emanating shear layer has still to develop and roll up to form it, and fail to predict the flow condition close to the leading edge.
The DDES results on the starboard side deteriorate in the rear part of the wing close to the trailing edge, as it can be seen from the c p curve in Fig. 8 at the chordwise location = 0.85 . The u plot shows two vortex cores, one on top of the other, and a large separation zone close to the leading edge. The reasons for this large separation region will be discussed in the next section and analyzed in future work as well. Besides, there could be an exchange of energy between the vortices in the turbulent structure and the second vortex, which loses kinetic energy (the velocity inside the core decreases), could feeds the first one.
The situation appears differently in the rear part of the portside (windward) wing. The vortices break down within the second half of the wing, as it has been explained in the unsteady flow features analysis. The surface coefficient of pressure in Fig. 8 at the chordwise locations = 0.75 and = 0.85 demonstrates that the vortex breakdown appears in experiment and the simulation results overestimate the suction. The two vortices (I.1 and II.1) do not break down at the same time. The first vortex (I.1) is also the first to burst and then subsequently the second vortex (II.1) breaks down as well. In fact, as it can be seen in the experimental data, the second vortex is still coherent at location = 0.75 , confirmed by the presence of the negative peak of c p for ≈ −0.6 . This suction footprint has vanished at location = 0.85 due to the second vortex breakdown. DDES results in Fig. 8 at the chordwise location = 0.75 predict the first vortex breakdown and demonstrate that it is always accompanied by an expansion of the vortex core and an abrupt reduction of axial (and rotational) velocity. At chordwise location = 0.85 , the DDES results correctly reproduce the breakdown of both the main vortices, whereas the onset of vortex breakdown starts to appear in the URANS ones. Table 3 summarizes the prediction of the breakdown onset position 2 for the two main vortices over the portside (windward) wing by comparing experimental data, URANS and DDES.
The URANS model fails to predict the correct flow physics. On the contrary, the DDES model still highlights some inconsistencies, but it shows a significant step forward for understanding these phenomena.

Alpha = 24 •
The same types of figures discussed in the previous section for the AoA= 20 • are now presented for AoA= 24 • : Figs. 9,  10 and 11 show instantaneous Q-criterion and x-density gradient, c p and field slices, respectively. Fig. 9, two different main vortices (I.1 and II.1) are present on the (leeward) starboard side wing and the burst vortex on the (windward) portside. As it can be seen from the normalized helicity in Fig. 9, on the (leeward) starboard side, the spanwise flow under the primary vortex subsequently separates a second time to form a counter-rotating secondary vortex (I.2) outboard of the primary vortex (I.1). On the other hand, downstream of the incoherent vortex present on the windward portside wing, the flow becomes chaotic and turbulent, as only the hybrid RANS/LES results in Fig. 9b show. An immediate consequence of the vortex breakdown on the windward portside is the increase of the pressure over the wing and, consequently, the reduction of the suction footprint on the wing surface. The prediction of the aforementioned surface coefficient of pressure affects in particular the values of the aerodynamic coefficients which will be analyzed in Sect. 4.2.

Unsteady (instantaneous) Flow Features As illustrated with the instantaneous Q-criterion iso-surface in
On the windward portside, the DDES results in Fig. 9b, d show the chaotic behaviour of the burst vortex and how the shear layer emanating from the leading edge does not roll up to from a coherent vortex. The chaotic structures captured by DDES show how the burst vortex affects also the starboard side wing in the rear part of the aircraft and this could affect the results accuracy.
The iso-surface drawn from ∇ x shows several shock waves over the starboard side of the wing. In this case, they do not affect the vortex evolution inducing a breakdown, but the interaction between shock wave and the first vortex core (I.1) is believed to be the reason for the wrong prediction of the suction footprint on the rear part of the leeward wing. The shock wave captured in Fig. 9d on the rear part of the wing should interact more strongly with the first vortex (I.1), by generating the reduction of velocity and suction over the wing, than with the second one (II.1), as it will be confirmed by analysing the mean flow features. Besides, Fig. 9b shows that the turbulent structures of the first vortex (I.1) tends to dissipate in correspondence of the aforementioned shock wave.
Mean flow features Figure 10 and 11 show the surface coefficient of pressure over the wing and the URANS results mispredict the flow pattern over the (windward) portside in particular close to the wing apex due to high turbulence and chaotic behaviour of the flow. On the contrary, the DDES approach captures the shear layer emanating from the leading edge and chaotically transported downstream over the wing even if the intensity of the suction footprint is slightly overestimated. This phenomenon abruptly changes the aerodynamic coefficients due to the drop in suction footprint behind the transported shear layer and the better simulation of burst vortex over the portside wing in the DDES results generates a significant improvement of the pitching moment coefficient, as it will be discussed in Sect. 4.2.
Regarding the (leeward) starboard side wing, in the front part of the aircraft it is evident how the hybrid method improves the results of the simulation. The u velocity contour plot in Fig. 11 at the chordwise location = 0.35 shows that only the DDES is able to capture the separation and the reversed flow over the front part of the starboard wing. The same consideration has been done for the test case with AoA= 20 • . This is assumed to be a relevant difference between URANS and hybrid RANS/LES methods to explain the different prediction of the suction footprint in the front part of the wing.
The DDES results capture the secondary vortex formation but the coefficient of pressure is still overestimated, as it can be seen in Fig. 11 at the chordwise station = 0.55 close to the leading edge of the starboard side ( > 0.35 ). On the other hand, the URANS results do not capture the secondary vortex formation as accurate. Taking a look at the experimental coefficient of pressure at the same location, it is interesting to note that the first vortex seems to be weaker and the second vortex stronger than the correspondents with the AoA= 20 • .
The trend of the DDES results of the surface c p distribution along the spanwise direction at different chordwise locations is always similar to the experimental data for < 0.55 . This demonstrates that the correct flow pattern is captured by the DDES method. It is useful to note that the slice plane at the chordwise location = 0.55 in Fig. 10a is placed before the second increment of the delta angle. The second vortex (II.1) is then generated, it merges with the first one and the DDES results become less reliable in the rear part of the (leeward) starboard wing, as it can be seen in Fig. 11. Taking the DDES results into account, two different vortices (I.1 and II.1) are still distinguished and interacting with each other at the chordwise location = 0.75 where the two peaks of u and c p are located. The experimental data confirms that at this location two different vortices are still present but the second one should be stronger than the first one, which should lose kinetic energy due to the shock waves captured in Fig. 9d.
The results slightly deteriorates going downstream over the starboard side wing, as it can be seen from the c p curve in Fig. 11 at the chordwise location = 0.85 and as it has been discussed for the test case with AoA= 20 • as well. For both methods (URANS and DDES), the u velocity contour plot shows a single vortex core and a large chaotic separation zone close to the leading edge (larger for the DDES). Actually, the experimental data indicate that the first and the second vortex are still separated over the wing (and not merged) and a strong reduction of suction should be located over the wing close to the leading edge.
Furthermore, as already briefly introduced, the experimental data in Fig. 11 at the chordwise location = 0.85 suggest that the second vortex (II.1, the outer one in Fig. 11) is again stronger than the first one (I.1, the inner one) because it generates a higher peak of suction. The first vortex (I.1) loses energy, while the second one (II.1) gains velocity.
Especially on the (windward) portside wing, the difference between the DDES and the experimental suction footprint, caused by the different prediction of the leading edge shear layer transportation, is almost constant. With the hybrid RANS/LES method important approximations and assumptions have been done on the resolved turbulent scales and this constant gap could be related to the relatively high energy content of the unresolved scales of turbulence. A grid-refinement study could be performed by improving the grid resolution in order to possibly achieve better hybrid RANS/LES results. Moreover, the so-called gray area between the two modes (RANS and LES) could generate a region close to the surface where the shear layer turbulence acts and the flow could not be correctly treated [7,23]. These aspects have been investigated by analysing the employed shielding function and comparing the turbulent eddy viscosity.

Turbulence-related variables and hybrid behaviour
In Fig. 12 the modelled turbulent eddy viscosity, t , is compared at chordwise locations = 0.35, 0.75 between URANS and DDES with = 24 • and = 5 • . The contours of R t = t ∕ is plotted, where is the molecular dynamic viscosity. The URANS SA-model has produced overall higher levels of turbulent eddy viscosity in the vortex cores on the leeward starboard side and in the burst vortex region on the windward portside. In general, the region with large t values in RANS computations corresponds to vortex motions with relatively large turbulence energy generation in relation to strong flow rotation and deformation [7]. 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 and in the burst vortex region, where the turbulent viscosity production is incredibly large and the RANS approach does not provide the required accuracy for the prediction of the flow physics at these flow conditions.
Regarding the hybrid RANS-LES run, Fig. 12 shows the instantaneous hybrid length scale over RANS length scale ( d ∕d ). It illustrates where the DDES approach switches from RANS to LES. The regions close to the wall are resolved by the RANS mode. The DDES approach employs the SGS eddy viscosity in the off-wall region and the DDES R t -contours show that the SGS eddy viscosity is much smaller than its RANS counterpart. The relatively low level of modelled t in LES is associated to the local fine grid resolution required for LES. On the contrary, the region with large t in LES mode indicate strong local flow rotation/deformation and/or coarse grid resolution, inducing usually intensive energy dissipation of resolved large-scale turbulence [7].
One of the most significant reasons behind the discrepancies in c p could be related to the so-called "grey-area" problem rooted in the DDES modelling, for which the resolved turbulence in the shear layer by the LES mode is much less saturated. This is due to the fact that the formation of the vortex is supported by a process of "rolling-up" and "wrapping" of near-wall layers that are modelled by the RANS mode and it generates a rather "stiff" resolved vortex motion by leading to a delayed vortex burst and breakdown. However, in this case, the good results achieved in the front part of the wing and the instantaneous plot in Fig. 12 do not lead directly to this conclusion, and the transition of the shear layer between the RANS and the LES mode should not be the main reason for some discrepancies highlighted in the results. In order to continue this analysis, further studies and articles will focus on the unsteadiness behaviour of the hybrid method results by investigating the Reynolds stresses and the resolved kinetic energy. On the other hand, the instantaneous hybrid length scale over RANS length scale in Fig. 12 suggest to perform a grid-refinement study by improving the grid resolution, in particular, close to the surfaces and to the fuselage in the leeward starboard side even considering that the chaotic structures of the burst vortex captured by DDES in the windward portside may affect the results in the rear part of the starboard side wing. Figure 13 shows the lift coefficient curve, the rolling and the pitching moment coefficient curve, respectively. The experimental data according to [9] are plotted in comparison with the URANS and the DDES results. It shows an interesting behaviour of force and moment coefficients and, based on that, several conclusions could be drawn. In principle, as Fig. 13 shows, four different regimes can be identified as follows:

Aerodynamic coefficients
1. AoA≤ 17.5 • , no vortex breakdown. Within the range below 17.5 • , everything behaves as expected, the lift rises almost linearly with AoA, the pitching moment is negative, stable and decreases slowly, the rolling moment is almost constant. 2. 17.5 • ≤ AoA ≤ 22.25 • , vortex breakdown on the portside (windward) wing. The vortex breaks down on the windward portside and moves upstream from the trailing to the leading edge with the increase of AoA. It generates a gradual lift reduction on the portside that is not very evident in Fig. 13a because meanwhile the lift increases on the leeward starboard side wing (the lift does not rise with the same slope as before). This double effect produces a positive increase of the rolling moment, as it can be seen in Fig. 13b. The pitching moment in Fig. 13c firstly increases due to the breakdown location close to the trailing edge that generates a lift reduction in that specific region and consequently a nose-up pitching of the aircraft. Then, when this phenomenon reaches more or less the x-coordinate of the aerodynamic center (a.c., see Fig. 1 The vortex breaks down on the leeward starboard side and moves upstream from the trailing to the leading edge with increasing AoA . All the comments made for the portside, can be reproduced for the starboard side in order to understand the plots in Fig. 13. The lift reduction in the rear part of the starboard side produces a nose-up pitching moment and a strong reduction of the rolling moment which tends to negative values.
The URANS and DDES results overestimate the experimental lift coefficient but it is worth noting that the DDES are closer to the experimental ones. In particular, Fig. 13a shows a smooth transition between the URANS and DDES points but what happens in between is not well documented in literature and has been discussed in the previous sections. Although the sharp drop of the curve for AoA = 24 • has not been clearly predicted within this study, the DDES results improve the prediction of the lift curve. The rolling and the pitching moment coefficient curves, plotted in Fig. 13b, c, respectively, are particularly interesting with the presence of a non-zero side slip angle. In fact, the integral moments react more sensitive to variations of the flow pattern than the force coefficients. The URANS and DDES results underestimate the experimental rolling moment. A significant deviation occurs between the results of the two approaches at AoA = 20 • , which is generated by the appearance of the vortex breakdown, captured only by the DDES run, on the windward portside wing. The DDES results show a good improvement of the pitching moment values as well. As it can be seen in Fig. 13c, they assume the correct sign of C m y , nose-up pitching for AoA = 20 and nose-down pitching for AoA = 24 . URANS totally mispredicts this coefficient which is caused by the wrong representation of the vortex breakdown over the windward portside wing.
Finally, in order to compare quantitatively the simulation results to the experimental data, the absolute and relative deviation have been computed for the URANS and DDES results with respect to the experimental data. The difference between the experimentally measured quantity and its numerical value from the simulations is the absolute deviation. The simulation results underestimate the experimental ones, if the absolute deviation is positive and vice-versa. Then the relative deviation is the absolute deviation divided by the magnitude of the experimental value. Fig. 14 shows the deviations which confirm that a better prediction of the flow physics affects more the integral moment coefficients than the lift coefficient.
The hybrid method improves the prediction of aerodynamic coefficients and a significant reduction of the Recalling the flow physics analysis in Sect. 4.1, it can be concluded that the still inaccurate prediction of integral moment coefficients (pitching and rolling) is mainly due to the discrepancies between the measured and the simulated vortex breakdown's strength or intensity, where it means the rate of change through it of the surface pressure, and onset position, which affect the suction footprint over the aircraft.

Conclusions
The vortex-dominated flow around the triple-delta wing ADS-NA2-W1 aircraft is investigated in order to achieve a better understanding of the flow physics phenomena that occur over the aircraft particularly at the transonic speed condition. Both URANS and scale-resolving DDES have been employed in order to explore the range of suitability of current CFD methods. The Spalart-Allmaras One-Equation Model with corrections for negative turbulent viscosity and Rotation/Curvature (SA-negRC) has been employed to close the RANS equations, whereas in the scale-resolving computations the SAneg-based DDES model has been applied. The transonic regime of M ∞ = 0.85 and Re ∞ = 12.53 × 10 6 has been selected. Different URANS simulations have been performed with constant side slip angle = 5 • , that emphasizes the asymmetry of the turbulent flow, varying the angle of attack between 12 • < < 28 • . DDES have been performed only for = 20 • and = 24 • due to the high computational costs. At = 20 • two well-distinguished vortices are present on the (leeward) starboard wing and two less-distinguished (more merged) vortices that break down on the rear part of the aircraft are captured on the (windward) portside. Taking instead = 24 • into account, two different vortices are present on the (leeward) starboard side wing and the burst or incoherent vortex fixed on the leading edge apex is located on the (windward) portside.
The sharp leading edge implies that the flow separation takes place well-defined along the entire edge. However, even without the necessity of predicting the flow separation the RANS model fails to predict several flow features correctly. The RANS turbulence model provides excessive eddy viscosity production in the vortex, as it has been discussed in Sect. 4.1.3, with implications on the unburst vortex size, type and velocities. Consequently, the suction peak and the pressure distribution differ from experiments. The breakdown is misrepresented by URANS solutions and, consequently, the surrounding flow and the post-breakdown region are also negatively affected.
Promising improvements have been achieved employing the SAneg-DDES numerical method. The DDES model improves the prediction of the aerodynamic coefficients and provides a significant reduction of the deviation from the experimental results compared with URANS, as it has been shown in Fig. 14. The accuracy of predicting the integral moment coefficients (pitching and rolling) is mainly related with the prediction of the vortex breakdown onset position and strength. It affects the suction footprint over the wings and consequently the surface coefficient of pressure behind the vortex breakdown all over the aircraft. Particularly taking the case with AoA = 20 • into account, the predicted vortex breakdown from the hybrid RANS/LES represents a very important improvement but it is still not strong enough and too close to the trailing edge, as it has been illustrated in Table 3. For this reason, the prediction of vortex breakdown has to be further improved in future. In general, the vortex breakdown phenomenon is of high interest as it changes abruptly the aerodynamic characteristics of a delta wing. A proper prediction of its position and strength is of fundamental importance during the design and development phase of a delta wing based aircraft, and hybrid RANS/LES could surely bring an advantage. Further improvements of the hybrid RANS/LES methods could lead for example to creation of a high-fidelity database in order to conduct aircraft design and could be used for improving cheaper RANS based models.
Good improvements have been obtained with the hybrid method in the front part of the aircraft for both test cases where only the first vortex is present. Only the DDES simulations are able to capture the separation and the reversed flow over the front part of the leeward wing. Based on the pressure gradients over the suction side, secondary vortices have been observed in particular in the DDES flow fields. Some discrepancies between SAneg-DDES results and experimental data are evident in the rear part of the leeward wing where the two generated fully developed vortices merge and interact with each other. Some hypothesis, such as improved grid resolution and grey-area mitigation, have been presented but the actual reasons for this misprediction are not fully understood, yet, and will be investigated in future work.
The flow physics over a delta wing gains further complexity at transonic conditions for the presence of shock waves and the interaction between shocks and vortices has been investigated as well, in particular for AoA = 20 • in Sect. 4.1.1. The interaction between vortex core and shock waves is fundamental for the understanding of the flow physics. It triggers the vortex breakdown on the windward side and only the DDES are able to capture this fundamental phenomenon.
The qualitatively and quantitatively illustrated results in Sect. 4 clearly show that all the computational time spent on DDES has been worth the effort. A significant advancement in the prediction of multiple-delta wing flow has been provided for the understanding of the multiple flow physics phenomena that occur over the aircraft, which to demonstrate was the purpose of the present paper. The hybrid method reveals several physics aspects which have not been seen before (vortex-vortex interaction, shock-vortex interaction) with the demonstrated accuracy.
In future work, the leading-edge vortex structure and shape will be analyzed and described in detail, by focussing in particular on the boundary layer separation process and the secondary vortex formation.
Moreover, other potential hybrid modelling approaches (alternative to the DDES model) have to be taken into account for further studies. For example, in the presented test case, the IDDES method may be better than the DDES model which essentially switches to URANS mode in the wall layer. Since the generation of the turbulence usually start quite early with a delta wing that features a sharp leading edge, it may help to drive the IDDES modelling towards wall-modelled LES by enabling better resolving capabilities in the wall layer downstream.
Finally, in order to overcome the deficiencies due to the Boussinesq linear assumption, the Reynolds-stress models and the Quadratic Constitutive Relation are promising examples of alternative RANS approaches since the standard RANS models are not capable to resolve the flow physics for these flow conditions.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.