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

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. 1 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. 2 At a sufficiently high angle of attack, the vortices undergo a sudden expansion known as vortex breakdown, which was first observed by Werle 3 on a slender delta wing. By using the Particle Image Velocimetry (PIV), Ol 4 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 5 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. 6 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. 7 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 8 investigated qualitatively and quantitatively the flow characteristics of a non-slender diamond wing by 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. 9 The case comprises a triangular wing with a symmetric NACA64A006 profile, see  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. 10,11 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.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.490759 m.
In the present calculations, a positive inflow angle α corresponds to an upward directed flow.
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 ρ is the fluid density, S ref the reference area, c r the chord of the wing at the root, V ∞ the inflow velocity (V ∞ = (u, v, w)) and ω the specific turbulence dissipation rate. The quantity Q is used to identify vortex regions, see e.g. Kamkar. 12 Quantity Equation

III.C. ISIS-CFD
The solver ISIS-CFD, available as a part of the FINE T M /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 13, 14 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. 15 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. 16 The technique included for the 6 degrees of freedom simulation of ship motion is described by Leroyer & Visonneau. 17 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, 18 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 flowrelated criteria. 19 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.  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.

III.C.1. Grid generation
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. Overviews of the surface grid on the wing for both meshes are given in Figure 4.

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. 20 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. 21 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. 21 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) 22 and ((I)D)DES approaches can be used, with Partially Averaged Navier Stokes (PANS) and LES approaches currently being studied, see Pereira et al. 23 The Spalart correction (proposed by Dacles-Mariani et al. 24 ) 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·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. 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 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.

IV.A. Prediction of forces and moments
Several turbulence models, namely Spalart-Allmaras (SA), k − ω 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 − ω 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.  . 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 × 0.58mm. This leads to a grid density of about 60 × 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. . 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 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. 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 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.