A HYBRID RANS-LES METHOD WITH COMPRESSIBLE k-omegaSSTSAS TURBULENCE MODEL FOR HIGH REYNOLDS NUMBER FLOW APPLICATIONS

Original scientific paper Three-dimensional, compressible, viscous and transient transonic turbulent flow over the wing was simulated by a hybrid RANS-LES modelling method, combined with the compressible k-omegaSSTSAS turbulence model. This approach is based on dividing the contribution of the fluctuating and the averaged velocity fields in the subgrid tensor, and modelling each of them with its corresponding turbulent viscosity. The "RANS mode" is used in flow field domains which can be treated with acceptable accuracy as relatively steady, such as in the boundary layer, while the "LES mode" is applied in the dominantly unsteady regions, far from the wing. Discretization of the governing equations is performed by Finite Volume Method on unstructured mesh. The method has been verified on the Onera M6 wing. The parallelization is achieved by decomposing the mesh into sub-domains and using the Open MPI technology. The implementation of turbulence model has been done using OpenFOAM. The flow simulation was also performed using ANSYS Fluent, and the results of the two methods were compared mutually, and with the Onera M6 experiment.


Introduction
Accurate analysis and simulation is essential for the design challenges in aerospace industry.Over the past several decades there has been intense research in the area of Computational Fluid Dynamics (CFD).The CFD has become a relevant supplement to the wind tunnel and flight tests.Majority of the research in this area comes out from the aerospace industry, with eternal requirement of the computational cost decrease.As computational power continuously increases, the requirements for the more realistic physical simulations also increase.According to this fact, the attention has moved from potential flow problems to problems governed by the compressible Euler and Navier-Stokes equations [1].
In this paper the OpenFOAM and built-in finite volume scheme are used to discretize the equations over the calculation domain, and a hybrid RANS-LES k-omegaSSTSAS (Shear-Stress Transport Scale-Adaptive Simulation) turbulence model, based on the k-omegaSSTSAS model, initially formulated by Menter [2,3,4], to close the system of equations.As is clear from the explanations of the RANS and LES strategies, the LES approach would be favourable in terms of simulating the unsteady flow of pure turbulence [5].At high Reynolds numbers and complex geometries, this approach leads to an unfavourably high computational effort.The RANS approach has lower demands on the mesh resolution, especially in the boundary layer region, but it is limited in its ability to resolve unsteady domains of the flow field.The general idea of here applied hybrid method is to provide a turbulence model which is capable of allowing LES-predictions to a certain extent with a lower computational effort, in comparison with full LES simulations.This is achieved by leaving a part of the flow to the RANS model.It is clear from theory of the RANS simulation, that the application of a RANS model in the entire domain is not favourable in order to resolve the unsteady flow in the free stream.For the LES applications a relatively high mesh resolution is required, which can be pointed to the fact that the mesh for LES simulations should ideally contain cubic cells only.According to this, similar mesh resolution in all dimensions of the mesh would be desirable, i.e. in direction normal to the wall, wall parallel and span-wise direction.While the "LESmode" of the hybrid model resolves the unsteady flow in the free stream, the relaxation of the mesh resolution requirements lowers the computational effort, in contrast to the full LES simulations.This relaxation falls back on the application of the "RANS-mode" in the boundary layer, which does not require cubic cells in the corresponding region [5].
Finite Volume Method (FVM) is an ideal method for computing discontinuous solutions arising in compressible flows.Any discontinuity must satisfy the Rankine-Hugoniot jump condition [6,7], which is a consequence of conservation.Since FVM are conservative, they automatically satisfy the jump conditions and hence give physically correct weak solutions.The significance of this property becomes apparent where computational simulations frequently involve complex geometrical domains.An additional feature is the local conservativity of the numerical fluxes, i.e. the numerical flux is conserved from one discretization cell to its neighbour.This last feature makes the FVM quite attractive when modelling problems for which the flux is of importance, such as in fluid mechanics, heat and mass transfer, etc. [8].
Mesh generation deserves a special study [9].In this paper the attention is devoted to the mesh quality and resolution.The mesh resolution is of special importance in the zone of shock wave, because the velocity and other properties of the fluid flow change almost instantaneously across the shock.Therefore, mesh should be dense enough to ensure the accurate generation of the shock line along the boundary.Since it is a transient case, shock wave/turbulent boundary layer interaction (SWBLI) is also an unavoidable physical phenomenon of transonic flow, and it significantly influences the load distribution over the wing.In this paper one of the primary aims is analysis of the solutions obtained using unstructured static mesh, while some future applications of here proposed k-omegaSSTSAS turbulence model will be focused on the employment of the dynamic meshes around the boundary.The ability of the turbulence model to catch normal shock wave has also been investigated.All analyses have been done for one of the test cases of the Onera M6 wing, which is well documented in the literature, and for which experimental and numerical data are available [10].
The Onera M6 wing is a classic CFD validation case for external flows, because of its simple geometry combined with complexities of transonic flow (i.e.local supersonic flow, shocks, and turbulent boundary layer separation).It has practically become a standard for CFD codes testing, and it has been used as a validation case in numerous CFD papers over the years [11].For the comparisons, the pressure coefficients at sections along the span of the wing, obtained by the experiment and published in the AGARD report [10], are used.
OpenFOAM employs domain decomposition [12] to split the mesh and fields into a number of sub-domains and allocate them to separate processors.Applications then run parallel [13] on separate sub-domains, with communication between processors established by software that uses the MPI communications protocol.High performance computing in ANSYS Fluent is also achieved by choosing the number of processors to be employed in parallel computation.

Fluid flow analysis 2.1 Governing equations of fluid flow
As mentioned before, a hybrid RANS-LES method employs switching from the RANS to the LES equations, by using the SAS term in the ω equation.Reynolds-Averaged Navier-Stokes equations are derived by applying time averaging on Navier-Stokes equations, so getting the additional terms, called the "Reynoldsstresses".Determination of these terms can be done with the help of the Boussinesq assumption [14].Boussinesq assumption defines eddy viscosity to link the Reynoldsstresses with the mean rate of deformation of the fluid.This assumption yields μ t = υ t /ρ, where μ t represents the turbulent eddy viscosity and υ t the turbulent kinematic eddy viscosity.The RANS models are, however, assumed to be limited in their capability of resolving unsteady flow structures.This deficiency can mainly be addressed to the determination of the turbulence length scale, which is used to compute the Reynolds-stresses.The length scale depends on the diffusion terms, which yields a correlation of the length scale and the shear-layer thickness.Standard two-equation models always return the shear-layer thickness as an appropriate length scale, independent of any resolved content.As a consequence, these models are tending to dampen the unsteady structures of the flow, which corresponds to an overly amount of turbulent eddy viscosity.To resolve this, by applying the spatial filtering function to Navier-Stokes equations we get Large Eddy Simulation (LES) equations [5,15].Analogous to the time filtering in the RANS approach, the spatial filtering introduces new terms in the LES equations.These terms represent the so-called Sub-Grid-Scale (SGS) stresses, which result from the filtering of the Navier-Stokes equations in space.The SGS stresses have to be modelled based on the same argument as for the Reynolds-stresses.This is done by the Sub-Grid-Scale (SGS) model.This finally allows resolving the turbulence in an unsteady manner, but requires very high mesh resolution in all directions, and so produces the high computational effort requirements.More detailed explanations of the RANS and the LES equations can be found, for example, in [13].
In the hybrid RANS-LES method, the unsteady Navier-Stokes equations are solved on a single mesh with the eddy viscosity computed from RANS model near the wall, and from LES model away from the wall.In the wall-stress modelling approach, the LES equations are formally defined everywhere in the domain, with RANS equations solved on mesh near the wall.The coupling between RANS and LES acts in a fashion similar to the widely used wall-functions in RANS, i.e. the RANS solver takes information from the computed LES flow field, and returns back the results in the form of wall fluxes, i.e. like the shear stress and heat transfer at the wall [16].

Subgrid scale modelling: k-omegaSSTSAS turbulence model
The solution of the filtered Navier-Stokes system of equations enables only the large eddies to be resolved, leaving the small eddies still unresolved.For compressible flows, particularly for transonic flows in which turbulent heat flux, turbulent diffusion, and viscous diffusion may become significant, the SGS modelling process is far from satisfactory [13].
A major boundary layer effect is the separation from a surface under adverse pressure gradient conditions.Separation has a strong effect on the near-wall turbulence [17].The idea behind the SST model is to combine the best elements of the k−ε and the k−ω model with the help of a blending function F 1 .The F 1 is equal to one near the surface, and zero in the outer part and for free shear flows.It activates the Wilcox model in the near-wall region and the k−ε model for the rest of the flow [17].The idea behind the SST-SAS k−ω model is to add an additional production term -the SAS term -in the ω equation, which is sensitive to resolved (i.e.unsteady) fluctuations.When the flow equations resolve turbulence, the length scale based on velocity gradients is much smaller than that based on time-averaged velocity gradients.Hence to the von Kárman length scale, L vk , is an appropriate quantity to be used as a sensor for detecting unsteadiness.In regions where the flow is on the limit of going unsteady, the objective of SAS term is to increase ω.The result is that k and μ t are reduced so that the dissipating (damping) effect of the turbulent viscosity on the resolved fluctuations is reduced, thereby promoting the momentum equations to switch from steady to unsteady mode [4].The SSTSAS model has also demonstrated the capability of accurate separation predictions in numerous incompressible cases that are described in [4], as well in the test cases reported in [5].
The SAS model allows the original flow instability to develop into a turbulent spectrum down to the resolution limit of the mesh.The behaviour of the SAS model is therefore similar to that of a DES model: the attached boundary layers are solved like in a RANS model and the detached unsteady-state flow behind a body results in a LES-like solution.Regarding to the "state of the art" DES models, the advantage of the SAS model compared to DES is that the mesh spacing does not explicitly affect the RANS model.Similar to the DES formulation, the SAS model also benefits from a switch in the numerical treatment between the steady and unsteady regions [5].
The compressible k-omegaSSTSAS turbulence model is evaluated on the Onera M6 test case and the results will be shown and discussed later.
The formulation of the k-omegaSSTSAS model is as follows, where detailed explanation of parameters used in equations is given in [3,4,5]: Turbulent frequency equation: where " u denotes the second velocity gradient, and is evaluated as the magnitude of velocity Laplacian .
The coefficients, ϕ, of the model are functions of The CD kω denotes the cross diffusion term: The subgrid-scale viscous tensor model can be written as: . 3 An additional feature of the SST model is the introduction of an upper limit for the turbulent shear stress in boundary layers in order to avoid excessive shear-stress levels typically predicted with Boussinesq eddy-viscosity models [2].The subgrid eddy viscosity is defined as: where α SGS denotes the subgrid turbulent heat flux and can be calculated as α SGS = μ SGS /Pr t .
The absolute value of the strain rate, , is used in definition of the eddy viscosity instead of vorticity in order to make a model suitable for aerodynamic applications [2].
There is a subtle but important aspect concerning the c s limiter.When used inside the SST-SAS term P SAS , it is without problems even for very coarse meshes.In this limit, L vk will be increased, which in turn will reduce the impact of P SAS Eq. ( 2) thereby keeping the model in RANS mode.In other words, the mesh limiter would not affect the RANS limit of the model.Production limiter k P  of the SST-SAS model is only used in the production term of k-equation.The purpose of this limiter is its activation in the adverse pressure gradient flows [2].It can be activated in regions where realizability is not at stake.In the proposed modified compressible k-omegaSSTSAS turbulence model in OpenFOAM, when k−ε model is activated, the ε is evaluated as follows The temperature equation can be written as (under assumption of the absence of source terms and under constant property assumptions): The Prandtl number, Pr , is the fluid property, and the turbulent Prandtl number, t Pr , is set to a constant value, assuming an analogy between turbulent heat and mass transfer; turbulent Prandtl number is set to 1,0 by experiment and analysis for compressible fluid.
The constants for this model are:

The fluid solver
This section describes how the numerical schemes for terms are specified, such as derivatives in equations [18], to achieve the optimum calculation process.Numerical treatment of the terms in equations mostly depends on the mesh and physics which is employed.As it is a demanding case, great attention is dedicated to the specification of numerical schemes, in order to get stable and accurate convergence of the solution.In this section the numerical schemes applied to the specific terms in the equations will be discussed.This section is mostly dedicated to OpenFOAM, as it offers a great freedom of choosing and combining the numerical schemes.
In this work a compressible density-based flow solver is used: rhoPisoFoam; it is a transient solver for subsonic/transonic, laminar or turbulent flow [19].Density-based solver is also chosen in Fluent.This type of simulation, i.e. the transonic turbulent flow solver which uses a hybrid RANS-LES method, requires very careful choice of the numerical schemes for mathematical terms such as gradient ∇ , Laplacian Δ, and especially the divergence ∇  .Specific numerical schemes make crucial influence to the accuracy of the results and stability of the calculation process.Wrong choice of the numerical schemes might lead to non-uniform and punctate pressure patterns on the walls, the displacement of the shock wave from its correct (experimental) position, and in most cases leads to the solution divergence.The most common phenomenon due to an inappropriate choice of a numerical scheme is a smooth pressure change, i.e. non-instantaneous change through the shock.
The study of the convection terms deserves most of the attention.In OpenFOAM applications, the convection term is commonly identified as div(phi; <field>), where phi refers to the flux ϕ = ρU.The <field> can be a vector or a scalar field.The typical choice of discretization in OpenFOAM is the Gauss scheme, and it requires a specification of interpolation scheme for the dependent field, i.e. <field>.For the simulation performed in this work, only improved version of upwind convection scheme with linear upwind differencing provided the reasonable results.There are also improved versions of some of the limited schemes for vector fields, in which the limit is formulated in a way to take into account the direction of the field [18].Hence, for the convection term which in the <field> section has a vector field, i.e. velocity U, the improved version of the limited secondorder bounded scheme had to be used.Upwind schemes [20] depend on flow direction, they are less dissipative than 1st order schemes, and there is great possibility of using the limiters.Upwind schemes seem to have gained, at least for the moment, much more popularity on unstructured meshes than the central schemes [21].
For the terms where a scalar field stands in the <field> section, the TVD (Total Variation Diminshing) scheme is employed, i.e. strictly bounded limited linear differencing.If the bounding was not applied, the solution would diverge.Some TVD/NVD (Normalized Variation Diagram) schemes [22,23,24,25] require a coefficient ψ, 0 ≤ ψ ≤ 1 where ψ = 1 corresponds to TVD conformance, usually giving best convergence and 0 ψ = corresponds to best accuracy.The coefficient ψ is usually given after chosen TVD/NVD interpolation scheme followed by the flux field.As the mesh used here appears to be relatively coarse (see Figs. 1 and 2), the best convergence was the primary objective.
The strength of the finite volume scheme is its applicability on unstructured meshes [26] where the rather heavy construction of data structures is compensated by the simpler and possibly automatic construction of meshes around complex geometries and an efficient implementation of the scheme (see [27,28,29,30,31,32]). The broader view of finite volume approximations of the Laplacian can be found in [32,33].In the most Laplacian terms in the equations, linear interpolation is chosen with explicit non-orthogonal correction for the surface normal gradients.
The discretization scheme for gradient terms can completely be specified by choosing least squares, or if the standard finite volume discretization of Gaussian integration [34,35] is chosen, it is necessary to perform interpolation of values from cell centres to face centres.It would be extremely unusual to select anything other than general interpolation schemes and in most cases the linear scheme is an effective choice, as was confirmed here.
It is well known that high-order spatial operators are much stiffer than lower-order ones.For time accurate problems, the allowable maximum of Courant number (Co) decreases with increasing the order of accuracy for explicit schemes.For viscous problems with highly clustered meshes to resolve the boundary layer, explicit high-order methods are severely limited by the time step size, and usually not competitive against low-order implicit methods in terms of efficiency [36].Here, for time derivatives the Euler scheme is used.Euler implicit discretization is the first order accurate in time, guarantees boundedness and is unconditionally stable [37].
The solver tolerance for pressure field was set to 1×e−08 and for the velocity and density fields it was set to 1×e−05.Relative tolerance is set to 0 due to the transient simulation, and so the solution was forced to converge to the solver tolerance in each time step.Under-relaxation technique was also used for improving stability of the computation.An optimum choice of relaxation factor α is a matter of compromiseit should be small enough to ensure stable computation, but also large enough to lead the iterative process towards the convergence quickly.
Pressure-implicit split-operator (PISO) algorithm [38] for solving equations for velocity and pressure was used since a transient problem is considered.This algorithm evaluates initial solution and then corrects further solutions.The number of corrections is set to 2. An additional correction used to account for mesh nonorthogonality is also applied, and the number of these corrections is set to 2. Numerical methods which have been employed during the calculation are summarized in Tab. 1.

The Onera M6 wing test case
In 1972, the ONERA Aerodynamic Department designed a swept back wing equipped with measuring systems, and it was used as an experimental support for basic studies of three-dimensional flows at high Reynolds numbers, from low to transonic speeds [10].
Wind tunnel data from this model, called the M6wing, present a good base both for computer program assessments, and for understanding of various flow phenomena, such as shock wave-boundary layer interaction and flow separation.
The selected data set was obtained in the ONERA S2MA wind tunnel at Much numbers of 0,7; 0,84; 0,88 and 0,92, for angles of attack up to 6° and for Reynolds numbers of about 12 million, computed with respect to the mean aerodynamic chord c [10].Solutions for the angle of attack of 3,06° and Mach number of 0,8395, have been computed and analysed in this paper.

Discussion of numerically obtained results
Figure 1 The RANS-LES mesh, y/b = 0 In this chapter, the obtained numerical solutions are presented and compared with the existing experimental results.Mesh used for numerical simulations is shown in Figs. 1 and 2, for cross-section y/b = 0 (i.e. the wing root, where y/b denotes relative span position on the wing model).This is a relatively coarse 3D mesh, which consists of 2,441,959 tetrahedrons.There are two main reasons for using such a mesh.First is an aim to obtain the best possible numerical solutions in aerodynamic sense on an unstructured and relatively coarse mesh, and the second is intention to apply similar meshes (not too complex but reliable) within much more complex, multidisciplinary fluidstructure interaction problem analyses, planned as future work.
In order to capture a specific range of turbulent scales, and enable determination of all relevant fluid parameters in the wake, it was necessary to generate a mesh that is appropriate for the combined RANS-LES method.If the typical RANS mesh was applied in the whole domain and used for the calculation with the proposed turbulence model, the fluid flow parameters would develop correctly in the vicinity of the wall (i.e.airfoil surface), but this would not allow those parameters to develop fully adequately within wake behind the wing (unsteady regions), and thus normally back-influence the parameters close to the wall.On the other hand, if the full LES mesh was applied (fluid cells are nearly of the same size in all directions within the entire computational domain), the parameters would develop correctly in every mesh cell, but that would not be acceptable for more complex fluid-structure interaction applications in future work, due to extremely large computational requirements.As the proposed method consists of the two parts, RANS and LES, the applied mesh had to have the "LES cells" of nearly the same size in the unsteady regions within the wake, and also acceptably large "RANS cells" near the wall.Such cells present no problem considering the overall accuracy, because of the robustness and low sensitivity of the k-omegaSSTSAS model to the cell sizes in the wall region.Although the best results should be expected when the cell sizes are of the same order in all directions within the entire calculation domain, the requirements for the reduction of computer resources requirements is always mandatory in operational engineering practice.In both calculation cases -by OpenFOAM code, custom-written for here presented analyses, and the ANSYS Fluent (commercial software, used for additional verifications), the normal shock wave, which is typical for the Onera M6 wing at these flow conditions is well defined, as can be seen in Figs. 3 and 4, as a quite sharp end of the red (local supersonic) zone.Example is given for relative semi-span position close to the tip, at y/b = 0,95.Also, another -front shock, in wing domain y/b ≈ 0,2 ÷ 0,9 is clearly captured (sharp transition from darkblue to light-blue zone), as can be seen on the wing top views in Figs. 5 and 6, where static pressure distribution is shown.In the sense of the pressure coefficient C p chordwise distributions, these shock waves can be identified in Figs.7 to 10, as sudden pressure coefficient jumps on upper wing surface.The two shocks are merged in one from y/b = 0,9 to the wing tip, as can be seen in Figs.11 and 12 (and in Figs. 5 and 6; also see [11,39]).Also, certain discrepancies between numerically obtained C p values behind the second shock, both by OpenFOAM and Fluent, compared with the experiment, can be seen on the presented diagrams.This could be improved either by increasing the here applied static mesh fineness, or by applying a dynamic mesh, capable of capturing the shock wave positions.The inclusion of realistic surface roughness could also give smoother C p transition just behind the second shock wave, but insisting on "perfection" in all those senses would certainly gradually increase the hardware resources requirements.Historically speaking, some of the best results for the Onera M6 test case were obtained using the Euler equations, applied on adaptive meshes [40].On the other hand, the planned application of the proposed method for the fluid-structure interaction problems of complete aircraft configurations requires the application of the Navier-Stokes equations with full turbulence modelling, where the Onera M6 wing has been analysed as a "pilot", or calibration case in this sense.
Specially, when high lift configurations are considered, the range of turbulent scales becomes enormous, and the full LES application would require huge computer resources.As a solution to this problem, a hybrid RANS-LES method, as the one demonstrated in this paper, can be an option.Also, in order to capture full range of the turbulent scales, an appropriate mesh must be generated, for example according to the recommendations given in [41].

Conclusion
In this paper the methods of numerical modelling and solutions of the turbulent flow over Onera M6 wing for angle of attack of 3,06° and Mach number of 0,8395, are presented.Completely unstructured mesh has been applied, with acceptably coarse elements in the vicinity of the wing.Nevertheless, the proposed method has confirmed its robustness, and ability to deliver the results of fair and satisfactory level of accuracy for the intended future applications, in the domain of much more complex fluid-structure interaction analyses of entire air vehicle configurations.
In this work, the compressible k-omegaSSTSAS turbulence model has been implemented in OpenFOAM for the first time.This model has shown slightly better abilities of capturing the shock wave positions, compared to the existing experimental results, than the commercial ANSYS Fluent package, used in this paper for additional verifications.
Calculations by both solvers were performed on the same mesh.It can be readily expected that the accuracy of all obtained results could have been improved by increasing the fineness and complexity of the mesh.On the other hand, if a mesh of such overall level of complexity should be applied for the analyses of the entire aircraft lofts in transonic speed domain, it should be done with special care, and as a very rational trade-off between the desired level of accuracy, and available hardware resources.

Figure 2
Figure 2The complete mesh-developing region behind the wing, y/b = 0