Dominik Sedlacek
email: dominik.sedlacek@tum.de.

Michel Visonneau
email: michel.visonneau@ec-nantes.fr

Emmanuel Guilmineau
email: emmanuel.guilmineau@ec-nantes.fr

Jeroen Wackers
email: jeroen.wackers@ec-nantes.fr

Assessment of Hybrid Delta Wing Vortex Flow Investigation -Part I at Subsonic Conditions

The dependencies of numerical simulations regarding grid sensitivity and modelling error for vortex dominated flow need particular investigation. In this study, an extensive investigation on the grid sensitivity is conducted and the influence of different turbulence closures is analyzed. For this purpose, a triple-and a double-delta-wing configuration are considered. In a collaboration between the French National Centre for Scientific Research, the Centrale Nantes Engineering School and the Technical University of Munich, simulations were conducted to identify grid and turbulence closure sensitivities. For the grid sensitivity, a comparison of the global force and moment coefficients and of the local flow in the vortex core is shown. From this, it can be concluded that the local flow is more sensitive to grid adaptation than forces and moments. Furthermore, the double-delta-wing configuration seems to be more affected by grid sensitivity. The comparison of k-𝜔 shear stress model, delayed detached eddy simulation and different grids and background closures indicates a lower influence of grid sensitivity to delayed detached eddy simulations than to k-𝜔 shear stress models. A grid convergence represents an important parameter but for vortex dominated flows, it can only be achieved by a local convergence for the vortex core flow is of primary importance as vortex core flow characteristics strongly determines the accuracy of the moment prediction.

I. Introduction

C omputational Fluid Dynamics (CFD) are widely spread in aerospace research and design of aircraft. However, for complex flow phenomena occurring in vortex dominated flow fields around highly-agile aircraft, the accuracy and the efficiency of numerical methods are promising but there is still room for further enhancement [1]. Hybrid-delta wings present many assets for the aerodynamic design of high agility aircraft. The flight envelope of such high performance aircraft has to include not only an efficient supersonic performance but also acceptable transonic cruise and high agile maneuverability at subsonic speeds. Therefore, high agility aircraft have to cover a wide range of different and partly contrary flight conditions. Such configurations typically feature moderately to highly swept wings with leading edges varying from small radii to sharp edges. This allows the exploitation of vortex related flow phenomena, such as non-linear lift increase and more resilient vortex systems [2].

Hybrid-delta wing configurations combine the advantages of nonslender and slender wings. These wings are split in sections with different leading-edge sweep. This leads to the development of multiple leading-edges vortices (LEV) which can interact with each other and display superior flight mechanics stability characteristics compared to single vortices [«]. In experimental investigations on similar configurations as the ones investigated here, it is shown that two primary vortices separate at double-and triple-delta-wing configurations [«]. The vortex developing in the front highly swept region is called the inboard vortex (IBV) and the vortex developing at the medium swept leading edge is called midboard vortex (MBV). These vortices interact with each other. With increasing angle of attack (AoA), a rising adverse pressure gradient towards the trailing edge results in a destabilization of the vortex up to vortex breakdown [»].

For vortices developing at slender leading-edge sections, the breakdown can be characterized by an abrupt change in the aerodynamic coefficients as well as reversed core flow and strong velocity fluctuations in the vortex core [START_REF] Gursul | Review of Unsteady Vortex Flows over Slender Delta Wings[END_REF]. Therefore,

II. Configurations

In this study, two hybrid-delta-wing configurations are investigated. The model is composed of a generic fuselage and a generic flat-plate wing with a sharp leading edge. The first test case, the so-called NA1_W1, is a triple-delta-wing plan form, shown in Figure 1a. This plan form consists of three areasȷ the Levcon, the first section, with a medium leading-edge sweep of 𝜑 1 = 52.5°, followed by the strake, a highly swept wing section, 𝜑 2 = 75°and the trapezoidal main wing section with a leading-edge sweep of 𝜑 1 = 𝜑 3 = 52.5°. The so-called NA1_W2 configuration consists of a highly swept section 𝜑 2 = 75°and a main wing with a leading-edge sweep of 𝜑 3 = 52.5°. The most important geometrical parameters of the two configurations are summarized in Table 1 [«].

III. Numerical methods

A. TUM

The numerical simulations at TUM are performed with the DLR-TAU-Code, a solver for the Euler or RANS equations on hybrid grids. DDES were also conducted. The idea of the DES, combining RANS and LES, is to model the areas close to the wall and areas with smaller element lengths than the turbulent length scales with a RANS model. The LES model is only used in regions fine enough for LES calculations. This significantly reduces the computing effort and increases the accuracy of the RANS model. The DDES is a modification introduced by Spalart [START_REF] Spalart | A New Version of Detached-eddy Simulation, Resistant to Ambiguous Grid Densities[END_REF][START_REF] Spalart | Detached-Eddy Simulation[END_REF]. It maintains the usage of RANS in the boundary layer even if the grid spacing parallel to the wall is smaller than the boundary layer, so the switching to LES is delayed in boundary layers. TAU uses a dual time step method for unsteady simulations. The physical time step is implicitly discretized. The next physical time step is determined iteratively using explicit pseudo time steps. For the RANS part, the turbulence was modelled by the one equation model of Spalart-Allmaras with rotation correction [START_REF] Shur | Turbulence Modeling in Rotating and Curved Channelsȷ Assessing the Spalart-Shur Correction[END_REF][START_REF] Allmaras | Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model[END_REF][START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]. The unsteady time step for DDES was defined as 𝛿 𝑡 = 2 • 10 -5 𝑠 and the computations converged well within 200 inner iterations per physical time step. As convergence criterion, the Cauchy convergence criterion was used with a control value of Δ = 1 • 10 -6 within 10 iterations for lift, drag and pitching moment coefficients as well as maximum eddy viscosity. A second order central scheme as discretization type was employed. The results of the converged DDES solutions are presented in this study. The shown results represent the averaged values of the flow field. The graphs and figures presented in the following study refer to the mean values of the time-accurate CFD computations.

B. ISIS-CFD

The in-house solver ISIS-CFD developed by the CNRS/ECN authors, also available as a part of the FINE™/Marine computing suite distributed worldwide by Cadence Design Systems, is an incompressible multiphase unsteady Navier-Stokes solver mainly devoted to marine hydrodynamics. It is based on a fully-unstructured (face-based) finite volume discretization with specific functionalities needed for multiphase flows and industrial applications (see Queutey [START_REF] Duvigneau | On the Role Played by Turbulence Closures in Hull Shape Optimization at Model and Full Scale[END_REF] and [START_REF] Cécora | Differential Reynolds-Stress Modeling for Aeronautics[END_REF]. All models are available with wall-function or low-Reynolds near wall formulations. Hybrid RANS/LES turbulence models based on Detached Eddy Simulation (DES SST, DDES SST, IDDES) are also implemented and have been thoroughly validated on automotive flows characterized by large separations, see Guilmineau et al. [START_REF] Guilmineau | Assessment of Hybrid RANS-LES Formulations for Flow Simulation around the Ahmed Body[END_REF] and ships at steady drift [START_REF] Visonneau | Local Flow around a Surface Combatant at Various Static Drift Conditionsȷ The Role Played by Turbulence Closures[END_REF]. Moreover, the solver accepts sliding and overset grids and features an anisotropic adaptive grid refinement (AGR) functionality ( [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF], [START_REF] Mozaffari | Average-based mesh adaptation for hybrid RANS/LES simulation of complex flows[END_REF]) applied to unstructured hexahedral meshes which will be used extensively here.

In this study, two turbulence models are usedȷ the k-𝜔 SST model [START_REF] Menter | Two-equation eddy viscosity turbulence models for engineering applications[END_REF] and the DDES-SST model [START_REF] Gritskevich | Development of DDES and IDDES Formulations for the k-𝜔 Shear Stress Transport Model[END_REF]. The time step used is 𝛿 𝑡 = 1.5 • 10 -4 𝑠 for the k-𝜔 SST model. The details for DDES-SST will be given further in the section devoted to this computation. The averaging time, 𝑡𝑈 ∞ /𝑐 𝑟 varies between «0 and 60 depending on the case. All simulations are performed using the anisotropic automatic grid refinement procedure.

«

IV. Experimental Methods

Experiments were conducted for the validation of the numerical methods. The wind tunnel A of TUM features an open test section with a cross section of 1.80 𝑚 × 2.40 𝑚 and a length of 4.80 𝑚. A maximum velocity of 𝑈 𝑊𝑇 = 65 𝑚/𝑠 can be reached and a turbulence level below 0.4% is provided. For the setting of the angle of attack 𝛼, the angle of sideslip 𝛽 and the roll angle 𝛾, a three component model support is used. The possible range for the angle of attack is

𝛼 = 0 • -40 • .
The experimental methods and results of these measurements are also published in [«, 2«]. For the measurements, a free stream velocity of 𝑈 ∞ = 48 𝑚/𝑠 corresponding to a Reynolds number of 𝑅𝑒 = 3.0 • 10 6 based on the reference length 𝑙 𝑅𝑒 = 1 𝑚 was used. The following will be a short overview of the most important parameters.

A. Force and moments

The aerodynamic forces and moments were acquired with an internal six-component strain-gauge balance. For axial, lateral, and normal forces, the maximum sustainable loads read 900 𝑁, 450 𝑁, 2500 𝑁, respectively. The maximum sustainable moments are 120 𝑁𝑚, 160 𝑁𝑚, 120 𝑁𝑚 for rolling, pitching, and yawing moments, respectively. The measurements were performed with a sampling frequency of 𝑓 𝑚 = 800 𝐻𝑧 for a total acquisition time of 𝑡 𝑚 = 10 𝑠.

Ref. [«] presents the accuracy of the aerodynamic coefficients for the applied test setup with respect to repeatability. The force and moment measurements presented in the following refer to the mean values of the acquired series of measurements per angle of attack.

B. Stereo-Particle Image Velocimetry

Stereo-Particle Image Velocimetry (PIV) is used to resolve the flow field above the wing and to obtain all three velocity components. The PIV system consists of two sCMOS cameras with a resolution of 2560 × 2160 pixels, a double pulsed NdȷYag laser with a maximum power of 325 𝑚𝐽 per pulse and a wave length of 532 𝑛𝑚 and a three axis traversing system next to the wind tunnel test section, see Figure 2

A. Grid sensitivity analysis

General context

In order to assess the respective weights of discretization and modelling errors and to evaluate the minimum requirements in terms of local and global grid density needed to maintain the numerical error under control, this study starts with a grid sensitivity analysis performed by CNRS/ECN. The French authors compare a first grid realized by TUM with a series of grids built by them with the help of the anisotropic adaptive grid refinement implemented in their flow solver. Both wing configurations, i.e. NA1_W1 and NA1_W2, are considered at the same angle of attack, 𝛼 = 16 • , and the computations are carried out with the k-𝜔 SST turbulence model. To build their series of flow-adapted grids, ECN/CNRS uses the procedure described in [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF], which consists in varying systematically the threshold which controls the adaptive grid procedure based on the so-called flux-component Hessian (FCH) criterion computed from Hessians of the pressure and velocity components.

Table 2 provides the values of thresholds, number of cells and faces on the aircraft for both wing configurations for the TUM grid and seven other grids built with the aforementioned grid adaptation. While the global number of cells places the TUM grid between G2 and G1 grids, the local density in the main IBV vortex rather ranks it between G« and G2. Figs. « and » illustrate the local crosswise grid densities for all these grids in the section 𝑥/𝑐 𝑟 = 0.592. For both wings, it is obvious that grids G6 to G» are locally too coarse to describe the main vortices but the last four grids G« to G0 appear to capture ideally the vortical structures of interest.

Forces and moments

Based on the previous grids, Figs. 5a and5b show the grid convergence of the drag coefficient for both wings. The four last adapted grids (in Figs. from right to left) lie in the asymptotic range, like the TUM grid. Therefore, the difference between the computed drag coefficient and the measurements (in blue) is unequivocally due to a modelling error associated with the linear isotropic 𝑘 -𝜔 SST model. The drag and lift are over-estimated for NA1_W1 and under-estimated for NA1_W2, a better agreement being observed for the NA1_W2 wing configuration. It is also reassuring to observe that the results obtained on the non-adapted TUM grid are in perfect agreement with that obtained in case of grid adaptation. Figs. 7 to 10, showing the grid convergence of the individual viscous and pressure components of the drag and lift for both wings, illustrate a well-known result, that the viscous component is significantly less sensitive to the grid discretization than the pressure component.

Local Flow

The previous section established that the drag and lift were not significantly affected by the grid density for this specific angle of attack, whatever the considered wing. This finding is confirmed by Figs. 1» and 15 which show dimensionless axial velocity at section 𝑥/𝑐 𝑟 = 0.593, once again for both wing types. For this angle of attack, a jet-like behavior is observed in the core for both IBV and MBV vortices which is already captured in the coarsest grid G6. The axial velocity is maximum in the core flow of both vortices, with higher values for NA1_W2 than for NA1_W1. The increase of the local grid density leads to higher values of the dimensionless axial velocity in the core of the primary IBV vortex and a better capturing of the secondary MBV vortex. The results obtained on the TUM grid look similar to those obtained on the G» or G« grids. They are clearly not grid converged in the core of the vortex at this specific section. This clearly proves that, contrary to the conclusion one might draw from the previous figures, the local core flow computed on grid G«, for instance, is not yet grid-independent, despite the high local density resulting from the use of an efficient local grid adaptation. Thus, it is important to assess the grid sensitivity on these local flow properties since we have observed that the behavior of the flow in the core of the vortex, notably its longitudinal fluid dynamics stability, will strongly determine the accuracy of the aerodynamic moment prediction. Moreover, to evaluate the ability of turbulence models to capture some of the vortex core flow physics, it is crucial to determine the minimal local grid density and number of points needed in the core of a vortex for which one can consider that the assessment of turbulence models is no more affected by the discretization error. This is the reason why we conclude this grid sensitivity study by a vortex core analysis (VCA). First of all, in order to comply with the method used in the experiments, the core of the vortex is identified as the maximum of the longitudinal vorticity. Since the vorticity and second invariant depend on the local flow gradient, these indicators might be very sensitive to the local grid density. When the vortex is intense and stable, both approaches provide a core location which is in good agreement. However, when the vortex gets more diffused and less stable, the two detection methods might provide a different core location, which is moreover quite different from the experimental core determined on an experimental grid which is, most of the time, significantly coarser. This vortex core analysis should therefore be conducted with caution. Then, the dimensionless axial velocity component is extracted at the core of each IBV and MBV vortices and compared with the available experimental data. The Vortex Core Analysis is shown for both wing types in Figs. 18 and turbulence kinetic energy (TKE) extracted at the core of the IBV vortex. For both wing types, the lateral and vertical positions of the IBV vortex (Figs. 19a and19b) are not sensitive to the local grid density, except very close to the trailing edge, where the vortex starts to destabilize. For NA1_W1, the axial velocity component is very much dependent on the local grid density during the initial progression of the IBV vortex up to 𝑥/𝑐 𝑟 = 0.5 where a larger spreading can be observed and results from the TUM grid depart from the others. The longitudinal vorticity appears to be very well converged all along the vortex trajectory for all the grids generated by automatic adaptation. Once again, the TUM grid appears to be too coarse to capture the strong velocity gradients present in the core of the IBV vortex, particularly in the initial zone of progression. Finally, the turbulence kinetic energy appears reasonably well converged on the series of adapted grids, which is not true for the TUM grid with a factor two between the two evaluations of the turbulence kinetic energy. However, the global shape of evolution of the turbulence kinetic energy -an initial decay up to 𝑥/𝑐 𝑟 = 0.5 and then an increase up to 𝑥/𝑐 𝑟 = 0.85 -is very similar for the two classes of grids. For NA1_W2, the longitudinal velocity component U (Fig. 19c) is much more affected by the grid sensitivity, especially close to the onset of the vortex. The maximum value is increased by 16% and converges only for the finest grids G1 and G0. Particularly, the solution computed on the TUM grid is far from being grid-converged before 𝑥/𝑐 𝑟 = 0.5. The same behavior is observed on the longitudinal vorticity (Fig. 19d) which is even more affected by the local grid density, as expected. Finally, one can also observe a less pronounced grid sensitivity for TKE (19e), although the solution on the TUM grid appears outside of the range of the solution obtained on the grids generated by local adaptation from the onset to 𝑥/𝑐 𝑟 = 0.5. The behavior of U, 𝜔 𝑥 and TKE reveals that the TUM grid is not well adapted to the onset of the IBV vortex and is probably not sufficiently refined in the region close to the leading edge. a significantly improved agreement of the DDES to experimental data. The solver and setup error is comparable due to the DDES SA and the DDES SST results on the TUM grid. The DDES SA is under-predicting and the DDES SST is over-predicting the lift and drag coefficients. The DDES SA exhibits a low agreement for the coefficients of side force, rolling and yawing moment for the NA1_W2 configuration.

Surface pressure distribution and skin friction lines. The surface pressure distributions provide a first assessment of the vortex behavior, interaction and breakdown and represent the connection between aerodynamic stability parameters and flow field patterns. Therefore, such data allow an inference of the origin of specific stability behaviors. Fig. 20 shows the surface pressure distribution and skin-friction lines obtained on the aforementioned combination of turbulence closures and grids for the NA1_W1 configuration. The influence of the grid density on the surface pressure distribution is obvious if one compares the two computations performed by the French authors on the two different grids with k-𝜔 SST.

The differences between the results of the hybrid RANS/LES model on the two grids are significantly alleviated. It is reassuring to notice that the computation performed with TAU on the TUM grid are similar to those made with ISIS-CFD on the same grid. Although the background RANS model is different (k-𝜔 SST for ISIS-CFD and Spalart-Allmaras for TUM), the surface pressure field does not seem to be much affected by that at this level of detail. The most obvious differences are the vortex breakdown in the k-𝜔 SST -TUM grid combination, see Fig. 20b at the windward wing. The surface pressures on the windward main wing are higher compared to the other setups, indicating a significantly more upstream vortex breakdown. This is the reason for the wrong prediction of the rolling moment in Table «. For a more detailed comparison, the surface pressure distributions at different cross-sections are shown in Fig. 22. The first cross section at 𝑥/𝑐 𝑟 = 0.125, see Fig. 22a features a different scale at the y-axis for an improved comparability. The distinct suction peaks of the IBV on each wing side are visible. The lower surface pressure at the windward side predicted by k-𝜔 SST on the TUM mesh indicates an already burst vortex. Compared to the stable IBV in the k-𝜔 SST -Hexpress+AGR, this behavior implies a high grid influence on the k-𝜔 SST closure. This behavior occurs in all shown surface pressure cross sections. The low surface pressures on the windward side result in the high positive rolling moment, shown in Tab. «. The DDES results are all exhibiting a similar suction peak, and are only varying in the development of the secondary vortex occurring close to the leading edge. Additionally, the IBV peak at the leeward side is slightly higher for DDES SA simulations. In the following sections, 0.475 < 𝑥/𝑐 𝑟 < 0.823, the suction peak reduces significantly for all numerical setups. The DDES SST and Hexpress+AGR combination features the highest negative surface pressure and the most pronounced secondary vortex peak. This can be attributed to the higher grid resolution compared to the TUM mesh and the improved prediction of the near-wall flow compared to RANS simulations. The surface pressure distributions at four cross-sections for the NA1_W2 configuration are shown in Fig. 2«. The different setups exhibit a similar behavior in the first slice, see Fig. 2«a. Compared to the NA1_W1 configuration, the k-𝜔 SST combined with the TUM mesh is more aligned with the other setups. The DDES SA and the DDES SST on the TUM mesh exhibit a similar behavior at the windward side, 𝑦/𝑠 𝑙𝑜𝑐 > 0, in all cross-sections. On the leeward side, the DDES SA features a continuously lower suction peak than the DDES SST. At 𝑥/𝑐 𝑟 = 0.592, the suction peak in DDES SA is flattened due to the movement of a fuselage vortex below the wing vortex system. In the other setups, this flattened peak does not occur. This is because the trajectory of the fuselage vortex crosses the trajectory of the 1» IBV more upstream. The differences between the k-𝜔 SST on TUM and Hexpress grid show a high sensitivity to grid resolution, especially, on the windward side and in regions where the vortex system gets more unstable. Fig. 22 (resp. Fig. 2«) shows the surface pressure distribution at four cross-sections for the NA1_W1 (resp. NA1_W2) wing configuration. For both wing forms, the main alteration concerns the windward wing for which the k-𝜔 SST cannot capture the local maximum observed on locally denser grids. The grid dependence is less marked on the leeward wing. Another unexpected observation is the grid sensitivity of the hybrid RANS/LES simulations. For both codes ISIS-CFD and TAU, the results seem to be less affected by the local grid density than to what is observed for the k-𝜔 SST closure. However, on the windward wing, the highest surface pressure maximum is still reached on the flow-adapted grid.

Crosswise axial velocity distribution at x/c r = 0.592. Fig. 2» (resp. Fig. 25) depicts the dimensionless axial velocity distribution at 𝑥/𝑐 𝑟 = 0.592 for the NA1_W1 (resp. NA1_W2) wing configuration. For NA1_W1, for the TUM grid, ISIS-CFD based on the k-𝜔 SST closure predicts a very large wake-type IBV vortex and no MBV vortex on the windward wing at this cross-section. This complete vortex breakdown results in a high surface pressure, a positive rolling moment coefficient as well as an alleviated pitching moment coefficient. On the flow-adapted grid, the solution is totally different with a less pronounced wake-type IBV vortex and a well define jet-type MBV vortex. On the leeward side, the two solutions are very similar. The same trend, but less pronounced, can be observed when the DDES-SST hybrid RANS/LES is usedȷ the wake-type IBV vortex is less diffused on the flow-adapted grid than on the TUM mesh.

The TAU code and its DDES model employed on the TUM grid leads to a dimensionless axial velocity component roughly similar to what is provided by ISIS-CFD with its DDES closure and the same TUM grid. However, the windward IBV vortex appears less diffused and it has a larger wake-type vortex core than what is found by ISIS-CFD on the same grid with DDES-SST. These differences should be attributed to the different background RANS model used by the two teams. Anyway, none of the computations is able to capture the intense jet-like IBV leeward vortex revealed by the TUM experiments. For NA1_W2, the simulation realized by ISIS-CFD on the TUM grid with k-𝜔 SST has more diffused IBV vortices on both windward and leeward wings, underlining the same grid sensitivity. For instance, the windward MBV vortex is not visible on the TUM grid and the leeward wake-type core of the IBV vortex is surrounded by a zone with a more marked overshoot on the axial velocity component. This overshoot in axial as well as crossflow velocities and a high annular vorticity level surrounding the vortex core seem to be the reason for the higher surface pressure in the k-𝜔 SST results. Moreover, the windward MBV vortex is clearly more visible on the flow-adapted grid. On the TUM grid, ISIS-CFD and TAU, using a hybrid RANS/LES model, predict a very similar velocity distribution, with more marked velocity overshoots for both vortices on the leeward side and for the MBV vortex on the windward side, for the TAU solution. This should be attributed to the differences on the RANS background closure. Both DDES models capture the flow characteristics similar to the PIV data.

Vortex core analysis. As discussed in subsubsection V.A.« and emphasized in paragraph V.B.2, the vortex core flow characteristics are crucial for an accurate prediction of the flight stability characteristics. Therefore, a VCA is shown in Figs. [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF] and 27. The location of the vortex core is obtained by the maximum axial vorticity for jet-type vortices and for wake-type vortices. The spanwise and vertical position as well as the non-dimensional axial velocity of the vortex core are presented in Fig. 2» for NA1_W1 configuration. The windward core characteristics are shown in Figs. 26a -26c and the leeward side is presented in Figs. 26d -26f. The y-position of the windward side is predicted well by the DDES simulations. The k-𝜔 SST trajectories with both grids are more off from the experimental data due to the burst vortex and highly transient fluctuations in the vortex core. Therefore, the maximum vorticity is not necessarily located in the center of the circular area with higher vorticity. The vertical position of the IBV core predicted by the k-𝜔 SST closure is located closer to the wing surface compared to experimental and DDES data. Particularly, the DDES SST with the flow-adapted grid and the DDES SA with the TUM mesh exhibit a similar trajectory. The velocities on the windward side indicate a vortex breakdown in the front wing area. All simulations exhibit reversed flow over the whole root chord, except the k-𝜔 SST with the flow-adapted mesh. The MBV exhibits a steady decrease in axial core velocity. This is captured by all setups reasonably well. The leeward side features a stable vortex system in the front region and an experimentally determined breakdown position at 𝑥/𝑐 𝑟 ≈ 0.62. The simulations predict the trajectory in spanwise and vertical direction quite well including the vortex interaction and the MBVs inboard as well as downward movement. However, the DDES SST with the flow-adapted grid is the only setup predicting the vortex position in the rear wing correctly. The other setups predict a downward movement of the IBV. In the front wing, the axial core velocity is well predicted by DDES SA until the vortex breakdown at 𝑥/𝑐 𝑟 ≈ 0.4 indicating a more upstream breakdown compared to the experimental data. The ISIS-CFD solutions show a rather wake-type core flow along the whole root chord. The MBV core velocity is underestimated by all setups, but by the k-𝜔 SST closures the most. The NA1_W2 core flow exhibits a totally different behavior, see Fig. 27. As for the NA1_W1 configuration the windward core characteristics are presented in Figs 27a -27c and the leeward characteristics in Figs 27d -27f. On the windward side the IBV breaks down at 𝑥/𝑐 𝑟 ≈ 0.5 and the MBV at 𝑥/𝑐 𝑟 ≈ 0.65. The windward trajectory in spanwise as well as vertical direction is mostly predicted well. The k-𝜔 SST closures show deficiencies regarding the trajectory after vortex breakdown. The axial core velocity is exceptionally well predicted by the DDES closures, especially the vortex breakdown is in very good agreement with the experimental data. The vortices on the leeward side interact strongly and the experimental data indicate a vortex merging at 𝑥/𝑐 𝑟 ≈ 0.9. The leeward side trajectory is again predicted in good agreement with the experimental results, especially the DDES SST with the TUM mesh is in very good agreement. The DDES SA trajectory diverges on the rear wing and the k-𝜔 SST closures predict a downward movement of the IBV in the rear wing region. The axial core velocity, see Fig. 27f, is significantly underestimated by the k-𝜔 SST closures. The hybrid RANS/LES models exhibit a better agreement with the experimental data. The DDES SA exhibits velocities close to the experimental results above the strake but a steep decrease above the main wing indicating vortex breakdown. This is also shown by the k-𝜔 SST models. The DDES SST exhibits a stable vortex system and the best agreement with the experimental data. The VCA shown in Fig. 26 (resp. Fig. 27) for the NA1_W1 configuration (resp. NA1_W2) exhibit a significant dependency to the turbulence closures and a higher grid sensitivity for the k-𝜔 SST closures. The DDES setups show differences at the leeward side for the NA1_W1 configuration with a more stable IBV in the DDES SA setup. For the NA1_W2 configuration, the velocities in the jet-type IBV are predicted higher by DDES SA compared to DDES SST which can be attributed to the influence of the background RANS closure. The SA with rotation correction leads to a less diffused vortex core and therefore, in most cases to a more stable vortex and higher core velocities.

VI. Conclusion

In computational fluid dynamics (CFD), grid sensitivity and the impact of turbulence models are significant for the accuracy of numerical simulations. In this study, an extensive grid sensitivity analysis is conducted and the influence of modelling error regarding global and local flow characteristics is investigated. The considered configurations are a generic triple-delta-wing and a double-delta-wing configuration. In a collaboration between Technical University of Munich (TUM), French National Centre for Scientific Research (CNRS) and the Centrale Nantes Engineering School (ECN), different grid types, turbulence models and solvers are compared. CNRS and ECN conducted simulations with the ISIS-CFD solver and performed grid adaptation which is based on the so-called flux-component Hessian criterion computed from Hessians of the pressure and velocity components. The grid sensitivity study is conducted on eight grids -seven grids with different resolutions generated by the French authors with an automatic grid refinement (AGR), and a hybrid grid introduced by the German authors. The TUM grid can be ranked between the second and fourth finest CNRS mesh depending on the scale used. The aerodynamic longitudinal stability characteristics resulting from simulations on different grids indicate well converged grids even for the coarser versions. However, the local vortex core flow characteristics exhibit major differences between

Fig. 1

 1 Fig. 1 Wing planforms

  . The cameras were placed up-and downstream of the measurement plane with an angle of 60 • . The sCMOS sensor planes were tilted by Scheimpflug adapters to meet the Scheimpflug criterion, cf. [2»]. The cameras have to be aligned with the laser sheet and were calibrated at an angle of attack of 𝛼 = 16 • . For a change in angle of attack, the cameras and laser sheet can be rotated around the Y-axis of the traversing system. Seeding particles with a diameter of ≈ 1 𝜇𝑚 were infused into the flow. The shown results are an average of the recorded 400 image pairs per cross section with a sampling frequency of 𝑓 𝑚 = 15 𝐻𝑧. These were cross-correlated and corrected with the program Davis. The uncertainties of the mean velocity components were quantified as |𝑢 𝑒𝑟𝑟 /𝑈 ∞ | < 0.06 and |𝑣 𝑒𝑟𝑟 /𝑈 ∞ | = |𝑤 𝑒𝑟𝑟 /𝑈 ∞ | < 0.035 [25].

  Fig. 2 Wind tunnel setup [3]

Fig. 3

 3 Fig. 3 Mesh at x/c r = 0.592 vs the size of mesh for the NA1_W1 geometry.

Fig. 4

 4 Fig. 4 Mesh at x/c r = 0.592 vs the size of mesh for the NA1_W2 geometry.

Fig. 5

 5 Fig. 5 Total drag coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 6

 6 Fig. 6 Total lift coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 7

 7 Fig. 7 Viscous drag coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 8

 8 Fig. 8 Viscous lift coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 9

 9 Fig. 9 Pressure drag coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 10

 10 Fig. 10 Pressure lift coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 11

 11 Fig. 11 Total pitching moment coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

  [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF] which display the lateral Y (resp. vertical Z) position of the core of the IBV vortex, the longitudinal velocity component U, the longitudinal vorticity 𝜔 𝑥 and the

Fig. 12

 12 Fig. 12 Viscous pitching moment coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 13

 13 Fig. 13 Pressure pitching moment coefficient vs the size of mesh for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 15 3D

 15 Fig. 15 3D longitudinal vortices represented by dimensionless Q*=50 surfaces vs the size of mesh for the NA1_W2 geometry for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 16

 16 Fig. 16 Dimensionless axial velocity at x/c r = 0.592 vs the size of mesh for the NA1_W1 geometry for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 17

 17 Fig. 17 Dimensionless axial velocity at x/c r = 0.592 vs the size of mesh for the NA1_W2 geometry for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 18 NA1_W1

 18 Fig. 18 NA1_W1 -Grid sensitivity and Vortex Core Analysis for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 19

 19 Fig. 19 NA1_W2 -Grid sensitivity and Vortex Core Analysis for 𝛼 = 16 • and 𝛽 = 0 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15. Table 3 NA1_W1 -Forces and moments at 𝛼 = 24 • and 𝛽 = 5 • . Model Grid C D C S C L C mx C my C mz Experiments 0.572» 0.0«»7 1.262«1 -0.002» 0.1096 -0.0151 ISIS-CFD k-𝜔 SST Hexpress+AGR 0.5791 0.0»26 1.28»» -0.0177 0.1011 -0.0120 ISIS-CFD k-𝜔 SST TUM 0.5219 0.0»08 1.155» 0.0187 0.0896 -0.0186 ISIS-CFD DDES SST Hexpress+AGR 0.5860 0.0«89 1.«026 -0.026« 0.1100 -0.0110 ISIS-CFD DDES SST TUM 0.5699 0.0«85 1.2720 -0.0161 0.1075 -0.0127 TAU-CFD DDES SA TUM 0.5»8« 0.0«9» 1.2252 -0.010« 0.11«9 -0.0117

1 «

 1 Additionally, the streamlines in all DDES simulations indicate reversed flow close to the leeward trailing edge caused by a flowTable 4 NA1_W2 -Forces and moments at 𝛼 = 24 • and 𝛽 = 5• . shows the same data for the NA1_W2 configuration and illustrates results obtained on the same grids and turbulence models already discussed for the NA1_W1 wing. The results show a similar flow behavior at the windward wing. The vortex breakdown is observable due to pressure increase in the rear wing region and at the main wing. This occurs significantly more upstream in Fig.21b. Similar to the NA1_W1 configuration, the DDES simulations indicate reversed flow in the leeward rear wing area. Additionally, the high suction area in the DDES results indicates a more stable vortex system on the leeward side compared to the RANS results.

Fig. 20 NA1_W1

 20 Fig. 20 NA1_W1 -Surface pressure distribution 𝐶 𝑝 and skin friction lines at 𝛼 = 24 • and 𝛽 = 5 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Fig. 21 NA1_W2

 21 Fig. 21 NA1_W2 -Surface pressure distribution 𝐶 𝑝 and skin friction lines at 𝛼 = 24 • and 𝛽 = 5 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

( a )Fig. 22 NA1_W1

 a22 Fig. 22 NA1_W1 -Surface pressure distribution 𝐶 𝑝 at different cross-sections at 𝛼 = 24 • and 𝛽 = 5 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

( a )

 a Fig. 23 NA1_W2 -Surface pressure distribution 𝐶 𝑝 at different cross sections at 𝛼 = 24 • and 𝛽 = 5 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

  Fig. 24 NA1_W1 -Dimensionless axial velocity distribution at 𝑥/𝑐 𝑟 = 0.592 for 𝛼 = 24 • and 𝛽 = 5 • ; 𝑅𝑒 = 3.0•10 -6 and 𝑀𝑎 = 0.15.

Fig. 25 NA1_W2

 25 Fig. 25 NA1_W2 -Dimensionless axial velocity distribution at 𝑥/𝑐 𝑟 = 0.592 for 𝛼 = 24 • and 𝛽 = 5 • ; 𝑅𝑒 = 3.0•10 -6 and 𝑀𝑎 = 0.15.

Fig. 26

 26 Fig. 26 NA1_W1 -Vortex core analysis by means of spanwise, 𝑦/𝑐 𝑟 , and vertical 𝑧/𝑐 𝑟 core position as well as nondimensional core-flow velocity 𝑢/𝑣 ∞ for 𝛼 = 24 • and 𝛽 = 5 • ; 𝑅𝑒 = 3.0 • 10 -6 and 𝑀𝑎 = 0.15.

Table 1 Geometrical parameters of the NA1_W1 and NA1_W2 configurations.

 1

			NA1_W1 (TC-01-Inc)	NA1_W2 (TC-0»-Inc)
	𝑐 𝑟	[𝑚]	0.802	0.802
	𝑠	[𝑚]	0.417	0.367
	𝑆 𝑅𝑒 𝑓	[𝑚 2 ]	0.329	0.266
	Λ	-	0.15	0.16
	𝑙 𝜇	[𝑚]	0.468	0.426
	𝑙 1 /𝑐 𝑟	-	0.125	-
	𝑙 2 /𝑐 𝑟	-	0.350	0.475
	𝑙 3 /𝑐 𝑟	-	0.475	0.475
	𝜑 1	•	52	-
	𝜑 2	•	75	75
	𝜑 3	•	52	52

  & Visonneau [12], Leroyer & Visonneau [1«]). The method features several sophisticated turbulence modelsȷ apart from the classical two-equation k-𝜖 and k-𝜔 models, the anisotropic two-equation Explicit Algebraic Reynolds Stress Model (EARSM), as well as Reynolds Stress Transport Models, are available, see Deng et al. [1»], Duvigneau & Visonneau

Table 2 Number of nodes and cells for the meshes of both configurations.

 2

				NA1_W1			NA1_W2
			Nb of cells	Nb of faces	Threshold Nb of cells	Nb of faces	Threshold
				on the aircraft			on the aircraft
			×10 6	×10 5		×10 6	×10 5
		G0	87.8	11.9	0.11»»8	7«.»	10.9	0.11»»8
		G1	76.5	11.0	0.120««	61.»	10.0	0.120««
	Hexpress+AGR	G2 G«	«2.8 19.9	8.2 5.6	0.160»» 0.22690	«».0 19.6	7.0 5.5	0.160»» 0.22690
		G»	18.9	5.»	0.2»066	18.9	5.2	0.2»066
		G5	16.7	5.0	0.«2088	16.1	5.0	0.«2088
		G6	1».«	».6	-	15.«	».8	-
	TUM	TUM	52.0	5.7	-	62.5	7.8	-

VII. Acknowledgements

The work performed by the ECN/CNRS authors was granted access to the HPC resources made available by GENCI of IDRIS and CINES under allocations in recurrent annual projects no. 00129. The project performed by the TUM authors is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)grant number BR1511/1«-1. The fruitful cooperation with Airbus Defense and Space is gratefully acknowledged. Furthermore, the authors thank the German Aerospace Center (DLR) for providing the DLR TAU-Code used for the numerical investigations. Moreover, the authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). The authors also highly appreciate the previous wind tunnel tests by Stefan Pfnür.

B. Simulation Sensitivity to Solver, Grid and Turbulence Model

Description of the hybrid RANS/LES models and computational setup

In the following, the simulation results for both considered configurations at 𝛽 = 5 • and 𝛼 = 24 • are compared. The French authors carried out simulations with the k-𝜔 SST and DDES-SST turbulence closures on the grid provided by TUM, named here "TUM" and the grid resulting from the automatic flow adaptation, named here "Hexpress+AGR". The time step used is 𝛿 𝑡 = 5 • 10 -5 𝑠 for the DDES model. The averaging time, 𝑡𝑈 ∞ /𝑐 𝑟 varies between «0 and 60 depending on the case, except for the case 𝛼 = 2»°and 𝛽 = 5°for the NA1_W2 geometry where the averaging time is 6.5. All simulations are performed using the automatic grid refinement procedure. The German authors provided in this part their hybrid RANS/LES simulations computed on their "TUM" grid. The considered hybrid RANS/LES model used the SA with rotation correction as background closure (DDES SA). For unsteady dual time stepping, a time step of 𝛿 𝑡 = 2 • 10 -5 𝑠 is used and the results are averaged within 𝑡𝑈 ∞ /𝑐 𝑟 = 15 for both configurations. The different setups are compared by means of forces and moments, surface pressure distribution, velocity crossflow sections and vortex core characteristics.

Analysis of the computations

Forces and moments. Table « shows the data of the force and moment predictions obtained on these different grids for the NA1_W1 configuration. The first two lines of this table illustrate the moderate but non-negligible influence of the local grid density on the forces and moments. The results obtained with the k-𝜔 SST closure on the flow-adapted grid being closer to the measurements, particularly, for the longitudinal characteristics 𝐶 𝐷 , 𝐶 𝐿 and 𝐶 𝑚𝑦 . In the lateral the grids. Especially, the results for the TUM mesh differs from the results on the AGR grids. Therefore, a thorough and extensive study on grid sensitivity has to include a vortex core flow analysis (VCA) because the behavior in the vortex core strongly determines the accuracy of the moment prediction especially for lateral stability characteristics. This is also shown in the comparison of simulations to identify grid sensitivity and modelling errors.

Here, simulations are performed for both configurations at an angle of attack of 𝛼 = 24 • and 𝛽 = 5 • . The simulations are performed with the two equation model k-𝜔 SST, and a hybrid RANS/LES method used with SST by the French authors and SA by the German authors as background RANS closure. The French authors perform the simulations on a flow-adapted grid as well as on a hybrid TUM mesh. Comparing the results on both grids, the k-𝜔 SST model shows a higher grid sensitivity than the DDES SST. The DDES results exhibit a better agreement to the experimental data and less differences between the meshes. However, the flow-adapted grid shows the highest resolution in the near-wall flow and the best agreement with the experimental data. For the triple-delta-wing configuration, the leeward inboard vortex (IBV) behavior is not captured by the simulations. On the one hand, experiments indicate a stable vortex system with a jet-type core. On the other hand, all simulations show a wake-type vortex core indicating a more upstream vortex breakdown onset. This leads to the previously discussed influence of the vortex core flow characteristics on the moment prediction. Due to the wake-type vortex core on the leeward side, the rolling moment exhibits a right wing up behavior. In the experimental data, the leeward IBV exhibits a jet-type core and therefore, induces a left wing up component to the rolling moment coefficient which is close to zero. The DDES simulations predict a similar flow field with differences regarding the level of diffusion in the vortex core caused by the background RANS closure.

For the doube-delta-wing configuration, flow field characteristics of the DDES are again in better agreement with the experimental data than the RANS model. But the k-𝜔 SST simulation with the flow-adapted grid provided the best agreement with the experimental force and moment coefficients except for the rolling moment. This results from the higher surface pressure coefficient on the windward wing compared to other simulation setups which induces a right wing up component to the rolling moment -𝐶 𝑚𝑥 is alleviated. The good agreement on the longitudinal force and moment coefficients may be associated with an inaccurate prediction of different flow field patterns. This strengthens the statement proposed before that the grid sensitivity study and the comparison of numerical simulation data on vortex dominated flow field with experimental data has to include a VCA along with surface pressure distributions. Consequently, a good agreement to experimental data in the flight mechanics stability characteristics cannot be considered as sufficient to conclude on a good agreement for the vortex core characteristics. Further investigations are needed to identify the reason for the more diffused vortex cores in the numerical results. This vortex core flow behavior results in an inaccurate prediction of the longitudinal and lateral stability characteristics and therefore, is crucial to requirements on the accuracy of the state-of-the-art CFD.


Michel Visonneau

Emmanuel Guilmineau

Jeroen Wackers

Computational Analysis of Subsonic Vortex Interaction on Multi Swept Delta Wing Configurations

This paper is part of the Special Session organized to report progress made by NATO AVT-316 on the topic Vortex Interaction Effects Relevant to Military Air Vehicle Performance. This collaborative research has been initiated to assess the capability of current CFD methods to predict vortex-interaction effects. In the present work, numerical simulations are performed with the ISIS-CFD flow solver by using a linear isotropic k-𝜔 SST model and an hybrid RANS/LES DDES-SST turbulence closure for several angles of attack at symmetric (𝛽 = 0°) and moderate drift (𝛽 = 5°) conditions. For every turbulence closure, the anisotropic adaptive grid refinement functionality is used. A global and local flow assessment is conducted on the basis of experiments performed by TUM.

I. Introduction

Due to their excellent high angle of attack aerodynamic performance, delta or double-delta wings have been implemented in several advanced fighters. The flow over such wings is dominated by a pair of primary vortices that play a significant role in generating an additional lift at high angles of attack compared with a low-sweep-angle wing [1]. High angle of attack for aircraft is frequently encountered during landing and take-off phases, and during combat manoeuvring. The separated flow field is characterized by unsteady and strong vortical flow structures that can interact with various components of the aircraft and among themselves. These complicated flow interactions are the primary cause of most flight dynamic instabilities. Recently, an example of the flow simulation of the interaction between the inner and outer vortices of F-16XL aircraft has been published [START_REF] Boelens | Comparison of Measured and Block Structured Results for the F-16XL Aircraft[END_REF].

A low-aspect-ratio multiple-swept wing fuselage configuration is currently studied in the NATO STO AVT-«16 task group on "Vortex Interaction Effects Relevant to Military Air Vehicle Performance". This collaborative research has been initiated to assess the capability of current CFD (Computational Fluid Dynamics) methods to predict vortex-interaction effects and extend our understanding of vortex-interaction flow physics for these problems through numerical and physical experimentation.

In this paper, we conduct a numerical simulation of the flow past this wing-fuselage combination at Mach number 0.15 for several angles of attack and two side slip angles, by means of RANS simulations and hybrid RANS-LES simulations. A comparison will be made between the numerical results and experimental data.

II. Test case

Figure 1 illustrates the test configurations [«, »]. The first one is a triple delta wing configuration, the so-called NA1_W1, see Figure 1a, with three different consecutive wing sections featuring a varying leading-edge sweep. The inboard and outboard wing sections exhibit a medium leading-edge sweep of 𝜑 1 = 𝜑 3 = 52.5°and the mid-board wing section is highly swept at 𝜑 2 = 75.0°. The second one is a double delta configuration, the so-called NA1_W2, see Figure 1b, with two different consecutive wing sections of varying leading-edge sweeps characterized by a high leading-edge sweep of 𝜑 2 = 75.0°in the inboard wing section and a medium leading-edge sweep of 𝜑 3 = 52.5°in the outboard wing section. Both configurations have the same total length of l tot = 1.160 m and the same root chord of c r = 0.802 m.

Experiments are conducted in the wind tunnel of Technical University of Munich (TUM) for a speed of 𝑉 ∞ = 51.97 m/s, which leads to a Reynolds number, based on c r , of 2.«6 × 10 6 and a Mach number of 0.15. Two side slip angles, 𝛽, are investigatedȷ 0°and 5°.

III. Flow Solver

The in-house solver ISIS-CFD developed by the authors, also available as a part of the FINE™/Marine computing suite worldwide distributed by Cadence Design Systems, is an incompressible multiphase unsteady Reynolds-averaged Navier-Stokes (URANS) solver mainly devoted to marine hydrodynamics. It is based on a fully-unstructured (face-based) finite volume discretization with specific functionalities needed for multiphase flows and industrial applications (see [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF] [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF], [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF] [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF]). The method features several sophisticated turbulence modelsȷ apart from the classical two-equation k-𝜖 and k-𝜔 models, the anisotropic two-equation Explicit Algebraic Reynolds Stress Model (EARSM), as well as Reynolds Stress Transport Models, are available, see [START_REF] Deng | Comparison of explicit algebraic stress models and second-order turbulence closures for steady flows around ships[END_REF] [START_REF] Deng | Comparison of explicit algebraic stress models and second-order turbulence closures for steady flows around ships[END_REF], Duvigneau & Visonneau (200«) [START_REF] Duvigneau | On the Role Played by Turbulence Closures in Hull Shape Optimization at Model and Full Scale[END_REF] and [START_REF] Cécora | Differential Reynolds-Stress Modeling for Aeronautics[END_REF]. All models are available with wall-function or low-Reynolds near wall formulations. Hybrid RANS/LES turbulence models based on Detached Eddy Simulation (DES-SST, DDES-SST, IDDES) are also implemented and have been thoroughly validated on automotive flows characterized by large separations, see [START_REF] Guilmineau | Assessment of Hybrid RANS-LES Formulations for Flow Simulation around the Ahmed Body[END_REF] [START_REF] Guilmineau | Assessment of Hybrid RANS-LES Formulations for Flow Simulation around the Ahmed Body[END_REF] and ships at steady drift [START_REF] Visonneau | Local Flow around a Surface Combatant at Various Static Drift Conditionsȷ The Role Played by Turbulence Closures[END_REF]. Moreover, the solver accepts sliding and overset grids and features an anisotropic adaptive grid refinement functionality ( [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF], [1«]) applied to unstructured hexahedral meshes which will be used extensively in this study.

IV. Anisotropic adaptive grid refinement

The grids employed in this study are generated with the help of the anisotropic adaptive grid functionality embedded in our flow solver ISIS-CFD, which is an essential element for complex-flow simulation due to the difficulty of creating high-resolution meshes for complex flows by hand. The choice of where to refine the grid is handled by the refinement criterion, formulated in the metric context [1», 15]. The criterion is a symmetric 3 × 3 tensor field computed from the flow, which indicates the desired local cell sizes in all directions. The actual mesh is then refined to fit this specification as closely as possible. The metric approach is very flexible; in this particular study, we use as refinement criteria, the so-called flux-component Hessian based on the second derivatives of the velocity fluxes and pressure. For all turbulence closures, computations are conducted in an unsteady way and the grid adaptation is based on averaged refinement criteria following the methodology described in [ 1«]. For resolved turbulence, mesh refinement is no longer only a means to reduce the numerical errorȷ the mesh resolution and the turbulence model influence each other and have to be mutually adjusted to obtain optimal results. For computations based on statistical (resp. hybrid RANS/LES) turbulence models, the grid is adapted every 10 (resp. 10) time steps in this specific study.

Figure 2 presents the converged flow-adapted grids used for the computations performed around the NA1_W1 with the k-𝜔 SST and hybrid RANS/LES DDES-SST closures at an angle of attack of 2»°for the symmetric and drift conditions. Only the cross-section x/c r =0.592 is shown in this figure. For the symmetric (resp. drift) flow condition, the global flow-adapted grid for the k-𝜔 SST simulation is comprised of 5».5M (resp. 89.6M) cells, 9.6 10 5 (resp. 17.1 10 5 ) faces on the wing surface, while the simulation realized with DDES-SST is comprised of 72.8M (resp. 97.9M) cells, 12.» 10 5 (resp. 28.» 10 5 ) faces on the wing surface. The typical size of cells inside the core of the main vortex is around 0.69mm to be compared with the local Taylor scale of 2.9 mm.

V. Results

A. Symmetric Free-stream Conditions, 𝛽 = 0°1

. Forces and moment Figure « shows the evolution of the drag coefficient, C D , the lift coefficient, C L , and the pitching moment coefficient, C my , versus the angle of attack, 𝛼, for the NA1_W1 and NA1_W2 configurations at 𝛽 = 0°. For the configuration NA1_W1, see Figure «a, the numerical results are in very good agreement with the experimental data, except for 𝛼 = «2°. It is a bit disappointing to see that the hybrid RANS/LES DDES-SST closure does not bring any significant improvement at this high angle of attack. The shape variation of the pitching moment curve around the angle of attack 12°, due to the onset of vortex breakdown at the trailing edge, is well predicted with the k-𝜔 SST model, which is very satisfactory. For the configuration NA1_W2, see Figure «b, the vortex breakdown, observed at 𝛼 = 2«°in the experiments, leads to a local reduction of the lift coefficient. This local reduction of the lift is also predicted by the numerical simulations but at a smaller angle, 𝛼 = 19°with the k-𝜔 SST 199» version (not shown here) and 𝛼 = 20°with the k-𝜔 SST 200« version. The same behaviour is observed on the pitch moment C my whose abrupt increase occurring in the experiments at 𝛼 = 2«°, is predicted too early at 𝛼 = 20°. For higher angles of attack, the lift predicted by 𝑘 -𝜔 SST continues to grow without being able to represent the plateau observed in the experiments from 𝛼 = 28°. Although the hybrid RANS-LES closure does not provide significantly better simulation in this range of moderate angles of attack, i.e. for 𝛼 comprised between 20°and 2»°, it is reassuring to see that its prediction of drag, lift and pitch moment at a high angle of attack 𝛼 = «2°, is now in perfect agreement with TUM measurements.

Angle of attack 𝛼=16°F

igure » presents the vortical structures obtained with both turbulence models for the angle of attack 𝛼 = 16°and both wing configurations. These figures show the isosurface of the dimensionless second invariant, 𝑄 * = 50, the blue translucent surface, and the zero-isosurface of the Cartesian axial component of the velocity, the red surface, which indicates a reverse flow. For NA1_W1, the flow is characterized by two primary leading-edge vortices. The first vortex develops at the most inboard wing section and is called inboard vortex (IBV). The second vortex develops at the kink from the highly swept to the medium swept wing section at x/c r = 0.»75 and is called mid-board vortex (MBV). All « these structures are predicted by every turbulence model. However, the size of these structures is thinner with the DDES-SST model, indicating more intense vortices. For the second wing configuration NA1_W2, k-𝜔 SST predicts a similar unique IBV vortex while DDES-SST simulation provides a collection of peripheral vortices spiraling around a thinner vortex. The interaction with the second vortex MBV, originating from the medium swept wing, leads to a local zone of reverse flow as indicated by the red region which occurs more upstream than in the NA1_W1 configuration for both turbulence models. However, for DDES-SST, the detached layer simulated by k-𝜔 SST is replaced by a collection of individual thinner longitudinal vortices originating from the leading edge of the medium swept wing (Fig. »d). Fig. 5 shows the surface pressure distribution and skin friction lines for both wings and turbulence closures. NA1_W1 has a lower surface pressure zone on the medium swept wing for 𝑘 -𝜔 SST which is not present in the DDES-SST simulation, indicating a more intense and stable longitudinal vortex, as expected from a DDES-SST turbulence closure. For NA1_W2, no major differences between the two simulations are visible in the surface pressure distribution and skin friction lines.

Figure 6 shows the dimensionless axial velocity component at the cross-section x/c r =0.592 for both wing profiles at an angle of attack 𝛼 = 16°. Results obtained with k-𝜔 SST and DDES-SST turbulence closures are compared with TUM experiments. First of all, the IBV vortex in the experiments has a wake-type for NA1_W1 and jet-type for NA1_W2. For NA1_W1, this characteristic is not captured by k-𝜔 SST which does not simulate the strong decrease of the longitudinal velocity, while DDES-SST provides a very accurate solution which captures this characteristic evolution in the core of the main vortex. The MBV vortex is also more intense with DDES-SST and in better agreement with TUM experiments. DDES-SST performs also better for NA1_W2 since it gives a stronger axial velocity component than 𝑘 -𝜔 SST in the core of the IBV, again in good agreement with TUM experiments. However, the axial velocity component in the MBV vortex is overestimated by DDES-SST.

For both wings, Figure 7 shows the dimensionless axial velocity, 𝑢/𝑈 ∞ , extracted at the cores of the inboard (IBV) and mid-board (MBV) vortices for the angle of attack 𝛼 = 16°. As mentioned previously, for NA1_W1, in the experiments, the vortex core in the IBV first shows a jet-type core flow with 𝑢/𝑈 ∞ ≃ 2. After x/c r ≃ 0.«6, the axial velocity decreases continuously until it reaches 0.1 at x/c r ≃ 0.65. Then the velocity increases slowly as it approaches the trailing edge of the wing. This flow change from a jet-type core to a wake-type core flow is well predicted by the DDES-SST turbulence model. With the k-𝜔 SST turbulence model, the high level of axial velocity at the onset of the vortex is not predicted. For the configuration NA1_W2, a jet-type core flow is observed for the IBV in the experiments for which an almost constant level of 𝑢/𝑈 ∞ is obtained with a slight tendency to increase. With the k-𝜔 SST turbulence model, the axial velocity in the core of the IBV is in agreement with the experimental data at the onset of the vortex but the velocity decreases continuously until x/c r ≃ 0.82, while DDES-SST is able to maintain a high level core flow plateau as shown in the TUM measurements.

Angle of attack 𝛼=24°F

igure 8 presents the vortical structures obtained by the two turbulence models for the angle of attack 𝛼 = 2»°and both wing configurations. For NA1_W1, both IBV and MBV vortices are present with a reverse flow in the core of IBV starting earlier than for the previous angle of attack. Both turbulence closures predict the existence of a third vortex parallel to the fuselage during the first part of its path and then deviated towards the outer part of the wing, leading to a potentially more complex vortex interaction than for the lower angle of attack. For NA1_W2, the two main vortices IBV and MBV are present with a reverse flow starting later. Moreover, the third vortex emanating close to the fuselage is also found for this wing but its path remains parallel to the fuselage. Fig. 9 shows the surface pressure distribution and skin friction lines for both wings and turbulence closures at this angle of attack. For NA1_W1, the blue zone associated with the local reverse flow is larger for DDES-SST than with k-𝜔 SST. For NA1_W2, this low pressure zone is less extended and, contrary to the previous wing, DDES-SST provides a slightly reduced low surface pressure region, compared to k-𝜔 SST.

Figure 10 shows the axial velocity component at the cross-section x/c r =0.592 for both wing profiles at an angle of attack 𝛼 = 2»°. For NA1_W1, as observed previously for 𝛼 = 16°, the IBV vortex has a marked wake-type illustrated by the large region of very low axial velocity in its core at this cross-section. k-𝜔 SST predicts qualitatively this trend but still over-estimates the axial velocity in the core of IBV. On the contrary, DDES-SST captures the right almost negative value of the axial velocity in the core and yields a more confined vortex in better agreement with TUM experiments. For NA1_W2, the IBV vortex has still a marked jet-type and DDES-SST is again the only model able to capture accurately the overshoot on the core axial velocity component, 𝑘 -𝜔 SST providing only a fair qualitative trend.

Figure 11 shows the dimensionless axial velocity, 𝑢/𝑈 ∞ , in the inboard vortex (IBV) and the mid-board vortex (MBV) cores for the angle of attack 𝛼 = 2»°. For NA1_W1, both models under-predict the core velocity at the onset but DDES-SST is the only one able to capture the low velocity, close to zero, observed from the middle of the wing. For NA1_W2, the axial velocity is overestimated by DDES-SST at the onset of IBV but the rapid decay of the core velocity around x/c r =0.6 is better captured by this hybrid RANS/LES model.

Angle of attack 𝛼=32°F

igure 12 presents the vortical structures obtained by the two turbulence models for the angle of attack 𝛼 = «2°and both wing configurations. For both turbulence closures, the main IBV vortex exhibits a core fully contaminated with a reverse flow from its onset to its destabilisation for NA1_W1. Moreover, the MBV vortex seems to disappear quickly, at least as shown by the iso-Q* surface although we will see further from the cross-sectional views that it is still present. Two vortices emanate from the fuselage, instead of one for the previous angle of attack. These vortices are deviated towards the outer part of the wing at a location situated more upwind compared to the previous angle of attack but their interaction with the IBV is not visible in the cross-sectional views which will be shown next. For the NA1_W2 configuration, the global topology is very similar but the reverse flow in the core of IBV appears only at some distance from the onset, closer to the leading edge for DDES-SST than for k-𝜔 SST.

Fig. 1« shows the surface pressure distribution and skin friction lines for both wings and turbulence closures. As expected, the blue low surface pressure region is more extended for NA1_W1 than for NA1_W2. Figure 1» shows the axial velocity component at the cross-section x/c r =0.592. For this high angle of attack, TUM experiments reveal that both wing profiles have a wake-type IBV vortex and slightly different longitudinal velocity distribution in the MBV vortex core at this cross-section. For NA1_W1 and NA1_W2, the radius of the low-velocity IBV core is under-estimated by both turbulence closures at this location, DDES-SST being the model which provides the thinnest core. For NA1_W2, k-𝜔 SST predicts a significantly higher axial velocity around the low-velocity core, a result which is not confirmed by the experiments. For the MBV vortex, 𝑘 -𝜔 SST predicts a jet type core while the prediction of DDES-SST, with a weaker longitudinal velocity appears to be in better agreement with TUM experiments. Figure 15 shows the dimensionless axial velocity, 𝑢/𝑈 ∞ , at the core of the inboard vortex (IBV) and the mid-board vortex (MBV) for the angle of attack 𝛼 = «2°. For NA1_W1, both k-𝜔 SST and DDES-SST are in very good agreement with TUM experiments since they predict the correct low and negative axial velocity plateau from the onset and during its progression to the trailing edge. For NA1_W2, the longitudinal evolution of the axial component of the velocity is characterised by a rapid decay occurring for x/c r between 0.» and 0.5. DDES-SST is able to capture the rapid decay,

Forces and moments

Figure 16 presents the evolution of the forces and moments coefficients versus the angle of attack for the configuration NA1_W1. For the longitudinal aerodynamic coefficients, see Figure 16a, the shape of the curves is comparable to that obtained for 𝛽 = 0°with however, some minor differences. The main difference is the shift of the angle where the vortex breakdown occurs. In the experiments, this phenomenon occurs at 𝛼 ≃ 15°with drift instead of 12°without drift. For this configuration, the numerical results obtained with the k-𝜔 SST and DDES-SST turbulence models are in good agreement with the experimental data for the drag and lift, the hybrid RANS/LES closure providing a slightly improved prediction of the pitch moment over the whole range of simulations. The agreement between measurement and simulation is less satisfactory for the lateral aerodynamic coefficients, see Figure 16b. On one hand, the side-slip coefficient C s and C mz coefficient are accurately predicted 𝑘 -𝜔 SST and DDES-SST, with a slight improvement brought by DDES-SST for the highest computed angle of attack, 𝛼=2»°. On the other hand, the numerical prediction of the roll moment coefficient C mx is inaccurate for 𝑘 -𝜔 SST which is not able to capture the decay and the change of slope occurring around 16°. Hybrid RANS/LES computations look more promising at 𝛼=16°with the right negative value and the sudden growth at 𝛼 = 18°where the roll moment coefficient is close to zero, in good agreement with the measurements. But for the higher angles of attack, i.e. Figure 17 presents the evolution of the forces and moments coefficients versus the angle of attack for the configuration NA1_W2. For the longitudinal aerodynamic coefficients, see Figure 17a, only three angles of attack are computed with both k-𝜔 SST and DDES-SST and the agreement is quite good on the three longitudinal force and moment coefficients for the first two angles of attack, i.e. 16°and 2»°. At «2°, the k-𝜔 model reaches its limits and provides too high values of drag, lift and pitch moment coefficients. As noticed for the NA_W1 wing profile, the situation is more complex for the lateral force and moment coefficients. At 16°, both turbulence models provide evaluations of C S and C mz in good agreement with TUM measurements. For higher angles of attack, k-𝜔 SST is unable to capture the monotonous decay of C mz up to «2°. Its prediction of C S is still good at 2»°, but strongly overestimates its growth at «2°. As noticed for the previous wing profile, the behaviour of the roll moment, with its change of sign around 20°, is more difficult to capture. Both k-𝜔 SST and DDES-SST predict a positive roll moment coefficient at a lower angle of attack of 16°but for higher angles of attack, i.e 2»°and «2°, the roll moment coefficient gets negative and decays, with values which are opposite to the experiments.

Angle of attack 𝛼=16°F

igs. 18 and 19 show the vortical structures, surface pressure distributions and skin friction lines for NA1_W1 and NA1_W2 at 16°. For both leeward and windward wings, the IBV vortex is thinner, more intense and more stable when computed with DDES-SST, leading to a local increase of the surface pressure under the trajectory of this vortex. Both wing configurations are concerned by this observation. Figure 20 shows the crosswise dimensionless axial velocity component at x/c r = 0.592 for an angle of attack of 16°. For the NA1_W1 wing (see Figs. 20a, 20c and20e), the TUM experiments reveal a strong asymmetry concerning the nature of the vortices developing on each wing. On the leeward wing, both IBV and MBV are jet-type vortices with a strong overshoot of axial velocity in the core of both vortices. On the windward wing, the IBV vortex is of wake type with a core axial velocity component close to zero at this section while the MBV vortex keeps its jet-type characteristic. Concerning the simulations, none of the turbulence closures is able to reproduce the right behaviour in the core of the IBV vortex on the leeward wing, 𝑘 -𝜔 SST providing a somewhat higher velocity than DDES-SST, but in both cases smaller than the measurements at this cross-section. Both simulations show two vortices very close to each other with opposite axial velocity distribution, while the TUM experiments show only one jet-type IBV vortex. On the windward wing, DDES-SST is the only model able to capture the strong decrease of the core axial velocity in the IBV vortex, k-𝜔 SST yielding too strong axial velocities. Moreover, DDES-SST computations of the MBV vortex is in better agreement than what is obtained with k-𝜔 SST. For the NA1_W2 (see Figs. 20b, 20d and20f), contrary to the previous wing configuration, the TUM experiments reveal more symmetric vortices developing on each wing. On the leeward and windward wings, both IBV and MBV are jet-type vortices with similar overshoots of axial velocity in the core of both vortices on leeward and windward sides. It is therefore interesting to see which turbulence closure is able to capture this different trend compared to NA1_W1, which will have a strong impact on the roll moment. Unfortunately, once again, there is no obvious winner. Although DDES-SST clearly captures the right overshoot in both IBV and MBV vortices on the leeward wing, it fails to reproduce the symmetric behaviour on the windward wing by predicting too weak an IBV vortex with an overshoot of 1.5 instead of the experimental of 2.5. On the other hand, the solution proposed by k-𝜔 SST is more symmetric but, on both leeward and windward wings, the overshoots in the IBV and MBV vortex core are underestimated by a factor 2 (1.2 instead of 2.5).

Figure 21 (resp. 22) shows the dimensionless axial velocity, 𝑢/𝑈 ∞ , in the inboard vortex (IBV) and mid-board vortex (MBV) cores for the angle of attack 𝛼 = 16°, on the leeward (resp. windward) wings. For NA1_W1, the leeward wing is characterised by an intense jet-type IBV vortex with a plateau of high velocity (around 2) from its onset and up to the trailing edge. 𝑘 -𝜔 SST predicts a plateau evolution but with a significantly lower core axial velocity around 1.0 while 1« the core velocity predicted by DDES-SST, decreases gradually from an initial accurate value of 2.0 at the onset to a lower value of 0.5 all along the second half of its trajectory. On the windward wing, TUM experiments show a rapid decay of the IBV core velocity which is very well predicted by DDES-SST while k-𝜔 SST does not simulate this decay. For NA1_W2, he same high core velocity plateau is observed on the leeward wing, but now DDES-SST is able to capture it and provides a result clearly better than k-𝜔 SST. On the windward wing, TUM experiments reveal the same high velocity plateau behaviour with a overshoot of from 2 to 2.5 during the main part of the progression while neither k-𝜔 SST not DDES-SST are able to maintain this high velocity plateau, although DDES-SST captures a higher core velocity close to the onset.

Angle of attack 𝛼=24°F

igs. 2« and 2» show the vortical structures, surface pressure distributions and skin friction lines for NA1_W1 and NA1_W2 at 2»°. For NA1_W1, DDES-SST predicts more intense vortical structures on both leeward and windward wings but the surface pressure and skin friction lines remain globally unchanged by the turbulence closures. For this angle of attack, a flow reversal is observed in the IBV vortex all along its trajectory on the windward wing while flow reversal occurs only close to the trailing edge of the leeward wing for both turbulence closures.

Figure 25 shows the crosswise dimensionless axial velocity component at x/c r = 0.592 for an angle of attack of 2»°. For the NA1_W1 wing (see Figs. 25a, 25c and25d), the TUM experiments reveal a stronger asymmetry than what was observed previously at 16°concerning the nature of the vortices developing on each wing. On the leeward wing, both IBV and MBV vortices are clearly jet-type vortices with a strong overshoot of axial velocity in the core of both 1» vortices (around 2.5 for IBV and 1.5 for MBV). On the windward wing as for 16°, the IBV vortex is of wake type with a core axial velocity component close to zero at this section while the MBV vortex keeps its jet-type characteristic.

Concerning the simulations at this higher angle of attack, none of the turbulence closures is able to reproduce the right behaviour in the core of the IBV vortex on the leeward wing, 𝑘 -𝜔 SST providing a very similar velocity as DDES-SST, in both cases smaller than the measurements at this cross-section. As observed for the lower angle of attack, one can notice the existence of two zones very close to each other corresponding to two co-rotating vortices, the upper one with a velocity overshoot and the lower one with a velocity undershoot. Such a complex vortex-vortex interaction is not confirmed by TUM measurements at this location which show a unique jet-type vortex. On the windward wing, both k-𝜔 SST and DDES-SST models are able to capture the strong decrease of the core axial velocity in the IBV vortex. Moreover, both models slightly overpredict the core velocity in the MBV vortex.

For the NA1_W2 (see Figs. 25b and25e), contrary to what is observed at 16°, the TUM experiments reveal now a strongly asymmetric situation characterised by a very intense jet-type (resp. wake-type) vortices on the leeward (resp. windward) wing, the MBV vortex keeping however its jet-type configuration. 𝑘 -𝜔 SST simulations seem to capture the right trend described in the measurements, but the overshoot in the leeward IBV vortex is strongly under-estimated (1.5 instead of «.0) and the negative axial component present in the core of the windward IBV vortex at this location is not reproduced by 𝑘 -𝜔 SST which predicts a value around 0.1 instead of -0.2. Moreover, the IBV vortex is surrounded by a region of high axial velocity in the k-𝜔 simulation, which is visible in the experiments.

Figure 26 (resp. 27) shows the dimensionless axial velocity, 𝑢/𝑈 ∞ , in the inboard vortex (IBV) and mid-board For NA1_W2, on the leeward wing, k-𝜔 SST predicts a jet-type IBV vortex but the axial core velocity is underestimated (around 2.« instead of «.0) but the trend along the wing is reasonably well captured with a similar high axial velocity plateau which is maintained up to x/c r = 0.6. On the windward wing, k-𝜔 SST underestimates the core axial velocity close to the onset (2.6 instead of «.2). Then, it appears to produce reasonably good results since the vortex breakdown is nicely reproduced for x/c r around 0.» and 0.6, and then followed by a region of zero core velocity up to the trailing edge.

DDES-SST -Analysis of the behaviour of the roll moment for 𝛼 between 16°and 24°for NA1_W1

In Fig. 16b, we have noticed the irregular evolution of the roll moment between 16°and 2»°. The roll moment computed by DDES-SST is negative like in the experiments at 16°. But at 20°, the roll moment keeps decreasing without being able to stay on the experimental plateau characterised by measured values of C mx close to zero. This decrease of C mx continues at 2»°with an even more negative roll moment coefficient in complete contradiction with the plateau observed in the measurements. To try to relate this wrong roll moment evolution to the local flow characteristics, it is interesting to analyze the local flow behaviour specifically for these three angles of attack and compare with the solution provided by k-𝜔 SST at 18°(less negative and therefore closer to the measurements) to understand what is missing at 20°and 2»°. Fig. 28 shows the axial velocity component at the usual cross-section x/c r = 0.592 for all the above-mentioned angles of attack. All the angles of attack are computed with DDES-SST except 18°which is computed with k-𝜔 SST. The velocity defect in the core of the IBV vortex on the windward wing increases with the incidence, which is the expected behaviour. But, for the leeward wing, one observes a non-monotonous evolution of the core axial velocity in the IBV vortex. The axial velocity grows from 16°to 18°and then decreases at 20°to a value very similar to the one observed at 16°and continued its decay at 2»°. The evolution observed at this cross-section x/c r = 0.592 is actually confirmed all along the main part of IBV trajectory on the leeward wing upwind of this specific cross-section up to a specific cross-section which depends on the incidence. It is also worthwhile to mention that, on the leeward wing, the simulations do not exhibit a single IBV vortex but, instead, two vortices close to each other characterised by different core flow characteristics, a wake-type vortex close to the wall and a jet-type vortex on top of it. This two-vortices structure is maintained for these four different angles of attack, although less marked at 18°obtained with k-𝜔 SST. This is in contradiction with the unique clear and coherent IBV vortex depicted in the measurements on the leeward wing.

In order to understand the origin of these two vortices observed the computed IBV vortex on the leeward wing, Figure 29 shows the evolution of the vortical structures from 16°to 2»°over the whole plane. Already at 16°but more marked at 20°and 2»°, one can see in the first quarter of the leeward wing the onset of one or several vortices (according to the incidence) parallel to the fuselage which are diverted towards the IBV vortex and interact with it according to various modalities. The interaction becomes more complex and occurs further downstream when the angle of attack increases. One can notice several local zones of flow reversal (in red) close to the location of this interaction. We believe that, downstream of this vortex-vortex interaction, the nature of the leeward IBV vortex changes from a jet-type (upwind of it) to a wake-type core (downwind of it). Interestingly, this interaction does not occur at 18°where k-𝜔 SST is able to predict an almost zero roll moment in agreement with TUM measurements while for the incidences for which this interaction occurs, the roll moment keeps its negative value when computed with DDES-SST or k-𝜔 SST. This local flow behaviour might be coincidental with no link of causality with the change of slope of the roll moment, but it would be valuable to pay more attention to this early vortex-vortex interaction, to assess if it is strongly dependent

VI. Conclusions

This article was devoted to a global and local study of the flow around two different wing configurations, the so-called NA1_W1 and NA1_W2 geometries, for several angles of attack and two drift conditions. Two statistic k-𝜔 SST and hybrid RANS/LES DDES-SST turbulence models were used on flow-adapted grids which capture automatically and efficiently the prominent vortical structures during their evolution. Systematic comparisons with TUM measurements were conducted on the forces and moments but also on the local flow details. Systematic vortex core analysis were provided, which can be used to assess the ability of CFD to represent the various jet-type or wake-type vortex behaviour all along the path of the main vortical structures present for each flow configuration. For the symmetric condition, the simulation performed with the help of DDES-SST, is globally equivalent to the one performed with k-𝜔 SST for low angles of attack although DDES-SST is the only model able to capture accurately the vortex breakdown. For high angles of attack, the hybrid RANS/LES clearly shows its superior reliability. We have also noticed that the configuration NA1_W2 appears easier to simulate than NA1_W1 where vortex interaction occurs at lower angles of attack. The side-slip condition appears more complex to simulate for both wing configurations. Particularly, the lateral roll moment change of sign occurring at a critical angle of attack is not satisfactorily modelled by the present computational approach, whatever the turbulence model used. These computations have at least established that this weakness is not related to a lack of local grid density but rather due to a modelling error affecting the vortical flow on the leeward wing. This is disappointing because the reliable prediction of the roll moment crisis is a key feature to assess the plane stability during maneuvers. By analyzing the local flow on the leeward wing for several incidences located close to the change of slope of C mx , we propose that this modelling error is due to the early interaction between the computed fuselage vortices and the IBV vortex, not far from its onset. Based on TUM experiments, we suppose that this interaction does not occur in the same way in the real flow. This interaction strongly affects the core flow velocity longitudinal evolution and therefore impacts the very sensitive roll moment coefficient. This hypothesis will be assessed more rigorously in future studies.

Fig. 1

 1 Fig. 1 Wing planforms

( a )

 a 𝑘 -𝜔 SST 𝛼 = 24°, 𝛽 = 0°(b) 𝑘 -𝜔 SST 𝛼 = 24°, 𝛽 = 5°( c) DDES-SST 𝛼 = 24°, 𝛽 = 0°(d) DDES-SST 𝛼 = 24°, 𝛽 = 5°F ig. 2 NA1_W1 -Views of the flow-adapted grids at the cross-section x/c r =0.592 with and without drift. Simulations performed with 𝑘 -𝜔 SST and DDES-SST.

Fig. 3 Fig. 4

 34 Fig. 3 Longitudinal aerodynamic coefficients for the NA1_W1 and NA1_W2 configurations versus the angle of attack 𝛼 and 𝛽 = 0°.

Fig. 5 Fig. 6 Fig. 7

 567 Fig. 5 Surface pressure distribution and skin friction lines on the NA1_W1 and NA1_W2 configurations for 𝛼 = 16°and 𝛽 = 0°.

Fig. 9

 9 Fig. 8 Vortical structures of the NA1_W1 and NA1_W2 configurations for 𝛼 = 24°and 𝛽 = 0°.

Fig. 11

 11 Fig. 11 Dimensionless axial velocity for IBV and MBV of the NA1_W1 and NA1_W2 configurations for 𝛼 = 24°a nd 𝛽 = 0°.

Fig. 12

 12 Fig. 12 Vortical structures of the NA1_W1 and NA1_W2 configurations for 𝛼 = 32°and 𝛽 = 0°.

  Fig. 13 Surface pressure distribution and skin friction lines on NA1_W1 and NA1_W2 configurations for 𝛼 = 32°and 𝛽 = 0°.

Fig. 14

 14 Fig. 14 Axial component of the velocity at x/c r = 0.592 for the NA1_W1 and NA1_W2 configurations at 𝛼 = 32°a nd 𝛽 = 0°.

Fig. 15

 15 Fig. 15 Dimensionless axial velocity for IBV and MBV of the NA1_W1 and NA1_W2 configurations for 𝛼 = 32°a nd 𝛽 = 0°.

  20°and 2»°, the roll moment coefficient decreases again instead of keeping close to zero, as indicated in the experiments. This change of trend between 16°, 18°and 20°needs to be correlated to the local behaviour of the leeward and windward vortices. Of course, the range of variation of the roll moment being very small, it is hard to simulate accurately and to track in the local flow characteristics.

Fig. 16

 16 Fig. 16 Longitudinal and lateral aerodynamic coefficients for the NA1_W1 configuration for 𝛽 = 5°versus the angle of attack 𝛼.

Fig. 17

 17 Fig. 17 Longitudinal and lateral aerodynamic coefficients for the NA1_W2 configuration for 𝛽 = 5°versus the angle of attack 𝛼.

Fig. 18

 18 Fig. 18 Vortical structures of the NA1_W1 and NA1_W2 configurations for 𝛼 = 16°and 𝛽 = 5°.

Fig. 20

 20 Fig. 20 Axial component of the velocity at x/c r = 0.592 for the NA1_W1 and NA1_W2 configurations at 𝛼 = 16°a nd 𝛽 = 5°.

Fig. 21 Fig. 22

 2122 Fig. 21 Leeward wing -Dimensionless axial velocity for the IBV and MBV of the NA1_W1 and NA1_W2 configurations for 𝛼 = 16°.

Fig. 23 Fig. 24

 2324 Fig. 23 Vortical structures of the NA1_W1 and NA1_W2 configurations for 𝛼 = 24°and 𝛽 = 5°.

Fig. 25

 25 Fig. 25 Axial component of the velocity at x/c r = 0.592 for the NA1_W1 and NA1_W2 configurations at 𝛼 = 24°a nd 𝛽 = 5°.

Fig. 26 Fig. 27

 2627 Fig. 26 Leeward wing -Dimensionless axial velocity for the IBV and MBV of the NA1_W1 and NA1_W2 configurations for 𝛼 = 24°and 𝛽 = 5°.

( a )

 a Fig. 28 NA1_W1 -Axial component of the velocity at x/c r = 0.592 between 𝛼 = 16°and 𝛼 = 24°for 𝛽 = 5°.

Acknowledgments

The computations were performed using HPC resources from GENCI (Grand Equipement National de Calcul Intensif) (Grant-A0082A00129, Grant-A0102A00129) which is gratefully acknowledged.