# Incompressible flow calculations of blunt leading edge separation for a 53 degree swept diamond wing

Michel Visonneau, Emmanuel Guilmineau, Serge Toxopeus

► To cite this version:

Michel Visonneau, Emmanuel Guilmineau, Serge Toxopeus. Incompressible flow calculations of blunt leading edge separation for a 53 degree swept diamond wing. 53rd AIAA Aerospace Sciences Meeting, AIAA Science and Technology Forum, Jan 2015, Orlando, United States. 10.2514/6.2015-0065. hal-01202526

HAL Id: hal-01202526

https://hal.science/hal-01202526

Submitted on 8 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Incompressible flow calculations of blunt leading edge separation for a 53 degree swept diamond wing

Michel Visonneau*

ECN/CNRS, Nantes, Cedex 3, France

Emmanuel Guilmineau†

ECN/CNRS, Nantes, Cedex 3, France

Serge Toxopeus‡

Maritime Research Institute Netherlands (MARIN), Wageningen, 6708 PM, The Netherlands

One of the challenges in simulating aerodynamic flows is the accurate prediction of onset and progression of separation. In earlier STO research, it was found that various predictions for three-dimensional airfoils at angle of attack showed different pitch moments due to differences in the predicted location of separation areas. Therefore, a diamond shaped wing was developed within the NATO STO AVT-183 research group, to isolate the blunt leading edge separation and understand the underlying physical mechanisms. In the present paper, two viscous-flow solvers dedicated to hydrodynamic applications are used to predict the flow around the diamond wing. Both codes solve the incompressible flow using unstructured grids and finite volume discretisation. The results comprise different grid set-ups, with/without peniche and different turbulence models ranging from isotropic or anisotropic statistical turbulence closures to hybrid LES turbulence models. Moreover, the role played by the discretisation error is assessed by using anisotropic automatic grid refinement procedures based on flow-related criteria. A detailed discussion is made, highlighting the similarities and differences of the results, the respective influence of modeling and discretisation errors and showing the potential of CFD tools to predict the type of flow under consideration.

## I. Introduction

Delta wings are used on many aircrafts, and as these aircrafts become more and more maneuverable, it is necessary to understand very well the vortex dynamics and the underlying flow physics.¹ The flow past a delta wing is dominated by the vortical structures that originate from its leading edges. Several variables like the angle of attack, leading-edge geometry, sweep angle, Reynolds number, ..., affect the vortex dynamics.² At a sufficiently high angle of attack, the vortices undergo a sudden expansion known as vortex breakdown, which was first observed by Werle³ on a slender delta wing. By using the Particle Image Velocimetry (PIV), Ol⁴ investigated the flow structure of the delta wing for two sweep angles $Λ = 50°$ and $Λ = 65°$ for low Reynolds number. He found that the $50°$ wing exhibited similar flow characteristics shared by slender delta wings and wings of high aspect ratio while results for the $65°$ wing are in accordance with the literature. Yaniktepe and Rockwell⁵ investigated aerodynamics of the delta wing with a sweep angle of $Λ = 38.7°$ by using PIV and reported that the vortex breakdown beginning from the wing apex occurs much more explicitly at a high angle of attack. Gursul et al.⁶ investigated unsteady aerodynamics of non-slender delta wings covering a broad spectrum of topics including shear layer instabilities, structure of non-slender breakdown. They reported that the separated flow region occurs before the formation of helical vortices.

Canpolat et al.⁷ observed the variation of flow structures on the delta wing surface with a sweep angle of $Λ = 40°$ as a function of the angle of attack and the yaw angle. With a yaw angle, the symmetry flow structure disappears, and the vortex breakdown occurs earlier on the windward side of the delta wing. Sahin et al.⁸ investigated qualitatively and quantitatively the flow characteristics of a non-slender diamond wing by

<footnote>*CNRS Research Director - Dynamics of Marine Propulsion Systems LHEEA Lab (ex LMF), CNRS UMR6598</footnote>

<footnote>†CNRS Researcher - Dynamics of Marine Propulsion Systems LHEEA Lab (ex LMF), CNRS UMR6598</footnote>

<footnote>‡Senior Researcher, R&D - Team CFD</footnote>
varying the yaw angle. They found that the formation of the vortex breakdown and vortical flow structures were sensitive to the yaw angle.

The objective of the collaborative action performed under the umbrella of NATO STO/AVT-183 group, entitled "Reliable Prediction of Separated Flow Onset and Progression for Air and Sea Vehicles" is to have a more detailed look at the basic physical problems associated with the prediction of onset and progression of vortical flow around rounded leading edges on swept configurations. By comparing various numerical and modeling approaches, one should be able to select the most promising approaches and give recommendations for further improvements in terms of numerical algorithms or turbulence modeling.

## II. Test case

In order to be able to focus on onset and progression of flow separation, a unit test case has been defined, see Luckring and Boelens.⁹ The case comprises a triangular wing with a symmetric NACA64A006 profile, see Figure 1. The leading edge sweep is 53° and the trailing edge sweep 26.5°. The chord $c_r$ of the wing is 1.2$m$, the half-span $b/2 = 0.66m$. The moment reference point is specified at 0.490759$m$ aft of the nose of the wing. The plan-form area $S_{ref}$ of the half-wing is 0.394385$m^2$.

<figure><img src="image_2.png" /><figcaption>Figure 1. Geometry of Diamond Wing: NACA64A006 profile, 53° leading edge sweep and 26.5° trailing edge sweep (left: profile, right: planform)</figcaption></figure>

Model tests at $Re=4.04 \ 10^6$ based on $C_r$, have been performed by Technische Universität Münich. Stereo PIV and hot-wire measurements were done to acquire the mean and fluctuating flow fields and the forces and moments on the wing for various inflow angles have been obtained, see Hövelmann et al.¹⁰,¹¹ To avoid influence of a boundary layer on the flow around the half-wing during the wind-tunnel tests, a peniche or stand-off with a height of 0.09 $m$ was used. The forces on the peniche are not taken into account in the forces presented in this report.

To be able to import the CAD geometry (IGES surface) as a Parasolid into Hexpress (see below), the geometry had to be slightly altered: the tip of the wing was slightly truncated and a small cap was placed at the truncation. Furthermore, the trailing edge was modified by rounding the last 1% of the profile.

## III. Numerical methods

### III.A. Coordinate system

The origin of the right-handed system of axes used in this study is located at the symmetry plane of the profile, at the nose, at the intersection of the wing with the peniche. The longitudinal axis $x$ is positive aft and the vertical axis $z$ is positive upward. The forces and moments presented in this report are given according to this coordinate system, see Figure 2. The pitch moment is given with respect to $x = 0.490759m$.

In the present calculations, a positive inflow angle $\alpha$ corresponds to an upward directed flow.
<figure><img src="image_3.png" /><figcaption>Figure 2. Nomenclature and sign convention</figcaption></figure>

### III.B. Data reduction and non-dimensionalisation

The lift $L$ and drag $D$ on the wing are derived from the computed forces, see also Figure 2:

$$L = -X \cdot \sin \alpha + Z \cdot \cos \alpha \quad D = X \cdot \cos \alpha + Z \cdot \sin \alpha$$

In this report, non-dimensional values will be given, unless otherwise indicated. The formula in the table below are used to make the results non-dimensional, where $\rho$ is the fluid density, $S_{ref}$ the reference area, $c_r$ the chord of the wing at the root, $V_\infty$ the inflow velocity ($\bar{V}_\infty = (u, v, w)$) and $\omega$ the specific turbulence dissipation rate. The quantity $Q$ is used to identify vortex regions, see e.g. Kamkar.¹²

<table><thead><tr><th>Quantity</th><th>Equation</th></tr></thead><tbody><tr><td>Forces</td><td>X', Z', C<sub>L</sub>, C<sub>D</sub> = <sup>X, Z, L, D</sup>&frasl;<sub><sup>1</sup>&frasl;<sub>2</sub>ρV<sub>∞</sub><sup>2</sup>S<sub>ref</sub></sub></td></tr><tr><td>Moment</td><td>C<sub>M</sub> = M' = <sup>M</sup>&frasl;<sub><sup>1</sup>&frasl;<sub>2</sub>ρV<sub>∞</sub><sup>2</sup>S<sub>ref</sub>c<sub>r</sub></sub></td></tr><tr><td>Pressure</td><td>c<sub>p</sub> = <sup>p - p<sub>∞</sub></sup>&frasl;<sub><sup>1</sup>&frasl;<sub>2</sub>ρV<sub>∞</sub><sup>2</sup></sub></td></tr><tr><td>Vorticity</td><td>ω'<sub>i</sub> = ω<sub>i</sub> &middot; <sup>c<sub>r</sub></sup>&frasl;<sub>V<sub>∞</sub></sub></td></tr><tr><td>Shear stress</td><td>c<sub>fi</sub> = <sup>f<sub>i</sub></sup>&frasl;<sub><sup>1</sup>&frasl;<sub>2</sub>ρV<sub>∞</sub><sup>2</sup></sub></td></tr><tr><td>Turbulent kinetic energy</td><td>k' = <sup>k</sup>&frasl;<sub>V<sub>∞</sub><sup>2</sup></sub></td></tr><tr><td>Turbulence dissipation rate</td><td>ε' = k &middot; ω <sup>c<sub>r</sub></sup>&frasl;<sub>V<sub>∞</sub><sup>3</sup></sub></td></tr><tr><td>Q</td><td>Q = <sup>1</sup>&frasl;<sub>2</sub> (|Ω|^2 - |S|^2) &middot; c<sub>r</sub><sup>2</sup>/V<sub>∞</sub><sup>2</sup></td></tr><tr><td>Helicity</td><td>H' = <sup>H</sup>&frasl;<sub>|V| &middot; |ω|</sub></td></tr></tbody></table>

### III.C. ISIS-CFD

The solver ISIS-CFD, available as a part of the FINE<sup>TM</sup>/Marine computing suite, is an incompressible unsteady Reynolds-averaged Navier-Stokes (URANS) method mainly devoted to marine hydrodynamics. 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¹³,¹⁴ with or without rotation corrections. All models are available with wall-function or low-Reynolds near wall formulations. Hybrid LES turbulence models based
on Detached Eddy Simulation (DES) are also implemented and have been validated on automotive flows characterized by large separations.¹⁵ Additionally, several cavitation models are available in the code.

The solver is based on the finite volume method to build the spatial discretization of the transport equations. The unstructured discretization is face-based. While all unknown state variables are cell-centered, the systems of equations used in the implicit time stepping procedure are constructed face by face. Fluxes are computed in a loop over the faces and the contribution of each face is then added to the two cells next to the face. This technique poses no specific requirements on the topology of the cells. Therefore, the grids can be completely unstructured; cells with an arbitrary number of arbitrarily-shaped faces are accepted. Pressure-velocity coupling is enforced through a Rhie & Chow SIMPLE type method: at each time step, the velocity updates come from the momentum equations and the pressure is given by the mass conservation law, transformed into a pressure equation. In the case of turbulent flows, transport equations for the variables in the turbulence model are added to the discretization. Free-surface flow is simulated with a multi-phase flow approach: the water surface is captured with a conservation equation for the volume fraction of water, discretized with specific compressive discretization schemes.¹⁶ The technique included for the 6 degrees of freedom simulation of ship motion is described by Leroyer & Visonneau.¹⁷ Time-integration of Newton's law for the ship motion is combined with analytical weighted or elastic analogy grid deformation to adapt the fluid mesh to the moving ship. To enable relative motions of appendages, propellers or bodies without having recourse to overlapping grids, a sliding grid approach has been implemented. Propellers can be modeled using actuator disc theory, by coupling with boundary element codes (RANS BEM coupling,¹⁸ or with direct discretization through e.g. the rotating frame method or sliding interface approaches. Finally, an anisotropic automatic grid refinement procedure has been developed which is controlled by various flow-related criteria.¹⁹

Parallelization is based on domain decomposition. The grid is divided into different partitions, which contain the cells. The interface faces on the boundaries between the partitions are shared between the partitions; information on these faces is exchanged with the MPI (Message Passing Interface) protocol. This method works with the sliding grid approach and the different sub-domains can be distributed arbitrarily over the processors without any loss of generality. Moreover, the automatic grid refinement procedure is fully parallelized with a dynamic load balancing working transparently with or without sliding grids.

### III.C.1. Grid generation

For the wing alone, the size of the computational domain is [-2.4 m; 7.2 m]×[0 m; 6.0 m]×[-6.0 m; 6.0 m] in X, Y and Z directions. Y = 0 corresponds to the symmetry plane of the wing. The origin of the computational domain is located at the trailing edge of the wing in the symmetry plane. The computational domain is shown in Figure 3. The mesh is generated by the unstructured automatic mesh generator Hexpress. This software generates meshes containing only hexahedral cells. The peniche is not taken into account in this mesh.

A second mesh is created with the anisotropic automatic grid refinement procedure. The criterion used is the pressure Hessian tensor with a minimum size of $8×10^{-4}$ m. Concerning the boundary conditions, a uniform velocity $V_∞$ is prescribed at the inlet plane. The same boundaries conditions are applied on the sides $Y_{max}$, $Z_{min}$ and $Z_{max}$. For the symmetry plane, a symmetry condition is used. At the outlet, the pressure is imposed. A no-slip condition is used on the wing. In this case, the averaged value of $y^+$ is 0.29 with a maximum value 1.16. The total number of cells and the number of faces on boundaries for both meshes are given in Table 1.

<table><caption>Table 1. Number of cells and faces in the grid for ISIS-CFD</caption><thead><tr><th></th><th>Initial mesh</th><th>Refined mesh</th></tr></thead><tbody><tr><td>Total (cells)</td><td>8,590,874</td><td>17,571,913</td></tr><tr><td>Wing</td><td>233,082</td><td>255,039</td></tr><tr><td>Outlet</td><td>160</td><td>160</td></tr><tr><td>Inlet</td><td>302</td><td>302</td></tr><tr><td>Far field</td><td>516</td><td>516</td></tr></tbody></table>

Overviews of the surface grid on the wing for both meshes are given in Figure 4.
<figure><img src="image_4.png" /><figcaption>Figure 3. Computational domain and boundaries for the flow solver ISIS-CFD</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Overview of the surface grid for ISIS-CFD</figcaption></figure>

### III.D. ReFRESCO

REFRESCO is a viscous-flow CFD code that solves multiphase (unsteady) incompressible flows with the RANS equations, complemented with turbulence closure models, cavitation models and volume-fraction transport equations for different phases, see Vaz et al.²⁰ The equations are discretized using a finite-volume approach with cell-centered collocated variables and in strong-conservation form. A pressure-correction equation based on the SIMPLE algorithm is used to ensure mass conservation as discussed by Klaij and Vuik.²¹ Time integration is performed implicitly with first or second-order backward schemes. At each implicit time step, the non-linear system of velocity and pressure is linearized with Picard's method and either a segregated or coupled approach is used. In the latter, the coupled linear system is solved with a matrix-free Krylov subspace method using a SIMPLE-type preconditioner.²¹ A segregated approach is always adopted for the solution of all other transport equations. The implementation is face-based, permitting grids with elements consisting of an arbitrary number of faces (hexahedra, tetrahedra, prisms, pyramids, etc.),
and, if needed, h-refinement i.e. hanging nodes. State-of-the-art CFD features such as moving, sliding and deforming grids, as well automatic grid refinement are also available in the code.

For turbulence modeling, RANS/URANS, Scale Adaptive Simulation (SAS)²² and ((I)D)DES approaches can be used, with Partially Averaged Navier Stokes (PANS) and LES approaches currently being studied, see Pereira et al.²³ The Spalart correction (proposed by Dacles-Mariani et al.²⁴) to limit the production of turbulence kinetic energy based on the stream-wise vorticity can be activated.

The code is parallelised using MPI and sub-domain decomposition, and runs on Linux workstations and HPC clusters. REFRESCO is currently being developed, verified and validated at MARIN in the Netherlands in collaboration with IST (Lisbon, Portugal), USP-TPN (University of São Paulo, Brazil), Delft University of Technology, the University of Groningen, the University of Southampton, the University of Twente and recently Chalmers University of Technology.

Automatic wall functions are available. However, for all cases presented in this study the $y_2^+$ values are below the threshold value for the automatic wall functions of 5, such that the equations are integrated down to the wall.

### III.D.1. Grid generation

The computational domain is formed by a horizontal cylinder with a radius of $4 \cdot c_r$, with the origin located at the trailing edge of the wing. In transverse direction, the cylinder extends from $y = -0.09 \, m$ to $y = 5.91 \, m$. In Figure 5, the computational domain is shown.

<figure><img src="image_6.png" /><figcaption>Figure 5. Computational domain and boundaries</figcaption></figure>

Unstructured grids with hexahedral elements were generated using Hexpress, in several steps. First, the top half of the domain was modeled. Then, when the quality of the grid was deemed good enough, the grid was mirrored to obtain the full domain.

For both the half and full domain grids, an inflation layer was added to the surfaces of the wing, peniche and tunnel wall to capture the gradients in the boundary layer. The thickness of the layer was based on the anticipated Reynolds number corresponding to wind tunnel conditions. The maximum $y^+$ value was below 1.7 for all cases (on the tunnel wall), indicating that the wall function is not activated. The maximum and average $y^+$ values on the wing surface were found to be about 0.34 and 0.12 respectively.

Two grids were generated, one initial grid (designated *coarse*) and a second grid (designated *fine*) in which the refinement close to the wing was increased. The total number of cells and the number of faces
on each boundary is given in Table 2. It is seen that in order to capture the details of the flow around the wing and peniche, a large increase of the number of cells in their vicinity is required.

Overviews of the surface grid on the wing, peniche and tunnel wall are given in Figure 6.

<table><caption>**Table 2.** Number of cells and faces in the grids</caption><thead><tr><th></th><th>Coarse<br>(faces)</th><th>Fine<br>(faces)</th></tr></thead><tbody><tr><td>Total (cells)</td><td>21,523,724</td><td>26,519,958</td></tr><tr><td>Wing</td><td>356,116</td><td>430,286</td></tr><tr><td>Peniche</td><td>32,300</td><td>57,344</td></tr><tr><td>Tunnel</td><td>47,816</td><td>53,476</td></tr><tr><td>Far-field</td><td>916</td><td>916</td></tr><tr><td>External</td><td>6,844</td><td>6,696</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Overview of the surface grid (left: on the wing, peniche and tunnel wall, right: on the wing and peniche)</figcaption></figure>

At the external boundary of the domain, the AutoDetect boundary condition is used. With this boundary condition, inflow or outflow faces are automatically detected based on the local flow direction. Symmetry conditions are posed on the far-field boundary. On the surfaces of the wing, peniche and tunnel wall no-slip wall boundaries are applied.

# IV. Analysis of the results

## IV.A. Prediction of forces and moments

Several turbulence models, namely Spalart-Allmaras (SA), $k-\omega$ SST, EARSM and DES-SST are compared. For ISIS-CFD, the simulations were done on very similar fine grids resulting from the use of the automatic grid refinement procedure based on the Hessian of the pressure. Drag and lift predictions provided by the statistical closures are very similar and only in fair agreement with the measurements (see Table 3. From these models, the SA results tend to be closest to the measurements, with an error in drag of about 4% and in lift of 2%. One notices that EARSM provides a larger drag, which is consistent with the fact that the modeled turbulence anisotropy tends to increase the longitudinal vorticity production, and consequently, the drag. Pitch moment predicted by these turbulence closures are very inaccurate with errors ranging from 23% for SA, 30% for $k-\omega$ SST and 39% for EARSM. On the other hand, the hybrid LES closure DES-SST provides far better predictions of the drag with an error of 4.44% and an almost perfect evaluation of the lift (0.08%). The pitch moment is also very accurately computed with an error of 4.05% with TUM's experiments. The first comparisons illustrate a very strong dependency of global quantities (like forces and moments) on the turbulence closures, which convincingly suggests the sensitivity of the simulation to modeling errors.
<table><caption>Table 3. Comparison of forces and moment</caption><thead><tr><th rowspan="2"></th><th colspan="3">ISIS-CFD</th><th colspan="2">ReFRESCO</th><th rowspan="2">Experiments</th></tr><tr><th>k-ω SST</th><th>EARSM</th><th>DES</th><th>k-ω SST</th><th>Spalart-Allmaras</th></tr></thead><tbody><tr><td>Cd</td><td>0.0840</td><td>0.0849</td><td>0.0834</td><td>0.0917</td><td>0.0829</td><td>0.0798</td></tr><tr><td>Cl</td><td>0.4833</td><td>0.4603</td><td>0.5164</td><td>0.5394</td><td>0.5266</td><td>0.5167</td></tr><tr><td>Cm</td><td>-0.0356</td><td>-0.0382</td><td>-0.0285</td><td>-0.0227</td><td>-0.0212</td><td>-0.0274</td></tr></tbody></table>

## IV.B. Pressure surface distribution and wall streamlines

Figures 7 show global views of the skin friction lines over upper side of the wing for the various turbulent closures considered in this study for ISIS-CFD and ReFRESCO results while, for ISIS-CFD only, figures 8 show, in addition the surface pressure distribution and the trajectory of a particle injected in the vicinity of the onset of the main vortex. First of all, an island of low pressure is visible from the leading edge to the inner part of the wing in correlation with the development of the main vortex. Its extent depends strongly on the turbulence closures, larger with DES-SST and smaller with $k-ω$ SST and EARSM. Although the wall streamlines present similarities at first glance, a closer look reveals some significant differences. The point of onset of the main open separation is located closer to the apex with EARSM than with $k-ω$ SST or DES-SST. All turbulence closures predict the existence of a weak secondary vortex associated with the convergent wall streamlines which are almost parallel to the root of the wing. Most important are the comparisons of the wall streamlines associated with the main vortical structures shed close to the leading edge. For $k-ω$ SST, we can distinguish, progressing from the leading edge to the interior of the upper side of the wing, a succession of three divergent, convergent and divergent streamlines. The same succession of divergent, convergent and divergent wall streamlines is visible in the EARSM results but these lines of accumulation are more marked and maintained almost up to the trailing edge where a destabilization of the vortex is indicated by the wall streamlines complex topology. DES-SST provides a somewhat different behavior since there are only two convergent and divergent streamlines which are visible. The convergent streamlines associated with the open separation are more visible and maintained without any sign of vortex destabilization up to the trailing edge. The trajectories of a particle injected close to the expected point of onset complete the comparisons and illustrate the vortex breakdown which occurs in $k-ω$ SST and EARSM computations close to the trailing edge, while a far more stable vortical motion is predicted by DES-SST. Figures 15 provide detailed views of the wall streamlines close to the expected onset region to clarify the previous analysis. We will proceed at the end of this study to a detailed analysis of the flow close the point of open separation. For ReFRESCO, the most significant difference between SST and SA is the earlier separation in SST and the stronger secondary inboard vortex which is hardly seen in the SA streamlines. Furthermore, SST shows strong flow separation and complex flow patterns at the wing tip, while the SA streamlines are more regular (showing reversed flow as well, however).

## IV.C. Cross section comparisons

TUM provided detailed PIV measurements of the velocity (three components), longitudinal vorticity and turbulent kinetic energy for cross-sections ranging from $x/c_r=0.100, 0.200, 0.295, 0.305, 0.395, 0.405, 0.500$ and $0.695$. Several roughness thicknesses were tested and we have decided, like all the participants of the AVT183 collaborative action, to retain the experiments performed with a thickness of $150μm$ since they appeared in best agreement with computations where the transition effects are not modeled. Thanks to the use of the automatic grid refinement procedure implemented in ISIS-CFD, a large number of grid points have been automatically included into the core region during the spatial development of the vortex, with a minimum cell size of $0.58mm \times 0.58mm$. This leads to a grid density of about $60 \times 30$ in the horizontal and vertical directions of a vortex cross-section, a density which is kept more or less unchanged all along the vortex progression (see, for instance, figure 9 showing the grid density at $x/c_r = 0.395$). Such a grid density appears sufficient to eliminate any significant influence of the numerical diffusion on the comparisons. We can therefore consider that the following comparisons can be used with confidence to assess the actual predictive capabilities of the considered turbulence models.
<figure><img src="image_8.png" /><figcaption>Figure 7. Wall streamlines</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Surface pressure distribution, skin friction lines and particle trajectories</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Mesh distribution in cross section $x/c_r = 0.395$</figcaption></figure>

IV.C.1. Wall pressure coefficient

Figures 10 show the wall pressure coefficient for cross-sections ranging from $x/c_r=0.100, 0.200, 0.295, 0.305, 0.395, 0.405, 0.500$ and $0.600$. At $x/c_r=0.100$, the peak of pressure appears to be better captured by the statistical closures. From $x/c_r=0.295$ to $x/c_r=0.305$, the Spalart-Allmaras model still predicts a unique peak of pressure occurring close to the leading edge while all the other turbulence models start to predict
a plateau followed by a peak located away from the leading edge. However, only the hybrid LES DES-SST is able to predict the right level of the pressure coefficient. From $x/c_r=0.395$ to $x/c_r=0.600$, the peak of pressure is gradually reduced but all the statistical closures strongly under-estimate this maximum while DES-SST appears to be in very good agreement with the measurements.

<figure><img src="image_11.png" /><figcaption>(a) $x/c_r=0.100$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) $x/c_r=0.200$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) $x/c_r=0.295$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) $x/c_r=0.305$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e) $x/c_r=0.395$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f) $x/c_r=0.405$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(g) $x/c_r=0.500$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(h) $x/c_r=0.600$</figcaption></figure>

<caption>Figure 10. Distribution of pressure coefficient in the cross sections</caption>

IV.C.2. Longitudinal velocity component and secondary streamlines

Figure 11 provides the longitudinal velocity component distribution and secondary streamlines for cross sections $x/c_r=0.100$, 0.200, 0.295, 0.395, 0.5 and 0.6 for the afore-mentioned turbulence closures, namely, Spalart-Allmaras (SA), $k-\omega$ SST, EARSM and DES-SST. For the statistical turbulence closures (i.e. SA,
$k-\omega$ SST and EARSM), we can see the development of a main vortex characterized by a low-streamwise component $U$ in the core of the vortex. A secondary co-rotating vortical structure is visible although less pronounced. Both results are very similar and differ only by the intensity of the vortical motion inside the core of the vortices and the location of their onset, as previously observed in the wall streamlines. The ReFRESCO SST results show a much more complex flow near $x/C_r = 0.6$ which does not agree well with the experiments. EARSM, thanks to the addition of non-linear terms and the resulting modeled anisotropy, predicts more intense vortical structures which are created earlier on the wing profile. On the contrary, DES-SST computations show a significantly different landscape. The streamwise velocity component distribution in the core of the vortex shows higher values from $U=0.2$ to $0.6$ in a smaller core while the statistical closures showed a more extended core with smaller values of $U$ ranging from to $0.0$ to $0.3$. Note that the secondary streamlines are not centered with the minimum of the streamwise velocity component because the Cartesian components of the velocity are used to draw these pictures instead of the velocity components projected on a local referential aligned with the local axis of the vortex path and a plane normal to this axis. The experiments performed at TUM are shown in figure 11(f). They are in better agreement with the DES-SST results both in terms of the streamwise component magnitude and for the extent of the region with low streamwise velocity component.

<figure><img src="image_19.png" /><figcaption>Figure 11. Streamwise component of the velocity and secondary flow</figcaption></figure>
IV.C.3. Streamwise vorticity

Figure 12 show the streamwise vorticity component for the same cross sections for the statistical turbulence closures. With EARSM and ReFRESCO SST, the longitudinal vorticity is created earlier than with the SA and ISIS-CFD k-ω SST model, which is associated with the earlier onset previously noticed. After $x/c_r=0.500$, the core of the vortex appears to be less coherent, which is an indirect indication of an imminent vortex breakdown. This destabilization is predicted earlier with the non-linear anisotropic EARSM model than with the linear isotropic k-ω SST closure, a fact which is obviously related with the lower level of turbulent viscosity predicted by EARSM in the core of the vortex during its progression along the wing. The ReFRESCO SST results also show significant breakdown, while the SA results show a coherent structure of relatively low x-vorticity. Once again, DES-SST results are characterized by a more consistent core of longitudinal vorticity all along the vortex progression. The first sign of destabilization appears only from $x/c_r=0.600$ where the unique island of longitudinal of vorticity starts to separate into two distinct parts. It is also interesting to notice the island of counter-rotating longitudinal vorticity which is visible close to the leading edge in the DES-SST results. The experiments are in better agreement with the DES-SST results since they show a persistent core of intense vorticity up to $x/c_r = 0.500$. At $x/c_r = 0.600$, some diffusion of the vorticity starts to occur with a less intense vortical core, compared to its periphery. However, no indication of breakdown can be noticed, contrary to what is predicted with all the statistical isotropic or anisotropic closures. From that point of view, the experiments are in better agreement with DES-SST prediction, although this hybrid LES model simulates a slightly too intense vortical flow after $x/c_r = 0.500$.

IV.C.4. Turbulence kinetic energy

It would be tempting to conclude from these comparisons that the hybrid LES DES-SST closure should be obligatorily used when we want to simulate the onset and progression of turbulent longitudinal vortices over such airfoils. However, one should remember that the diffusion of a turbulent vortex depends on the turbulence intensity in its core from the point of onset to its advection along the wing. From that point of view, assessing the longitudinal vorticity intensity and the turbulence kinetic energy (TKE) in the core of a vortex is like looking at the two faces of the same coin. Consequently, to draw more solid conclusions, figure 13 provides the turbulence kinetic energy in the same cross sections for the same turbulence closures previously compared. The turbulence kinetic energy distributions predicted by the statistical closures ISIS-CFD k-ω SST and EARSM are globally similar in terms of intensity while DES-SST predicts vortices with significantly higher values of the turbulence kinetic energy in the core of the main vortex. The ReFRESCO results shows levels of TKE somewhat lower than those found in ISIS-CFD, but still reasonably close to the experiments. The ReFRESCO SST results, however, are much too low, probably due to the Spalart correction on turbulence production based on the longitudinal vorticity. In most of these sections, the ratio between the simulated and modeled TKE is around 6, meaning that these high levels of TKE can be mainly attributed to the high temporal variations of the velocity field at these locations. At first glance, we see that EARSM and k-ω SST are not able to capture the right level of TKE at $x/c_r=0.395$ while DES-SST predicts a significantly higher value at this cross-section. However, DES-SST fails to predict the gradual decay of TKE which is observed in the experiments from $x/c_r=0.5$. If we examine more closely the results for the cross-sections located closer to the onset (figure 14), we notice that DES-SST is the only model able to predict the correct production of TKE from $x/c_r=0.295$ to 0.450. The peak of TKE which is visible in the experiments for these cross-sections is also captured by DES-SST, although slightly over-estimated, while neither k-ω SST nor EARSM are able to detect it. These plots also clearly show that the location of the onset of leading edge flow separation is predicted early by EARSM and ReFRESCO SST and too far aft by SA.

IV.D. Detailed analysis of the open separation

Figures 15 provide detailed views of the wall streamlines and surface pressure distribution to identify the flow characteristics in the immediate vicinity of the 3D open separation. Globally, except ReFRESCO SA model which provides a topology of the skin friction lines which is very different with no visible line of convergence, all the other statistical turbulence models predict wall flow behaviors with many similarities. For ISIS-CFD and ReFRESCO k-ω SST models, we notice a succession of three divergent, convergent and divergent lines. The 3D separation takes place probably at the location of the minimum pressure, which
<figure><img src="image_20.png" /><figcaption>(a) ISIS-CFD - k-ω SST</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) ISIS-CFD - EARSM</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) ISIS-CFD - DES</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(d) ReFRESCO - k-ω SST</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(e) ReFRESCO - Spalart-Allmaras</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(f) Experiments</figcaption></figure>

<caption>Figure 12. Vorticity in X-direction</caption>

is quite pronounced in both computations. ISIS-CFD EARSM results are also very similar although the onset point occurs earlier as indicated previously. DES-SST shows a somewhat different skin friction lines topology with only two pronounced lines of convergence and divergence. The line of convergence located close to the leading edge is a bit sinuous, indicating probably a lack of temporal convergence of the results close to the wall. Similarly, the wall pressure distribution show three different minima, the first one being the most pronounced, which is also an indirect indication of a slight lack of convergence. However, despite these minor irregularities, the skin friction lines topology is typical of a 3D open separation with a line of separation located probably in the vicinity of the first minimum of pressure around $x/c_r = 0.295$.

V. Conclusions

This paper has provided a detailed study of the influence of turbulence closures on the 3D separation occurring on 3D wings with rounded leading edge. This study was carried out with two different incompressible flow solvers to give more solid foundations to the conclusions. Comparisons with low Mach detailed experiments conducted at TUM were used to assess the computations. The discretization error was maintained under control thanks to the use of grids which are locally very fine in the vortex core region. It appears that the statistical turbulence closures are not able to predict the right evolution of the main vortex. The location of the 3D separation line is highly dependent on the turbulence closure, although the topology of
<figure><img src="image_26.png" /><figcaption>Figure 13. Turbulence kinetic energy</figcaption></figure>

the wall streamlines in the vicinity of the onset appears less dependent on the turbulence closures, as soon as a fine enough mesh is used. However, the development of the 3D vortex, as described in the measurements performed in several cross-sections, is not correctly predicted by the statistical turbulence closures. Moving from linear isotropic turbulence closures to anisotropic non-linear does not bring any significant improvement. The hybrid LES turbulence model (here DES-SST) is the only closure able to predict both global and local flow characteristics with a remarkable consistent accuracy.

## Acknowledgments

The experimental data were provided in the framework of STO AVT task group 183 by the Institute of Aerodynamics and Fluid Mechanics of Technische Universität München, Germany. The study performed by the French authors was granted access to the HPC resources of IDRIS and CINES under the allocation 2a0129 made by GENCI.

## References

<footnote n="1">¹Heron, I. and Myose, R. Y., “On the impingement of a Von Kàrmàn vortex street on a delta wing,” *22nd Applied Aerodynamics Conference and Exhibit*, Providence, Rhode Island, 16-19 August 2004, AIAA 2004-4731.</footnote>

<footnote n="2">²Yaniktepe, B., *Origin and control of vortex breakdown of unmanned combat air vehicles*, Ph.D. thesis, Çukurova Uni-</footnote>
<figure><img src="image_27.png" /><figcaption>Figure 14. Turbulence kinetic energy close to the onset</figcaption></figure>

versity, Turkey, 2006.

³Werlé, H., "Quelques résultats expérimentaux sur les ailes en flèche, aux faibles vitesses, obtenus en tunnel hydrodynamique," *La Recherche Aéronautique*, Vol. 41, 1954, pp. 15–21.

⁴Ol, M. V., *The passage toward stall of nonslender delta wings at low Reynolds number*, Ph.D. thesis, California Institute of Technology, 2001.

⁵Yaniktepe, B. and Rockwell, D., "Flow structure on a delta wing of low sweep angle," *AIAA Journal*, Vol. 42, No. 3, 2004, pp. 513–523.

⁶Gursul, I., Gordnier, R., and Visbal, M., "Unsteady aerodynamics of nonslender delta wings," *Progress in Aerospace Sciences*, Vol. 41, 2005, pp. 515–557.

⁷Canpotat, C., Yayla, S., Sahin, B., and Akilli, H., "Dye visualization of the flow structure over a yawed nonslender delta wing," *Journal of Aircraft*, Vol. 46, No. 5, 2009, pp. 1818–1822.

⁸Sahin, B., Yayla, S., Canpolat, C., and Akilli, H., "Flow structure over the yawed nonslender diamond wing," *Aeropsace Science and Technology*, Vol. 23, 2012, pp. 108–119.

⁹Luckring, J. M. and Boelens, O. J., "A Unit-Problem Investigation of Blunt Leading-Edge Separation Motivated by AVT-161 SACCON Research," *NATO RTO AVT-189 Specialists Meeting on Assessment of Stability and Control Prediction Methods for NATO Air and Sea Vehicles*, No. RTO-MP-AVT-189-27, Portsdown West, UK, Oct. 2011.

¹⁰Hövelmann, A., Grawunder, M., Buzica, A., and Breitsamter, C., "Experimental Analyses on the Flow Field Characteristics of the AVT-183 Diamond Wing Configuration," *53rd AIAA Aerospace Sciences Meeting (SciTech)*, American Institute of Aeronautics and Astronautics, Kissimmee, Florida, Jan. 2015.

¹¹Hövelmann, A., Knoth, F., and Breitsamter, C., "Leading-Edge Roughness Effects on the Flow Separation Onset of the AVT-183 Diamond Wing Configuration," *53rd AIAA Aerospace Sciences Meeting (SciTech)*, American Institute of Aeronautics and Astronautics, Kissimmee, Florida, Jan. 2015.
<figure><img src="image_28.png" /><figcaption>Figure 15. Wall pressure distribution and skin friction lines</figcaption></figure>

<footnote n="12">¹²Kamkar, S. J., Jameson, A., Wissink, A. M., and Sankaran, V., “Feature-Driven Cartesian Adaptive Mesh Refinement in the Helios Code,” *48*th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, No. AIAA 2010-171, American Institute of Aeronautics and Astronautics, Orlando, Florida, Jan. 2010.</footnote>

<footnote n="13">¹³Deng, G. D. and Visonneau, M., “Comparison of explicit algebraic stress models and second-order turbulence closures for steady flow around ships,” *7*th Symposium on Numerical Ship Hydrodynamics, Nantes, France, July 1999, pp. 4.4–1–4.4–15.</footnote>

<footnote n="14">¹⁴Duvigneau, R., Visonneau, M., and Deng, G. B., “On the role played by turbulence closures in hull ship optimization at model and full scale,” *Journal of Marine Science and Technology*, Vol. 8, 2003, pp. 11–25.</footnote>

<footnote n="15">¹⁵Guilmineau, E., Deng, G. B., and Wackers, J., “Numerical simulation with a DES approach for automotive flows,” *Journal of Fluids and Structures*, Vol. 27, 2011, pp. 807–816.</footnote>

<footnote n="16">¹⁶Queutey, P. and Visonneau, M., “An Interface Capturing Method for Free-Surface Hydrodynamic Flows,” *Computers & Fluids*, Vol. 36, 2007, pp. 1481–1510.</footnote>

<footnote n="17">¹⁷Leroyer, A. and Visonneau, M., “Numerical methods for RANSE simulations of a self-propelled fish-like body,” *Journal of Fluids and Structures*, Vol. 20, 2005, pp. 975–991.</footnote>

<footnote n="18">¹⁸Deng, G. B., Queutey, P., Visonneau, M., and Salvatore, F., “Ship propulsion prediction with a coupled RANSE-BEM approach,” *V International Conference on Computational Methods in Marine Engineering, MARINE2013*, Hamburg, Germany, May 2013.</footnote>

<footnote n="19">¹⁹Wackers, J., Deng, G. B., Guilmineau, E., Leroyer, A., Queutey, P., and Visonneau, M., “Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation,” *Computers & Fluids*, Vol. 92, 2014, pp. 209–222.</footnote>

<footnote n="20">²⁰Vaz, G., Jaouen, F. A. P., and Hoekstra, M., “Free-Surface Viscous Flow Computations. Validation of URANS Code FRESCO,” *28*th International Conference on Ocean, Offshore and Arctic Engineering (OMAE), No. OMAE2009-79398, Honolulu, Hawaii, May 31–June 5 2009.</footnote>

<footnote n="21">²¹Klaij, C. M. and Vuik, C., “SIMPLE-type preconditioners for cell-centered, colocated finite volume discretization of incompressible Reynolds-averaged Navier-Stokes equations,” *International Journal for Numerical Methods in Fluids*, Vol. 71, No. 7, March 2013, pp. 830–849.</footnote>

<footnote n="22">²²Menter, F. and Egorov, Y., “A Scale-Adaptive Simulation Model using Two-Equation Models,” No. AIAA 2005-1095, American Institute of Aeronautics and Astronautics, 2005.</footnote>

<footnote n="23">²³Pereira, F. S., Vaz, G., and Eça, L., “On the use of Hybrid Turbulence Models,” *17*th Numerical Towing Tank Symposium (NuTTS), Marstrand, Sweden, Sept. 2014, pp. 135–140.</footnote>

<footnote n="24">²⁴Dacles-Mariani, J., Zilliac, G. G., Chow, J. S., and Bradshaw, P., “Numerical/experimental study of a wingtip vortex in the near field,” *AIAA Journal*, Vol. 33, Sept. 1995, pp. 1561–1568.</footnote>