Comparative assessment of SAS,IDDES and hybrid filtering RANS/LES models based on second-moment closure

Abstract

The question of which turbulence model is better for a given class of applications is always confusing for the CFD researchers and users. Comparative assessments of scale-adaptive simulation (SAS), improved delay detached-eddy simulation (IDDES) and other hybrid RANS/LES models based on eddy-viscosity models (EVMs) are thoroughly investigated. But how well they perform based on a second-moment closure needs to be answered. In this paper, a widely acclaimed Reynolds-stress model (RSM) in aeronautical engineering, SSG/LRR- $\omega$ model, is carried out. The relevant test cases include the NACA0012 airfoil stalled flows and turret separated flow. In order to make a more reasonable comparison, a seventh-order scheme WCNS-E8T7 is adopted to reduce the influence of the numerical dissipation and a symmetrical conservative metric method is used to ensure the robustness. By comparing with the relevant experimental data and the solutions by original SSG/LRR- $\omega$ model (etc. URANS), it shows that all of the three hybrid methods (SAS, IDDES and hybrid filtering methods) based on the SSG/LRR- $\omega$ model have a good ability to simulate unsteady turbulence. Among them, the IDDES correction has the most potential.

Introduction

For computational fluid dynamics (CFD), the predicted result of arbitrary flow simulation depends on the appropriation of the underlying representation of flow physics and the accuracy of the numerical method solving the corresponding equations. ¹ To accurately predict the flow physics occupied by turbulence. The most perfect way is to directly resolve all turbulent fluctuations ranging from integral scale to Kolmogorov one (i.e. Direct Numerical Simulation, DNS), and the ideology of resolving large scale and modelling small scale can also be an alternative way (i.e. Large-eddy Simulation, LES). However, at a realistic flight Reynolds number, the turbulence scale in attached boundary layer is usually very small, which leads the computational cost too large to be afforded in current situation. Therefore, turbulence models based on the Reynolds-averaged Navier-Stokes (RANS) equation are still the backbone of industrial applied CFD methods. ²

Hybrid RANS/LES approaches, which aim to combine the superior accuracy of LES in the detached region with the efficiency of RANS in attached boundary layers, have drawn much attention, especially in separated flow in the past two decades. ³ Various hybrid methods have emerged, such as detached eddy simulation (DES),^4,5 scale adaptive simulation (SAS), ⁶ hybrid filtering method,^7,8 and others.

DES is the most popular one of Hybrid RANS/LES methods, which was proposed by Spalart ⁴ based on the one-equation SA model and then was extended to the two-equation SST model by Strelets. ⁹ Due to its simplicity and good predictions in massively separated flows, DES has gained a wide attention. However, as the application went deep, researchers found that the DES limiter can be falsely activated by grid refinement inside attached boundary layers, which is called Grid-Induced Separation (GIS). ⁵ In order to avoid this, the DES concept was extended to Delayed-DES (DDES), ¹⁰ following the proposal of Menter et al. ¹¹ to “shield” the boundary layer from the DES limiter. However like the DES, DDES approach is still only suitable for computing of massively separated turbulent flows. Shur et al. ¹² proposed a novel DES model, called Improved-DDES or IDDES. It has two branches, DDES and WMLES, including a set of empirical functions of subgrid length-scales designed to achieve good performance from these branches themselves and their coupling. By switching the activation of RANS and LES in different flow regions, IDDES significantly expands the scope of application of DDES with a well-balanced and powerful numerical approach to complex turbulent flows at high Reynolds number. ¹³

SAS is proposed by implanting the second derivative of velocity into turbulence scale equation. The first version is based on the one-equation KE1E model ¹⁴ with von Karman length scale to adapt to the underlying turbulent structures. It is a great attempt to achieve Hybrid RANS/LES performance without explicit grid dependency. Afterwards, a two-equation SAS model was proposed by Menter and Egorov ⁶ through reformulating Rotta’s equation and then transferred to SST model. One of the features of SAS model is that the limiter does not affect RANS behavior of the model. In other words, even the grid scale is close to Kolmogorov scale, the result cannot converge to DNS solution. That is the reason why some researchers classified SAS model as an advanced URANS method. ¹⁵

The hybrid filtering approach is regarded as one of the most rigorous hybrid methods. ¹⁶ The theoretical framework is based on the similarity of mathematical formation between Reynolds-average and filtered Navier-Stokes equations, ¹⁷ and then revised for compressible flow. ⁷ The hybrid filtering method includes two critical factors: turbulence model (including RANS model and LES model) and blending function. Currently, more and more advanced RANS models are used in hybrid way,^11,18–20 with a large number of blending functions proposed and applied.^21–23 Thus the hybrid filtering approach may also be a suitable choice for industrial application.

Above hybrid RANS/LES methods have their certain advantages, but from an engineering standpoint, the question of which model is better for a given class of applications is always confusing for the CFD users. In the EVM framework, comparative assessments are quite enough. ²⁴ But how well they perform based on a second-moment closure needs to be answered. As pointed out in NASA’s CFD Vision 2030 Study Report ²⁵ : the prediction of any separation that is initiated in the boundary layer will still require improvements in RANS-based methods. Among the RANS turbulence models, Reynolds-stress models are perceived as the most advanced ones and in principle will be potential in capturing the flow separation for a wider range of flows. For instance, Chaouat ²⁶ combined a subgrid-scale stress model with a partial integral transport model (PITM). Probst et al. ²⁷ introduced the idea of DES into the ${\epsilon }^h$ -RSM model and demonstrated the advantage of it over the DES based on an EVM in the simulation of small separation of airfoil wake. Wang et al. ²⁸ introduced IDDES correction into SSG/LRR- $\omega$ model and showed good results in simulations of airfoil stall and delta wing. The SSG/LRR- $\omega$ model, proposed by Eisfeld and Brodersen, ²⁹ is a combination of the Speziale-Sarkar-Gatski (SSG) model ³⁰ with the Launder-Reece-Rodi model ³¹ toward the wall, where the length scale is supplied by Menters baseline $\omega$ equation. In this paper, this widely acclaimed RSM is also carried out. Meanwhile, three kinds of hybrid RANS/LES methods are constructed, which are named SSG/LRR-IDDES, SSG/LRR-SAS and hybrid SSG/LRR, respectively.

Besides physical model, numerical scheme is also an important aspect when considering turbulent structures accurately resolved. ³² Although the numerical error can be mitigated by refining computational grid scale in principle when using low-order numerical scheme, it will certainly make the computation unaffordable. ³³ Therefore, we consider high-order scheme, which has more potential in delivering higher accuracy with less computational costs, ³⁴ and could be recommended to LES or Hybrid RANS/LES. ³⁵ Meanwhile, the coupling of turbulence model and numerical scheme is of great significance and deserve profound discussion. In present study, a seventh-order weighed compact nonlinear scheme (WCNS)^36,37 is applied to the comparative assessment of IDDES, SAS and hybrid filtering method based on the SSG/LRR- $\omega$ model. The relevant test cases include the NACA0012 airfoil stalled flow and turret separated flow. Among these simulations, a symmetrical conservative metric method (SCMM)^38,39 is adopted to ensure the accuracy and robustness of these high-order finite difference methods (FDMs) in the curvilinear coordinates.

The framework of this paper is organized as follows. In Section 2, we present a brief description of the turbulence models, including SSG/LRR- $\omega$ model, SSG/LRR-IDDES model, SSG/LRR-SAS model and Hybrid SSG/LRR model. In subsequent section, the numerical methods, including the WCNS and the SCMM, are introduced. Some numerical results are provided In Section 4, to make a comparative assessment of these turbulence models. Finally, some conclusions and remarks are drawn in Section 5.

Turbulence model

Original SSG/LRR- $\omega$ model

The reference for the standard implementation of the SSG/LRR- $\omega$ model is in NASA Turbulence Modeling Resource (TMR) website. ⁴⁰ A version with simple diffusion (called SSG/LRR-RSM-w2012-SD in the TMR website) is adopted in this paper, which was proposed by Eisfeld ⁴¹ first. The seven turbulence equations (written in conservation form) is given by the following:

$$\frac{\partial \overline{\rho }{\widehat{R}}_{ij}}{\partial t}+\frac{\partial \left(\overline{\rho }{\widehat{u}}_k{\widehat{R}}_{ij}\right)}{\partial x_k}=\overline{\rho }{P}_{ij}+\overline{\rho }{\Pi }_{ij}-\overline{\rho }{\epsilon }_{ij}+\overline{\rho }{D}_{ij}$$

(1)

$$\begin{array}{l}\frac{\partial \left(\overline{\rho }\omega \right)}{\partial t}+\frac{\partial \left(\overline{\rho }{\widehat{u}}_k\omega \right)}{\partial x_k}=\frac{{\alpha }_{\omega }\omega }{\widehat{k}}\frac{\overline{\rho }{P}_{kk}}{2}-{\beta }_{\omega }\overline{\rho }{\omega }^2\\ +\frac{\partial }{\partial x_k}\left[\left(\overline{\mu }+{\sigma }_{\omega }\frac{\overline{\rho }k}{\omega }\right)\frac{\partial \omega }{\partial x_k}\right]+{\sigma }_d\frac{\overline{\rho }}{\omega }\max \left(\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j},0\right)\end{array}$$

(2)

The Reynolds stresses are solved directly by $\overline{\rho }{\hat{R}}_{\mathrm{ij}}=-{\tau }_{\mathrm{ij}}=\overline{\rho {u\prime }_i{u\prime }_j}$ , where the hat represents the Favre average ( $\hat{\varphi }=\overline{\rho \varphi }/\overline{\rho }$ ) and $\overline{\rho }$ is mean density.The production term in equation (1) is exactly given by

$$\overline{\rho }{P}_{ij}=-\overline{\rho }{\widehat{R}}_{ik}\frac{\partial {\widehat{u}}_j}{\partial x_k}-\overline{\rho }{\widehat{R}}_{jk}\frac{\partial {\widehat{u}}_i}{\partial x_k}$$

(3)

The dissipation term is modeled by

$$\overline{\rho }{\epsilon }_{ij}=\frac{2}{3}\overline{\rho }\epsilon {\delta }_{ij}$$

(4)

where $\epsilon =C_{\mu }k\omega$ , $C_{\mu }=0.09$ and $k={\hat{R}}_{\mathrm{ii}}/2$ . The pressure-strain correlation is given by combination of the SSG model and the LRR model, as

$$\begin{array}{l}\overline{\rho }{\Pi }_{ij}=-\left(C_1\overline{\rho }\epsilon +\frac{1}{2}{C}_1^{\ast }\overline{\rho }{P}_{kk}\right){\widehat{a}}_{ij}\\ +C_2\overline{\rho }\epsilon \left({\widehat{a}}_{ik}{\widehat{a}}_{kj}-\frac{1}{3}{\widehat{a}}_{kl}{\widehat{a}}_{kl}{\delta }_{ij}\right)\\ +\left(C_3-C_3^{\ast }\sqrt{{\widehat{a}}_{kl}{\widehat{a}}_{kl}}\right)\overline{\rho }k{{\widehat{S}}^{\ast }}_{ij}\\ +C_4\overline{\rho }k{\left({\widehat{a}}_{ik}{\widehat{S}}_{jk}+{\widehat{a}}_{jk}{\widehat{S}}_{ik}-\frac{2}{3}{\widehat{a}}_{kl}{\widehat{S}}_{kl}{\delta }_{ij}\right)}_{ij}\\ +C_5\overline{\rho }k\left({\widehat{a}}_{ik}{\widehat{W}}_{jk}+{\widehat{a}}_{jk}{\widehat{W}}_{ik}\right)\end{array}$$

(5)

where pressure dilatation is neglected, and the anisotropy tensor is ${\hat{a}}_{\mathrm{ij}}={\hat{R}}_{\mathrm{ij}}/k-2{\delta }_{\mathrm{ij}}/3$ . Also,

$$\begin{array}{l}{\widehat{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\widehat{u}}_i}{\partial x_j}+\frac{\partial {\widehat{u}}_j}{\partial x_i}\right),\\ {{\widehat{S}}^{\ast }}_{ij}={\widehat{S}}_{ij}-\frac{1}{3}{\widehat{S}}_{kk}{\delta }_{ij},\\ {\widehat{W}}_{ij}=\frac{1}{2}\left(\frac{\partial {\widehat{u}}_i}{\partial x_j}-\frac{\partial {\widehat{u}}_j}{\partial x_i}\right)\end{array}$$

(6)

The diffusion term is modeled via:

$$\overline{\rho }{D}_{ij}=\frac{\partial }{\partial x_k}\left[\left(\overline{\mu }+\frac{D}{C_{\mu }}{\mu }_t\right)\frac{\partial {\widehat{R}}_{ij}}{\partial x_k}\right]$$

(7)

where ${\mu }_t$ is obtained by ${\mu }_t=\overline{\rho }k/\omega$ .

This model follows Menter’s ⁴² approach, combining the SSG- $\epsilon$ model at the outer edge of the boundary layer with Wilcox’s ⁴³ LRR- $\omega$ model near the wall by blending the coefficients,

$$\varphi =F_1{\varphi }^{(\mathrm{LRR}-\omega )}+(1-F_1){\varphi }^{(\mathrm{SSG}-\epsilon )}$$

(8)

The bounding values of the $\omega$ equation are listed in Table 1, while the Reynolds-stress equations are listed in Table 2.

Table 1.

Bounding values of $\omega$ equation coefficients.

	${\alpha }_{\omega }$	${\beta }_{\omega }$	${\sigma }_{\omega }$	${\sigma }_d$
${\varphi }^{\left(\epsilon \right)}$	0.44	0.0828	0.856	1.712
${\varphi }^{\left(\omega \right)}$	0.5556	0.075	0.5	0

Table 2.

Bounding values of the Reynolds-stress equation coefficients, $C_2^{\left(\mathrm{LRR}\right)}=0.52$ .

	${\hat{C}}_1$	${{\hat{C}}_1}^{*}$	${\hat{C}}_2$	${\hat{C}}_3$	${{\hat{C}}_3}^{*}$	$D$
${\varphi }^{\left(\mathrm{SSG}\right)}$	1.7	0.9	1.05	0.8	0.65	0.22 × 2/3
${\varphi }^{\left(\mathrm{LRR}\right)}$	1.8	0	0	0.8	0	$0.5C_{\mu }$
		${\hat{C}}_4$			${\hat{C}}_5$
${\varphi }^{\left(\mathrm{SSG}\right)}$			0.625		0.2
${\varphi }^{\left(\mathrm{LRR}\right)}$			$(9C_2^{\left(\mathrm{LRR}\right)}+6)/11$		$(-7C_2^{\left(\mathrm{LRR}\right)}+10)/11$

SSG/LRR-IDDES method

Detached eddy simulation was first constructed on one-equation SA model. ⁴ According to the grid resolution, it can switch between RANS and LES mode. Afterwards, a similar formulation is applied to two-equation SST model, ⁹ which based on a criterion between local grid scale $L_g=max(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)$ and turbulence length scale $L_{t,\mathrm{DES}}=\sqrt{k}/\left(C_{\mu }\omega \right)$ .

The original DES model acquires widespread acceptance in industrial applications, but there are still some deficiencies have been found, such as model stress depletion (MSD) ¹⁰ problem and grid induced separation (GIS) phenomenon. ⁵ In order to solve these problems, many advanced revisions were proposed, such as Delayed-DES(DDES) model ¹⁰ and Improved-DDES (IDDES) model. ¹² The IDDES model is constructed by coupling wall-modeled LES (WMLES) and DDES to eliminate the log layer mismatch (LLM) problem, and maintain the compatibility of general DES model. The major improvement of IDDES model mainly reflects in the near-wall modification of LES filter, and more rapid transition between RANS and LES than original DES model.

Referring Shur et al.’s ¹² work, the detailed IDDES formulation for the SSG/LRR- $\omega$ model is made by modifying the isotropic dissipation $\epsilon$ in equation (1),

$$\epsilon =k^{3/2}/L_{\mathrm{IDDES}}$$

(9)

Noted that the $\epsilon$ does not only appear in the dissipation term, but also in the pressure-strain term. The limiter $L_{\mathrm{IDDES}}$ is

$$L_{\mathrm{IDDES}}={\hat{f}}_d(1+f_e)L_{t,\mathrm{DES}}+(1-{\hat{f}}_d)C_{\mathrm{DES}}{L}_g,$$

(10)

where the blending function ${\hat{f}}_d$ is defined as

$${\hat{f}}_d=max((1-f_d),f_B).$$

(11)

The $f_B$ is the empirical function in WMLES part, which is solved by

$$f_B=min(2\exp (-9{\alpha }^2),1.0)$$

(12)

where

$$\alpha =0.25-d_w/max(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)$$

(13)

$f_d$ is the empirical function in DDES, which defined as

$$f_d=1-\tanh {\left(8r_d\right)}^3$$

(14)

with

$$r_d=(\nu +k/\omega )/\left(\sqrt{\frac{\partial {\hat{u}}_i}{\partial x_j}\frac{\partial {\hat{u}}_i}{\partial x_j}}{\kappa }^2{d_w}^2\right)$$

(15)

$d_w$ is the distance to the nearest wall. $C_{\mathrm{DES}}$ is LES length scale coefficient and obtained by the shielding function $C_{\mathrm{DES}}=(1-F_1)C_{\mathrm{DES}}^{\left(\epsilon \right)}+F_1{C}_{\mathrm{DES}}^{\left(\omega \right)}$ , same as that in SST-DES model. Slight difference is existed between $C_{\mathrm{DES}}^{\left(\epsilon \right)}$ and $C_{\mathrm{DES}}^{\left(\omega \right)}$ for different codes. In this paper, $C_{\mathrm{DES}}^{\left(\epsilon \right)}$ is set to 0.61 and $C_{\mathrm{DES}}^{\left(\omega \right)}$ equals to 0.78.

Elevating function $f_e$ is another empirical function involved in equation (10). It can avoid excessive reduction of Reynolds stress, observed in the sandwich of RANS and LES regions. The elevating function is intended to eliminate the LLM problem.

$$f_e=max((f_{e1}-1),0)f_{e2},$$

(16)

where the function $f_{e1}$ is defined as

$$f_{e1}=\{\begin{array}{c}2\exp (-11.09{\alpha }^2),\alpha \ge 0\hfill \\ \multicolumn{1}{c}{2\exp (-9.0{\alpha }^2),\alpha <0}\end{array},$$

(17)

and the function $f_{e2}$ is

$$f_{e2}=1.0-max(\tanh \left({\left(c_t^2{r}_{\mathrm{dt}}\right)}^3\right),\tanh \left({\left(c_l^2{r}_{\mathrm{dl}}\right)}^{10}\right)).$$

(18)

where the quantities $r_{\mathrm{dt}}$ and $r_{\mathrm{dl}}$ are the “turbulent” and “laminar” analogues of $r_d$ (equation (15)).

In addition, a modified grid size length scale was adopted by Shur et al., ¹² which is used here:

$$L_g=min(max(C_w{d}_w,C_w{h}_{max},h_{\mathrm{wn}}),h_{max}).$$

(19)

Where, $h_{max}=max(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z)$ is the original grid length scale in DES model, $h_{\mathrm{wn}}$ is the grid step in wall-normal direction. Besides, all of the empirical constants are given here: $\kappa =0.41$ , $C_w=0.15$ , ${\mathrm{c}}_t=1.87$ and ${\mathrm{c}}_l=5.0$ . Finally, it should be point out that the equations (11)–(19) are same with those in SST-IDDES, however, the empirical coefficientes need to be recalibrated in the RSM framework.

SSG/LRR-SAS method

Scale adaptive simulation is proposed by implementing second derivative of velocity into turbulence scale equation. The derivation is based on the theory of Rotta, ⁴⁴ and then introduced to two-equation SST model. ⁶ Unlike the hybrid framework of DES model, SAS model bases on a criterion between von Karman length scale $L_{\mathrm{vK}}=\kappa \sqrt{2{\hat{S}}_{\mathrm{ij}}{\hat{S}}_{\mathrm{ij}}}/\sqrt{{\left({\nabla }^2\hat{u}\right)}^2+{\left({\nabla }^2\hat{v}\right)}^2+{\left({\nabla }^2\hat{w}\right)}^2}$ and turbulence length scale $L_{t,\mathrm{SAS}}=\sqrt{k}/\left({c_{\mu }}^{1/4}\omega \right)$ .

Considering the similarity of scale-eqution between SSG/LRR- $\omega$ and SST models, the SAS bases on SSG/LRR- $\omega$ model differs from the original RANS formulation only by the additional SAS source term $Q_{\mathrm{SAS}}$ in the $\omega$ -equation.

$$\begin{array}{c}{Q}_{\mathrm{SAS}}=max[\overline{\rho }{\zeta }_2\kappa {\hat{S}}^2{\left(\frac{L_{t,\mathrm{SAS}}}{L_{\mathrm{vK}}}\right)}^2\hfill \\ \multicolumn{1}{c}{-C_{\mathrm{SAS}}\frac{2\overline{\rho }k}{{\sigma }_{\Phi }}max(\frac{1}{k^2}\frac{\partial k}{\partial x_j}\frac{\partial k}{\partial x_j},\frac{1}{{\omega }^2}\frac{\partial \omega }{\partial x_j}\frac{\partial \omega }{\partial x_j}),0.0]}\end{array}$$

(20)

where $\hat{S}=\sqrt{2{\hat{S}}_{\mathrm{ij}}{\hat{S}}_{\mathrm{ij}}}$ .

In addition, high wave number is damped to avoid accumulating energy of small scales. It can be realized by imposing a lower limit of the von Karman length scale ⁶ :

$$\begin{array}{c}{L}_{\mathrm{vK}}=max(\frac{\kappa \sqrt{2{\hat{S}}_{\mathrm{ij}}{\hat{S}}_{\mathrm{ij}}}}{\sqrt{{\left({\nabla }^2\hat{u}\right)}^2+{\left({\nabla }^2\hat{v}\right)}^2+{\left({\nabla }^2\hat{w}\right)}^2}}\hfill \\ \multicolumn{1}{c}{,C_s\sqrt{\frac{\kappa {\zeta }_2}{{\beta }_{\mathrm{SAS}}/C_{\mu }-{\alpha }_{\mathrm{SAS}}}}\mathrm{\Delta })}\end{array}$$

(21)

where $\mathrm{\Delta }=\sqrt[3]{V}$ and $V$ is the volume of the grid cell. Meanwhile, all of the empirical constants are given here: ${\zeta }_2=3.51$ , ${\sigma }_{\Phi }=2/3$ , $C_{\mathrm{SAS}}=2$ , $C_s=0.11$ , $\kappa =0.41$ , ${\alpha }_{\mathrm{SAS}}=0.44$ and ${\beta }_{\mathrm{SAS}}=0.0828$ .

Hybrid Reynolds-stresses and subgrid-stresses method

The proposed hybrid filter is a linear combination of URANS and LES operators:

$${〈\varphi 〉}_{\mathrm{Hybrid}}=F{〈\varphi 〉}_{\mathrm{RANS}}+(1-F){〈\varphi 〉}_{\mathrm{LES}}$$

(22)

where $F$ is a normalized blending function ( $0\le F\le 1$ ), it can be both spatially and temporally dependent. The hybrid operator defined in equation (22) is applied to the original Navier-Stokes equations, and generate extra equations for hybrid variables. These extra terms were found to ensure a proper transfer of momentum in RANS to LES transition region, hence avoid incorporating artificial stochastic forcing. Due to their relatively complex formulation, these extra terms have always been neglected to date. The simplified hybrid filter has performed over a wide range of turbulent configurations with a satisfied performance. ¹⁶ Meanwhile, the method was extended to compressible flow as well. ⁷

Seven partial different equations for the Reynolds stresses ${\hat{R}}_{\mathrm{ij}}$ , specific dissipation rate $\omega$ are required in the SSG/LRR- $\omega$ model, whereas only six equations for the subgrid-scale (or subfilter-scale) stresses (SGS) ${\hat{R}}_{\mathrm{ij}}^{\mathrm{sgs}}$ is employed in the LES. Combining these within the hybrid RANS/LES framework, a new hybrid turbulence-stress term ${\hat{R}}_{\mathrm{ij}}^{\mathrm{hyb}}$ can be defined as follows:

$${\hat{R}}_{\mathrm{ij}}^{\mathrm{hyb}}=F{\hat{R}}_{\mathrm{ij}}^{\mathrm{rans}}+(1-F){\hat{R}}_{\mathrm{ij}}^{\mathrm{sgs}}$$

(23)

Applying the hybrid filtering approach to the turbulence models, new transport equations can be derived for the hybrid turbulence-stress ${\hat{R}}_{\mathrm{ij}}^{\mathrm{hyb}}$ :

$$\begin{array}{c}\frac{\partial \overline{\rho }{\hat{R}}_{\mathrm{ij}}^{\mathrm{hyb}}}{\partial t}+\frac{\partial \left(\overline{\rho }{\hat{u}}_k{\hat{R}}_{\mathrm{ij}}^{\mathrm{hyb}}\right)}{\partial x_k}=\overline{\rho }[F({\Pi }_{\mathrm{ij}}^{\mathrm{rans}}-{\epsilon }_{\mathrm{ij}}^{\mathrm{rans}}+D_{\mathrm{ij}}^{\mathrm{rans}})\hfill \\ +(1-F)({\Pi }_{\mathrm{ij}}^{\mathrm{sgs}}-{\epsilon }_{\mathrm{ij}}^{\mathrm{sgs}}+D_{\mathrm{ij}}^{\mathrm{sgs}})]+\overline{\rho }{P}_{\mathrm{ij}}^{\mathrm{hyb}}\hfill \end{array}$$

(24)

The URANS and SGS eddy viscosities are formulated as:

$$\begin{array}{c}{\mu }_{\mathrm{rans}}=\frac{\rho k^{\mathrm{hyb}}}{\omega }\hfill \\ {\mu }_{\mathrm{sgs}}=\rho C_{\nu }\mathrm{\Delta }\sqrt{k^{\mathrm{hyb}}}\hfill \\ k^{\mathrm{hyb}}={\hat{R}}_{\mathrm{ii}}^{\mathrm{hyb}}/2\hfill \end{array}$$

(25)

The modeling of the dissipation-rate in SGS model is not given by its transport equation as in RANS models but it is then explicitly computed by means of the grid size $\mathrm{\Delta }$ as

$$\overline{\rho }{{\epsilon }_{\mathrm{ij}}}^{\mathrm{sgs}}=\frac{2}{3}\overline{\rho }{\epsilon }_{\mathrm{sgs}}{\delta }_{\mathrm{ij}}$$

(26)

where ${\epsilon }_{\mathrm{sgs}}=C_{\epsilon }{k}_{\mathrm{sgs}}^{3/2}/\mathrm{\Delta }$ . The redistribution term proposed by Deardorff^45,46 is adopted in the SGS part:

$$\overline{\rho }{{\Pi }_{\mathrm{ij}}}^{\mathrm{sgs}}=-c_m\frac{k_{\mathrm{sgs}}^{1/2}}{\mathrm{\Delta }}(\overline{\rho }{\hat{R}}_{\mathrm{ij}}^{\mathrm{sgs}}-\frac{2}{3}{k}_{\mathrm{sgs}}{\delta }_{\mathrm{ij}})+\frac{2}{5}{k}_{\mathrm{sgs}}{\overline{S}}_{\mathrm{ij}}$$

(27)

The diffusion term is modeled by analogy with the RANS modeling as

$$\overline{\rho }{D_{\mathrm{ij}}}^{\mathrm{sgs}}=\frac{\partial }{\partial x_k}\left[c_m\overline{\rho }\mathrm{\Delta }{k}_{\mathrm{sgs}}^{1/2}\frac{\partial {\hat{R}}_{\mathrm{ij}}^{\mathrm{sgs}}}{\partial x_k}\right]$$

(28)

Finally, a proper blending function $F$ that will switch from the URANS model to the LES model when appropriate must be selected. In the present effort, the $F_2$ blending function initially proposed by Menter for the SST model has been successfully applied by Sanchez-Rocha and Menon ⁷ and is repeated here:

$$F=F_2=\tanh [max{(2\frac{\sqrt{k^{\mathrm{hyb}}}}{C_{\mu }\omega d_w},\frac{500\nu }{{d_w}^2\omega })}^2]$$

(29)

Numerical methods

Coordinate transformation and symmetrical conservative metric method

In Cartesian coordinate, the three-dimensional RANS equations are

$$\frac{\partial Q}{\partial t^{*}}+\frac{\partial E}{\partial x^{*}}+\frac{\partial F}{\partial y^{*}}+\frac{\partial G}{\partial z^{*}}=\frac{1}{\mathrm{Re}}(\frac{\partial E_v}{\partial x^{*}}+\frac{\partial F_v}{\partial y^{*}}+\frac{\partial G_v}{\partial z^{*}})$$

(30)

where the supscript “*” represents the non-dimensional variable, which yields $Q={[{\overline{\rho }}^{*},{\overline{\rho }}^{*}{\hat{u}}^{*},{\overline{\rho }}^{*}{\hat{v}}^{*},{\overline{\rho }}^{*}{\hat{w}}^{*},{\overline{\rho }}^{*}{\hat{e}}^{*}]}^T$ . $E$ , $F$ , $G$ and $E_v$ , $F_v$ , $G_v$ are inviscid and viscous flux, respectively. The above form can be found in many textbooks, ⁴³ where the supscript “*” is always omitted as following.

As the transformation from Cartesian coordinates to curvilinear coordinates is applied, that is,

$$x=x(\xi ,\eta ,\zeta ),y=y(\xi ,\eta ,\zeta ),z=z(\xi ,\eta ,\zeta )$$

(31)

equation (26) becomes

$$\frac{\partial \tilde{Q}}{\partial t}+\frac{\partial \tilde{E}}{\partial \xi }+\frac{\partial \tilde{F}}{\partial \eta }+\frac{\partial \tilde{G}}{\partial \zeta }=\frac{1}{\mathrm{Re}}(\frac{\partial {\tilde{E}}_v}{\partial \xi }+\frac{\partial {\tilde{F}}_v}{\partial \eta }+\frac{\partial {\tilde{G}}_v}{\partial \zeta })$$

(32)

where

$$\begin{array}{c}\tilde{Q}=\mathrm{JQ},\hfill \\ \tilde{E}={\tilde{\xi }}_{x}E+{\tilde{\xi }}_{y}F+{\tilde{\xi }}_{z}G,{\tilde{E}}_v={\tilde{\xi }}_x{E}_v+{\tilde{\xi }}_y{F}_v+{\tilde{\xi }}_z{G}_v\hfill \\ \tilde{F}={\tilde{\eta }}_{x}E+{\tilde{\eta }}_{y}F+{\tilde{\eta }}_{z}G,{\tilde{F}}_v={\tilde{\eta }}_x{E}_v+{\tilde{\eta }}_y{F}_v+{\tilde{\eta }}_z{G}_v\hfill \\ \tilde{G}={\tilde{\zeta }}_{x}E+{\tilde{\zeta }}_{y}F+{\tilde{\zeta }}_{z}G,{\tilde{G}}_v={\tilde{\zeta }}_x{E}_v+{\tilde{\zeta }}_y{F}_v+{\tilde{\zeta }}_z{G}_v\hfill \end{array}$$

(33)

$J$ is the Jacobian of grid transformation and ${\tilde{\xi }}_x$ , ${\tilde{\xi }}_y$ , ${\tilde{\xi }}_z$ , ${\tilde{\eta }}_x$ , ${\tilde{\eta }}_y$ , ${\tilde{\eta }}_z$ , ${\tilde{\zeta }}_x$ , ${\tilde{\zeta }}_y$ , ${\tilde{\zeta }}_z$ are the grid metrics, which are solved by using a symmetrical conservative metric method (SCMM)^38,39 to satisfy the geometric conservation law (GCL).^47,48

High-order weighted compact nonlinear scheme

After the weighted compact nonlinear scheme (WCNS) ⁴⁹ was proposed, a series of WCNS were developed and applied to a wide range of applications.^38,50–52 The WCNS scheme consists of three components: (i) cell-edge to cell-node central flux difference; (ii) flux evaluation at the cell-edge, and (iii) cell-node to cell-edge weighted nonlinear interpolation of flow variables. In this paper, a 7th-order tri-diagonal compact one WCNS-E8T7^36,37 is adopted for the comparative assessment of IDDES, SAS and hybrid iltering models based on the SSG/LRR- $\omega$ turbulence model.

Considering the WCNS-E8T7, the inviscid term in equation (28) is discretizated by an explicit eighth-order central flux differencing:

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial {\tilde{E}}_i}{\partial \xi }=\frac{1225}{1024\mathrm{\Delta }\xi }({\tilde{E}}_{i+1/2}-{\tilde{E}}_{i-1/2})}\\ \multicolumn{1}{c}{-\frac{245}{3072\mathrm{\Delta }\xi }({\tilde{E}}_{i+3/2}-{\tilde{E}}_{i-3/2})}\\ \multicolumn{1}{c}{+\frac{49}{5120\mathrm{\Delta }\xi }({\tilde{E}}_{i+5/2}-{\tilde{E}}_{i-5/2})}\\ \multicolumn{1}{c}{-\frac{5}{7168\mathrm{\Delta }\xi }({\tilde{E}}_{i+7/2}-{\tilde{E}}_{i-7/2})}\end{array}$$

(34)

where

$$\begin{array}{c}{\tilde{E}}_{i+1/2}=\tilde{E}[Q_{i+1/2}^L,Q_{i+1/2}^R,\hfill \\ \multicolumn{1}{c}{{\left({\tilde{\xi }}_x\right)}_{i+1/2},{\left({\tilde{\xi }}_y\right)}_{i+1/2},{\left({\tilde{\xi }}_z\right)}_{i+1/2}]}\end{array}$$

(35)

Only the discretization in $\xi$ direction is given here, and a similar procedure for other directions. Equation (31) is the flux at cell edges which can be evaluated by various flux schemes. In current paper, Roe’s ⁵³ flux difference scheme is used. $Q_{i+1/2}^L$ , $Q_{i+1/2}^R$ are cell-edge left and right variable values, which are calculated by weighted nonlinear interpolation methods, only the $Q_{i+1/2}^L$ is described here as an example.

For the Step (iii), a seventh-order weighted compact nonlinear interpolation is adopted as

$$\begin{array}{c}{\left(\begin{array}{c}\frac{1}{2}{\omega }_5\hfill \\ \multicolumn{1}{c}{{\omega }_0+{\omega }_1+{\omega }_2}\\ \multicolumn{1}{c}{\frac{3}{14}{\omega }_5}\end{array}\right)}^T\left(\begin{array}{c}{Q}_{j-1/2}^L\hfill \\ \multicolumn{1}{c}{Q_{j+1/2}^L}\\ \multicolumn{1}{c}{Q_{j+3/2}^L}\end{array}\right)\hfill \\ ={\left(\begin{array}{c}\frac{3}{8}{\omega }_0-\frac{25}{896}{\omega }_5\hfill \\ \multicolumn{1}{c}{\frac{1}{32}(-40{\omega }_0-4{\omega }_1+9{\omega }_5)}\\ \multicolumn{1}{c}{\frac{3}{64}(40{\omega }_0+16{\omega }_1+8{\omega }_2+5{\omega }_5)}\\ \multicolumn{1}{c}{\frac{1}{8}(3{\omega }_1+6{\omega }_2+\frac{5}{4}{\omega }_5)}\\ \multicolumn{1}{c}{\frac{1}{8}(-{\omega }_2+\frac{9}{16}{\omega }_5)}\end{array}\right)}^T\left(\begin{array}{c}{Q}_{j-2}\hfill \\ \multicolumn{1}{c}{Q_{j-1}}\\ \multicolumn{1}{c}{Q_j}\\ \multicolumn{1}{c}{Q_{j+1}}\\ \multicolumn{1}{c}{Q_{j+2}}\end{array}\right)\hfill \end{array}$$

(36)

The weights ${\omega }_k$ is given:

$${\omega }_k=D_k(1+\frac{{\tau }_7}{{\beta }_k+{\epsilon }_{\beta }}),k=0,1,2,5$$

(37)

where

$${\tau }_7=|{\beta }_5-\frac{{\beta }_0+4{\beta }_1+{\beta }_2}{6}|{\left(\frac{{\beta }_0-{\beta }_2}{{\beta }_0+{\beta }_2+{\epsilon }_{\beta }}\right)}^2$$

(38)

The ${\beta }_0,{\beta }_1$ and ${\beta }_2$ a are the smooth indicators, which was presented in Liu et al., ³⁶ Deng et al. ⁵⁰ Besides, $D_0=1/16,D_1=5/8,D_2=5/16$ and $D_5=1$ are the optimal weights to recover seventh-order interpolation.

Comparative results

NACA0012 airfoil stalled flows

Flow configuration

NACA0012 airfoil is one of the benchmark cases in aeronautical separated flows. There are plenty of experimental studies,^54,55 covering data from $0^{°}$ to ${90}^{°}$ angles of attack (AoA). The experimental flow parameters, needed to set up appropriate numerical simulations are $R{e}_c=1.3\times {10}^6$ and $M{a}_{\infty }=0.5$ . The chord length can be set as $c=1.0$ m.

Simulations of the NACA0012 airfoil at four different angle of attack of $5^{°}$ , ${17}^{°}$ , ${45}^{°}$ , ${60}^{°}$ are carried out. An O-type grid is generated, which contains $192\times 102\times 30$ cells in the streamwise, wall-normal, and spanwise directions, respectively. Figure 1 shows the grid on the wall and X-Y plane, whose size is equal to $100c$ . In the spanwise direction, all the grids are uniform and the span-size of the domain is set to $1.0c$ , which has been adopted in Yang and Zha ¹³ and proven to be adequate for three dimensional numerical studies.

Figure 1.

Computational grid in the X-Y plane and wall, NACA0012 airfoil.

All of the turbulence models described above (SSG/LRR-IDDES, SSG/LRR-SAS and Hybrid SSG/LRR models) are employed here and compared with the traditional SST-IDDES approach. The seventh-order compact scheme WCNS-E8T7 is carried out to reduce the effect of numerical dissipation. Time iterative is performed using the dual-time stepping technique with 0.01 times the non-dimensional time steps ( $\mathrm{\Delta }{\mathrm{t}}^{*}=0.01c/U_{\infty }$ ). It is reasonable to compute the time-averaged solution over an interval of 50 convective times ( $T=c/U_{\infty }$ ) after the initial transients were eliminated. The stationary state (in a time-averaged sense) is reached after 50 convective times; therefore, the simulation was performed until $100T$ , for a total of approximately 10,000 non-dimensional time steps.

Instantaneous flow field at ${45}^{°}$ angle of attack

At ${45}^{°}$ angle of attack, it is a stall flow with massive flow separation. Figure 2 presents the comparison of flow structures in the form of an iso-surface of the instantaneous Q-criterion. Only a two-dimensional form of the body-shedding vortex is predicted by the SST model without any scale correction. The SSG/LRR-URANS approach produces relatively lower eddy viscosity than the SST-URANS, and thus a little three-dimensional vortex can still be found near the wall. By contrast, all of the scale-resolving approaches reveal the shedding of vortical structures from the leading to trailing edges. The vortex behind the airfoil is highly chaotic, which conforms to the stochastic behavior of the turbulence flow.

Figure 2.

Iso-surface of the Q-criterion at 100T, colored by Mach number.

Figure 3 displays the instantaneous vorticity contours of 50% span at non-dimensional time 100T. Whether for SST or SSG/LRR- $\omega$ model, the URANS approach will capture the phase locked vortex shedding that usually occurs at laminar flow, while IDDES model obtain more realistic and turbulent vortical flow structures in the regions of massive separations. The SSG/LRR-SAS and hybrid SSG/LRR approaches can also simulate unsteady turbulence, but the preserved details are not rich enough compared with SSG/LRR-IDDES.

Figure 3.

Vorticity contours of 50% span at 100T.

Time-averaged result

The averaged lift and drag coefficients obtained by different scale-resolving approaches are listed in Figure 4. At the AoA of $5^{°}$ , the flow around the airfoil is approximately an attached turbulence (only minor separation occurs near the trailing edge). All of the preditions are naturally close.

Figure 4.

Time-averaged lift and drag coefficients of all NACA0012 computations and comparison with experimental data.

At the AoA of ${17}^{°}$ , the airfoil is in the critical region of stall, where the flow is complex. The SSG/LRR-IDDES method obtains a slightly lower lift coefficient result than the experimental data, while others each predict a higher lift coefficient.

At the AoA of ${45}^{°}$ and ${60}^{°}$ , the massive separated flows occur on the leeward side of the airfoil. The SST-URANS computation naturally over-predicts the lift and drag coefficients due to its too much physical dissipation, but the improvement by SSG/LRR-URANS is still limited. A scale-resolving approach is necessary at these high AoA with massive flow separations. Like the SST-IDDES, the SSG/LRR-SAS, SSG/LRR-IDDES and hybrid filtering models can successfully predict the lift and drag. Figure 5 compares the pressure coefficient predictions at the AoA of ${45}^{°}$ with the experimental data for a normal flat plate.^9,56 This further confirms the above results and trends.

Figure 5.

Time-averaged pressure coefficient distribution over NACA0012 airfoil at the AoA of ${45}^{°}$ .

Taken all cases together, the degree of coincidence among RSM-based approaches is SSG/LRR-IDDES > Hybrid SSG/LRR > SSG/LRR-SAS > SSG/LRR-URANS, compared with the experimental data. On the other hand, the results based on SSG/LRR-IDDES are slightly better than the SST-IDDES. It shows that the RSM has a good potential as the basic turbulence model in the hybrid RANS/LES framework.

Flow over a turret

Flow configuration

The turret is a common type of protrusion outside the aircraft. For example, the camera pod of the UAV (unmanned aerial vehicle), which always affects the aerodynamic characteristics of the belly of UAV. In present study, the turret geometry and flow conditions are selected to model an experiment performed at the U.S. Air Force Academy, ⁵⁷ see Figure 6. The turret consists of a half-foot radius hemisphere on top of a circular cylinder attached to the wind-tunnel wall. The flow conditions for the computation and experiment are a Reynolds number based on the turret diameter of $R{e}_D=2.4\times {10}^6$ and Mach number $M{a}_{\infty }=0.4$ .

Figure 6.

Computational wind-tunnel cross section. ⁵⁷

The numerical domain extended over $-3.75\le x\le 10.0$ in the streamwise direction, $-0.00\le y\le 3.0$ in the normal direction, and $-1.7083\le z\le 1.29167$ in the spanwise direction. The above settings refer to the Morgan and Visbal ⁵⁸ and is presented in Figure 7. A H-O type grid is generated, which contains $121\times 241\times 129$ in the body normal, body tangential, and circumferential directions, respectively. No-slip boundary conditions are used in the simulations for the turret surface and wind-tunnel floor, while inviscid boundary conditions are enforced for the top and side wind-tunnel walls.

Figure 7.

The numerical domain and surface mesh for turret and wind tunnel.

The SSG/LRR-IDDES, SSG/LRR-SAS and hybrid SSG/LRR models are employed here and compared with the SSG/LRR-URANS approach. At the same time, SST-URANS and SST-IDDES are carried out as references. The seventh-order scheme WCNS-E8T7 is adopted to reduce the effect of numerical dissipation. Time iterative is performed using the dual time stepping technique with the non-dimensional time steps $\mathrm{\Delta }{t}^{*}=0.005T_D$ ( $T_D=D/U_{\infty }$ ). A total of 10000 $\mathrm{\Delta }{t}^{*}$ needs to be implemented and the time-averaged statistics begins after 5000 $\mathrm{\Delta }{t}^{*}$ .

• Referring to Morgan and Visbal’s setting, the inflow profile was generated from an individual flat-plate simulation. Figure 8 shows the time-mean velocity profile. This ensures that all methods are coincident at inlet.

Figure 8.

Time-mean streamwise inflow velocity profile.

Instantaneous flow field

A 3D view of the whole flow region is given in Figure 9. This figure compares the instantaneous flow features at 40T predicted by the six approaches. It is predictable that the SST-URANS solutions are absent of multiple levels of fine scale structures. The SSG/LRR-URANS gives a relatively obvious improvement, but the retained flow details are still not enough. Based on SSG/LRR- $\omega$ model, SAS approach and IDDES model can effectively reduce dissipation in the separated zone. However for the hybrid filtering approach, there is no obvious improvement than SSG/LRR-URANS. It is indicated that Deardorff’s subgrid-stress model has too much physical viscosity.

Figure 9.

Iso-surface of the Q-criterion at 40T, colored by Mach number.

Time-averaged result

Next, we focus on the time-averaged solutions. Figure 10 compares the surface streamlines predicted by the six approaches to experimental oil-surface-flow visualizations (Figure 11; The copyrights of Figure 11 belong to the AIAA Journal and the relevant authors). In this view, we are looking from a perspective above and downstream of the turret. The surface-flow visualization from the experiment was enhanced by the original authors with schematic streamlines to highlight flow patterns. The flow on the wind tunnel floor is very similar to that seen in the experiment. All solutions display very similar features in the region ahead of the turret influenced by the necklace vortex. However in the turret wake, especially in the backside surface of turret, the difference appears. The SSG/LRR-SAS and hybrid SSG/LRR models both predict a pair of strong vortices in the shape of the ears. This flow structure is similar with that obtained by SSG/LRR-URANS. The SSG/LRR-IDDES predict a vortex structure formed from the shoulders of the cylinder, which agrees better with the experimental oil-surface-flow visualizations (bold arrow in Figure 11) than that of SSG/LRR-URANS. It is consistent with the trend of the instantaneous solutions. Similar differences also exist in SST-IDDES and SST-URANS, which indicates that the type of scale correction has a major influence on the wake region.

Figure 10.

Comparison of time mean near surface flow topology.

Figure 11.

The oil-surface-flow visualization from the experiment. ⁵⁷

A comparison of the numerical and experimental mean-surface pressure coefficient is shown in Figure 12. All RSM-based solutions exhibit the same pressure drop as the flow accelerates over the front portion of the turret dome. However in the separated flow region, all solutions are slightly deviated from the experimental data. Morgan and Visbal ⁵⁸ implemented the gird convergence analysis on this case using hybrid RANS/ILES method and indicted that the simulation accuracy can be improved by reducing the grid size in this region. As far as this grid is concerned, the RSM-based solutions are better than SST-based solutions.

Figure 12.

Time-mean surface-pressure coefficient comparison for turret dome.

A comparison of mean u-velocity profiles at four locations in the wake flow are shown in Figure 13. The location of the profiles is given in Figure 14. Velocity profiles 2 and 4 are located at the $z=-0.5$ plane near where the horseshoe vortex and wake vortices interact. All methods show a good agreement with experimental data, mainly in the outer portion of the boundary layer. Velocity profiles 3 lies on the symmetry plane, where the SST-IDDES gives the best u-velocity profile. Velocity profiles 5 is located at the downstream of the velocity profiles 3. There is no essential difference in all methods.

Figure 13.

Comparison of time-mean u-velocity profile.

Figure 14.

u-velocity profile locations.

Conclusions

Based on the high-moment turbulence model SSG/LRR- $\omega$ , three hybrid RANS/LES methods (SAS, IDDES and hybrid filtering models) are carried out. The IDDES and hybrid filtering models have the same aims to establish the correlation between turbulent dissipation and the local grid scale in the separation regions. Both approaches behave like a subgrid-stress transport model when the turbulent length scale is greater than the grid length scale, and the dissipation equation is decoupled from the turbulence-stress equations in this region. The SAS is achieved by the introduction of the second derivative of the velocity field into the turbulence scale equation. Compared with URANS, the difference is the smallest.

A comparative assessment of SSG/LRR-SAS, SSG/LRR-IDDES and hybrid SSG/LRR models is made with the high-order WCNS-E8T7 scheme. At the same time, the traditional SST-IDDES is also used as a reference. The relevant cases include NACA0012 airfoil stalled flows and flow over a turret. The former is a benchmark case for hybrid RANS/LES models, while the latter is a small challenge to the turbulence models and discretization methods. In above cases, all of the three new methods show a good ability to simulate unsteady turbulence, while the original SSG/LRR-URANS method is insufficient. Among the hybrid RANS/LES methods, the calculation results obtained by SSG/LRR-IDDES are generally the best. The hybrid SSG/LRR model shows more physical viscosity than SSG/LRR-IDDES. It is due to the fact that Deardorff’s pressure-strain correlation gives an insufficient dissipation in the separated region. The SAS modification can also cause the physical viscosity to be reduced in the separation regions. But unlike the IDDES or hybrid framework, it is based on the use of a von Karman length scale, which always relies on the accuracy of the calculated velocity gradient.

This paper preliminarily shows the potential of hybrid RANS/LES methods base on a Reynolds-stress model by comparing with the SST-IDDES. Like other hybrid methods based on the EVM, there may still be some defects in above methods. Therefore, future work needs to consider more test cases.