A new formulation of f_k for the PANS model

ABSTRACT

The partially averaged Navier–Stokes (PANS) model, proposed in Girimaji [Partially-averaged Navier-Stokes model for turbulence: a Reynolds-averaged Navier-Stokes to direct numerical simulation bridging method. ASME J Appl Mech. 2006;73(3):413–421], can be used to simulate turbulent flows either as a RANS, LES or DNS. The PANS model includes $f_k$ which denotes the ratio of modelled to total kinetic energy. In RANS, $f_k=1$, and in DNS it goes to zero. In the present study we propose a new formulation for $f_k$ based on the H-equivalence introduced by Friess et al. [Toward an equivalence criterion for hybrid RANS/LES methods. Int J Heat Fluid Flow. 2015;122:233–246]. In this formulation the expression of $f_k$ is derived to mimic DES. This new formulation behaves very much like classic DES, even though the two formulations use different mechanisms to separate modelled and resolved scales. They show very similar performance in separated flows as well as in attached boundary layers. Moreover, the new formulation exhibits similar robustness features as DES.

KEYWORDS:

• LES
• PANS
• DES
• channel flow
• hill flow
• hump flow

1. Introduction

The PANS model was proposed by Girimaji [1] and the PITM (Partially Integrated Transport Model) was proposed by Schiestel and Dejoan [2], Chaouat and Schiestel [3]. The critical parameter in both models is $f_k$ (it is sometimes called r in PITM). It is defined as the ratio of the modelled to the total turbulent kinetic energy: (1) $f_k=\frac{k}{k_{\mathrm{t}\mathrm{o}\mathrm{t}}}.$(1) In DNS (Direct Numerical Solution), it tends to zero and in RANS (Reynolds-Averaged Navier–Stokes) it should be one. In LES (Large Eddy Simulation), the parameter takes values between zero and one. It is natural to link $f_k$ to the mesh resolution and in that aim several proposals have been made on how to compute $f_k$. By seeking the smallest resolved length scale for a given $f_k$ through dimensional analysis, Girimaji and Abdul-Hamid [4] proposed one way to compute $f_k$: (2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) using $\mathrm{\Delta }={\mathrm{\Delta }}_{\mathrm{m}\mathrm{i}\mathrm{n}}$, the smallest grid cell size. Basara et al. also used $f_k$ prescribed from Equation (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ), however taking the geometric average $\mathrm{\Delta }=(\mathrm{\Delta }V{)}^{1/3}$. Kenjeres and Hanjalic [5] have made a slightly different proposal which reads (3) $f_k=\frac{\mathrm{\Delta }}{L_t}$(3) Another way is to compute $f_k$ from its very definition (1(1) $f_k=\frac{k}{k_{\mathrm{t}\mathrm{o}\mathrm{t}}}.$(1) ). $k_{\mathrm{t}\mathrm{o}\mathrm{t}}$ is then computed using the running average. The expression (1(1) $f_k=\frac{k}{k_{\mathrm{t}\mathrm{o}\mathrm{t}}}.$(1) ) shall hereafter be refered to as $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$, the observed $f_k$. More recently, Foroutan and Yavuzkurt [6] derived an expression from the Kolmogorov energy spectrum which reads (4) $f_k=1-{\left[\frac{(\mathrm{\Lambda }/\mathrm{\Delta }{)}^{2/3}}{0.23+(\mathrm{\Lambda }/\mathrm{\Delta }{)}^{2/3}}\right]}^{9/2}.$(4) In [7], the expression in Equation (3(3) $f_k=\frac{\mathrm{\Delta }}{L_t}$(3) ) was found to give far too small $f_k$. The form in Equation (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ) and (1(1) $f_k=\frac{k}{k_{\mathrm{t}\mathrm{o}\mathrm{t}}}.$(1) ) were evaluated but it was found that a constant $f_k=0.4$ in the LES region is superior.

The present paper is based on the work in [8] where they derived a relation between DES and PITM. They showed that the DES model could be formulated using $f_k$. They call this model an equivalent DES model. The relation between DES and PANS is used in the present work, but it is used the other way around: a new form of $f_k$ is derived based on the DES model. In [9] the form of $f_k$ (new PANS, detailed in Section 2.1) was compared to DES. and the two models were found to give more or less identical results, particularly regarding first-order moments: velocity, skin friction coefficient $C_f$. In the present work the new PANS will be compared to the expression in Equation (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ) (the old PANS using Δ=(ΔV)^1/3) which in the literature is the most common way to compute $f_k$. 

The paper is organised as follows. First, the form of $f_k$ is derived and analysed analytically. In the following section we present the numerical method. Then the new and old PANS are compared in three different flows (channel flow, hump flow and hill flow). Some conclusions are drawn in the final section.

2. The PANS model

The low-Reynolds number Partially-Averaged Navier–Stokes (LRN PANS, see [10]) turbulence model reads (5) $\begin{array}{rl}\frac{\mathrm{d}k}{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{ku}}\right)\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\right]+P_k-\epsilon \\ \frac{\mathrm{d}\epsilon }{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon u}}\right)\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{x}_j}\right]+C_{\epsilon 1}{P}_k\frac{\epsilon }{k}-C_{\epsilon 2}^{\ast }\frac{{\epsilon }^2}{k}\\ {\nu }_t& =C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon },P_k=2{\nu }_t{\overline{s}}_{ij}{\overline{s}}_{ij},{\overline{s}}_{ij}=\frac{1}{2}\left(\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{\overline{v}}_j}{\mathrm{\partial }{x}_i}\right)\\ C_{\epsilon 2}^{\ast }& =C_{\epsilon 1}+\frac{f_k}{f_{\epsilon }}(C_{\epsilon 2}{f}_2-C_{\epsilon 1}),{\sigma }_{ku}\equiv {\sigma }_k\frac{f_k^2}{f_{\epsilon }},{\sigma }_{\epsilon u}\equiv {\sigma }_{\epsilon }\frac{f_k^2}{f_{\epsilon }}\\ {\sigma }_k& =1.4,{\sigma }_{\epsilon }=1.4,C_{\epsilon 1}=1.5,C_{\epsilon 2}=1.9,C_{\mu }=0.09,\end{array}$(5) where $\mathrm{d}/\mathrm{d}t=\mathrm{\partial }/\mathrm{\partial }t+{\overline{v}}_j\mathrm{\partial }/\mathrm{\partial }{x}_j$ denotes the material derivative. The damping functions are given by (6) $\begin{array}{rl}{f}_2& ={\left[1-\exp (-\frac{y^{\ast }}{3.1})\right]}^2\left\{1-0.3\mathrm{e}\mathrm{x}\mathrm{p}[-(\frac{R_t}{6.5}{)}^2]\right\}\\ f_{\mu }& ={\left[1-\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{y^{\ast }}{14})\right]}^2\left\{1+\frac{5}{R_t^{3/4}}\mathrm{e}\mathrm{x}\mathrm{p}[-(\frac{R_t}{200}{)}^2]\right\}\\ R_t& =\frac{k^2}{\nu \epsilon },y^{\ast }=\frac{U_{\epsilon }y}{\nu },U_{\epsilon }=(\epsilon \nu {)}^{1/4}.\end{array}$(6) The functions $f_k$ and $f_{\epsilon }$ denote the ratio of modelled to total kinetic energy and modelled to total dissipation, respectively. For flows at high Reynolds numbers (as in the present work), the dissipation is modelled which means that $f_{\epsilon }=1$. In the PITM model, ${\sigma }_{ku}\equiv {\sigma }_k$ and ${\sigma }_{\epsilon u}\equiv {\sigma }_{\epsilon }$.

2.1. $f_k$ derived from the equivalence criterion

In [8] a relation between $f_k$ and the grid step is derived, through the establishment of a statistical equivalence between DES and PITM. To that aim, they performed perturbation analyses about the equilibrium states, representing small variation of the energy partition. They did the analysis with and without considering inhomogeneity. That derivation is summarised here in a homogeneous framework, as a first step. Let us first consider the PANS/PITM equations. For equilibrium turbulence $d\tau /dt=0$ where $\tau =k/\epsilon$, Equation (5(5) $\begin{array}{rl}\frac{\mathrm{d}k}{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{ku}}\right)\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\right]+P_k-\epsilon \\ \frac{\mathrm{d}\epsilon }{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon u}}\right)\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{x}_j}\right]+C_{\epsilon 1}{P}_k\frac{\epsilon }{k}-C_{\epsilon 2}^{\ast }\frac{{\epsilon }^2}{k}\\ {\nu }_t& =C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon },P_k=2{\nu }_t{\overline{s}}_{ij}{\overline{s}}_{ij},{\overline{s}}_{ij}=\frac{1}{2}\left(\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{\overline{v}}_j}{\mathrm{\partial }{x}_i}\right)\\ C_{\epsilon 2}^{\ast }& =C_{\epsilon 1}+\frac{f_k}{f_{\epsilon }}(C_{\epsilon 2}{f}_2-C_{\epsilon 1}),{\sigma }_{ku}\equiv {\sigma }_k\frac{f_k^2}{f_{\epsilon }},{\sigma }_{\epsilon u}\equiv {\sigma }_{\epsilon }\frac{f_k^2}{f_{\epsilon }}\\ {\sigma }_k& =1.4,{\sigma }_{\epsilon }=1.4,C_{\epsilon 1}=1.5,C_{\epsilon 2}=1.9,C_{\mu }=0.09,\end{array}$(5) ) gives (7) $\begin{array}{rl}\frac{\mathrm{d}\tau }{\mathrm{d}t}& =\frac{1}{\epsilon }\frac{\mathrm{d}k}{\mathrm{d}t}-\frac{k}{{\epsilon }^2}\frac{\mathrm{d}\epsilon }{\mathrm{d}t}=\frac{1}{\epsilon }\left(P^k+D^k-\epsilon \right)\\ & -\frac{k}{{\epsilon }^2}\left(C_{\epsilon 1}\frac{\epsilon }{k}{P}^k+D^{\epsilon }-C_{\epsilon 2}^{\ast }\frac{{\epsilon }^2}{k}\right)=0,\end{array}$(7) where $D^k$ and $D^{\epsilon }$ denote the diffusion term for k and ϵ, respectively. For local homogeneous turbulence (i.e. $D^k=D^{\epsilon }=0$), it can be written (8) $\begin{array}{rl}& \gamma (C_{\epsilon 1}-1)Sk=(C_{\epsilon 2}^{\ast }-1)\epsilon \\ & \gamma =\frac{P^k}{Sk},S=(2{\overline{s}}_{ij}{\overline{s}}_{ij}{)}^{1/2}.\end{array}$(8) The quantities that are affected by the partition between modelled and resolved turbulence (i.e. $f_k$) in Equation (8(8) $\begin{array}{rl}& \gamma (C_{\epsilon 1}-1)Sk=(C_{\epsilon 2}^{\ast }-1)\epsilon \\ & \gamma =\frac{P^k}{Sk},S=(2{\overline{s}}_{ij}{\overline{s}}_{ij}{)}^{1/2}.\end{array}$(8) ) are γ, S, k and $C_{\epsilon 2}^{\ast }$.¹ Differentiation of Equation (8(8) $\begin{array}{rl}& \gamma (C_{\epsilon 1}-1)Sk=(C_{\epsilon 2}^{\ast }-1)\epsilon \\ & \gamma =\frac{P^k}{Sk},S=(2{\overline{s}}_{ij}{\overline{s}}_{ij}{)}^{1/2}.\end{array}$(8) ), by considering infinitesimal perturbations $\delta \gamma$, $\delta S$, $\delta k$ and $\delta C_{\epsilon 2}^{\ast }$ of the variables, yields: (9) $\delta \gamma Sk+\delta S\gamma k+\delta k\gamma S=\frac{\delta C_{\epsilon 2}^{\ast }\epsilon }{C_{\epsilon 1}-1}$(9) so that (10) $\frac{\delta \gamma }{\gamma }+\frac{\delta S}{S}+\frac{\delta k}{k}=\frac{\delta C_{\epsilon 2}^{\ast }\epsilon }{(C_{\epsilon 1}-1)\gamma Sk}=\frac{\delta C_{\epsilon 2}^{\ast }}{C_{\epsilon 2}^{\ast }-1}.$(10) Equation (10(10) $\frac{\delta \gamma }{\gamma }+\frac{\delta S}{S}+\frac{\delta k}{k}=\frac{\delta C_{\epsilon 2}^{\ast }\epsilon }{(C_{\epsilon 1}-1)\gamma Sk}=\frac{\delta C_{\epsilon 2}^{\ast }}{C_{\epsilon 2}^{\ast }-1}.$(10) ) was derived for the PANS/PITM equations. Now we repeat the derivation for the DES equations. The differences between DES and PITM/PANS are that in DES (i) $C_{\epsilon 2}^{\ast }=C_{\epsilon 2}$ is constant and (ii) the dissipation term in the equation for modelled energy k is replaced with $\psi \epsilon$, i.e. (11) $\begin{array}{rl}\frac{\mathrm{d}k}{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\right]+P_k-\psi \epsilon \end{array}$(11) (12) $\begin{array}{rl}\psi & =max\left(1,\frac{k^{3/2}/\epsilon }{C_{DES}{\mathrm{\Delta }}_{max}}\right),{\mathrm{\Delta }}_{max}=max(\mathrm{\Delta }{x}_1,\mathrm{\Delta }{x}_2,\mathrm{\Delta }{x}_3)\\ \frac{\mathrm{d}\epsilon }{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{x}_j}\right]+C_{\epsilon 1}{P}_k\frac{\epsilon }{k}-C_{\epsilon 2}\frac{{\epsilon }^2}{k}.\end{array}$(12) Assuming $\mathrm{d}\tau /\mathrm{d}t=0$ and local homogeneous turbulence gives (13) $\gamma (C_{\epsilon 1}-1)Sk=(C_{\epsilon 2}-\psi )\epsilon .$(13) We differentiate so that (14) $\frac{\delta \gamma }{\gamma }+\frac{\delta S}{S}+\frac{\delta k}{k}=-\frac{\mathrm{d}\psi \epsilon }{(C_{\epsilon 1}-1)Sk\gamma }=-\frac{\mathrm{d}\psi }{C_{\epsilon 2}-\psi }.$(14) Equations (9(9) $\delta \gamma Sk+\delta S\gamma k+\delta k\gamma S=\frac{\delta C_{\epsilon 2}^{\ast }\epsilon }{C_{\epsilon 1}-1}$(9) ) and (14(14) $\frac{\delta \gamma }{\gamma }+\frac{\delta S}{S}+\frac{\delta k}{k}=-\frac{\mathrm{d}\psi \epsilon }{(C_{\epsilon 1}-1)Sk\gamma }=-\frac{\mathrm{d}\psi }{C_{\epsilon 2}-\psi }.$(14) ) describe how $C_{\epsilon 2}^{\ast }$ and ψ depend on variations in γ, S and k. The parameters $C_{\epsilon 2}^{\ast }$ and ψ vary from $C_{\epsilon 2}$ and 1 (RANS values), respectively, to $C_{\epsilon 2}^{\ast }$ and $\psi \left(\mathrm{\Delta }\right)$ (LES values). Combining Equations (9(9) $\delta \gamma Sk+\delta S\gamma k+\delta k\gamma S=\frac{\delta C_{\epsilon 2}^{\ast }\epsilon }{C_{\epsilon 1}-1}$(9) ) and (14(14) $\frac{\delta \gamma }{\gamma }+\frac{\delta S}{S}+\frac{\delta k}{k}=-\frac{\mathrm{d}\psi \epsilon }{(C_{\epsilon 1}-1)Sk\gamma }=-\frac{\mathrm{d}\psi }{C_{\epsilon 2}-\psi }.$(14) ) and integrating from RANS to LES conditions ($C_{\epsilon 2}^{\ast }$ and ψ) (15) ${\int }_{C_{\epsilon 2}}^{C_{\epsilon 2}^{\ast }}\frac{\mathrm{d}{C}_{\epsilon 2}^{\ast }}{C_{\epsilon 2}^{\ast }-1}={\int }_1^{\psi }-\frac{\mathrm{d}\psi }{C_{\epsilon 2}-\psi }\Rightarrow \ln \left(\frac{C_{\epsilon 2}^{\ast }-1}{C_{\epsilon 2}-1}\right)=\ln \left(\frac{C_{\epsilon 2}-\psi }{C_{\epsilon 2}-1}\right).$(15) By using the expression for $C_{\epsilon 2}^{\ast }$ in Equation (5(5) $\begin{array}{rl}\frac{\mathrm{d}k}{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{ku}}\right)\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_j}\right]+P_k-\epsilon \\ \frac{\mathrm{d}\epsilon }{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon u}}\right)\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{x}_j}\right]+C_{\epsilon 1}{P}_k\frac{\epsilon }{k}-C_{\epsilon 2}^{\ast }\frac{{\epsilon }^2}{k}\\ {\nu }_t& =C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon },P_k=2{\nu }_t{\overline{s}}_{ij}{\overline{s}}_{ij},{\overline{s}}_{ij}=\frac{1}{2}\left(\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{\overline{v}}_j}{\mathrm{\partial }{x}_i}\right)\\ C_{\epsilon 2}^{\ast }& =C_{\epsilon 1}+\frac{f_k}{f_{\epsilon }}(C_{\epsilon 2}{f}_2-C_{\epsilon 1}),{\sigma }_{ku}\equiv {\sigma }_k\frac{f_k^2}{f_{\epsilon }},{\sigma }_{\epsilon u}\equiv {\sigma }_{\epsilon }\frac{f_k^2}{f_{\epsilon }}\\ {\sigma }_k& =1.4,{\sigma }_{\epsilon }=1.4,C_{\epsilon 1}=1.5,C_{\epsilon 2}=1.9,C_{\mu }=0.09,\end{array}$(5) ) (with $f_2=f_{\epsilon }=1$), and ensuring that $0<f_k\le 1$ we finally get (16) $f_k=max\left[0,min\left(1,1-\frac{\psi -1}{C_{\epsilon 2}-C_{\epsilon 1}}\right)\right],$(16) where ψ is given by Equation (12(12) $\begin{array}{rl}\psi & =max\left(1,\frac{k^{3/2}/\epsilon }{C_{DES}{\mathrm{\Delta }}_{max}}\right),{\mathrm{\Delta }}_{max}=max(\mathrm{\Delta }{x}_1,\mathrm{\Delta }{x}_2,\mathrm{\Delta }{x}_3)\\ \frac{\mathrm{d}\epsilon }{\mathrm{d}t}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{x}_j}\right]+C_{\epsilon 1}{P}_k\frac{\epsilon }{k}-C_{\epsilon 2}\frac{{\epsilon }^2}{k}.\end{array}$(12) ). This may be compared with the old formulations in Equations (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ), (3(3) $f_k=\frac{\mathrm{\Delta }}{L_t}$(3) ) and (4(4) $f_k=1-{\left[\frac{(\mathrm{\Lambda }/\mathrm{\Delta }{)}^{2/3}}{0.23+(\mathrm{\Lambda }/\mathrm{\Delta }{)}^{2/3}}\right]}^{9/2}.$(4) ). However, the present study focuses on a comparison between Equations (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ) and (16(16) $f_k=max\left[0,min\left(1,1-\frac{\psi -1}{C_{\epsilon 2}-C_{\epsilon 1}}\right)\right],$(16) ).

2.2. Self-adaptivity properties of the new formulation of $f_k$ for PANS

In this section, we will distinguish $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ (see Equation (1(1) $f_k=\frac{k}{k_{\mathrm{t}\mathrm{o}\mathrm{t}}}.$(1) )), the observed energy ratio, from $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$, the targeted (or prescribed) energy ratio. In the turbulent closure equations (i.e.in PANS, computing $C_{\epsilon 2}^{\ast }$), $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ is used, but nothing ever guarantees that this level of modelled energy will be exactly reached, see e.g. [11]. This can be explained by the fact that $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ is a rough estimate, obtained under assumptions that are not always fully valid. Meanwhile, $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ can be obtained from postprocessing, using its definition in Equation (1(1) $f_k=\frac{k}{k_{\mathrm{t}\mathrm{o}\mathrm{t}}}.$(1) ).

If the targeted energy ratio, $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$, is computed following Equation (16(16) $f_k=max\left[0,min\left(1,1-\frac{\psi -1}{C_{\epsilon 2}-C_{\epsilon 1}}\right)\right],$(16) ), and if we use the definition (which is rigorous in average) (17) $k\simeq f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}{k}_{\mathrm{t}\mathrm{o}\mathrm{t}},\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{k}_{\mathrm{t}\mathrm{o}\mathrm{t}}=⟨k⟩+⟨k_{\mathrm{r}\mathrm{e}\mathrm{s}}⟩$(17) and the assumption that the whole dissipation is modelled, i.e. $f_{\epsilon }=1$ (at high Reynolds number): (18) $\epsilon \approx {\epsilon }_{\mathrm{t}\mathrm{o}\mathrm{t}},$(18) then Equation (16(16) $f_k=max\left[0,min\left(1,1-\frac{\psi -1}{C_{\epsilon 2}-C_{\epsilon 1}}\right)\right],$(16) ) can be rewritten as: (19) $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}=1-\frac{1}{C_{\epsilon 2}-C_{\epsilon 1}}\left[max\left(1,\frac{f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}^{3/2}}{C_{\mathrm{D}\mathrm{E}\mathrm{S}}{\mathrm{\Delta }}_{max}}\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{{\epsilon }_{\mathrm{t}\mathrm{o}\mathrm{t}}}\right)-1\right].$(19) (note that, for sake of clarity, the form of $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ is here allowed to go negative). This implies that $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ is implicitly linked to $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$, in the following way: if $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ is lower (resp. higher) than a certain threshold value, $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ will increase (resp. decrease), leading dynamically to an increase (resp. decrease) of k and thus $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$.

Of course, this rough reasoning assumes that the resolved part of k ‘responds’ to those changes quickly enough to leave $k_{\mathrm{t}\mathrm{o}\mathrm{t}}$ almost unaffected.

Thus, one can conclude that defining $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ according to Equation (16(16) $f_k=max\left[0,min\left(1,1-\frac{\psi -1}{C_{\epsilon 2}-C_{\epsilon 1}}\right)\right],$(16) ) drives $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ towards a certain threshold (which may be more or less close to $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$). Such a feature is not intrinsic to approaches such as PITM or PANS, if some spectral law or Equations (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ), (3(3) $f_k=\frac{\mathrm{\Delta }}{L_t}$(3) ) or (4(4) $f_k=1-{\left[\frac{(\mathrm{\Lambda }/\mathrm{\Delta }{)}^{2/3}}{0.23+(\mathrm{\Lambda }/\mathrm{\Delta }{)}^{2/3}}\right]}^{9/2}.$(4) ) are used to determine $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$. It is actually intrinsic to Detached Eddy Simulation (DES), from which the present new formulation of $f_k$ inherits this interesting property. This strategy may be interpreted as passive control, since it does not require the explicit computation of any extra quantity, such as $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$.

3. Numerical solver

An incompressible, finite volume code is used [12]. The convective terms in the momentum equations are discretised using central differencing except for the hump flow where 95% is taken from CDS and 5% from the second-order upwind MUSCL scheme. Hybrid central/upwind is used for the k and ϵ equations. The Crank–Nicolson scheme is used for time discretisation of all equations. The numerical procedure is based on an implicit, fractional step technique with a multigrid pressure Poisson solver [13] and a non-staggered grid arrangement.

The filtered momentum equations with an added turbulent viscosity read (20) $\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\overline{v}}_j{\overline{v}}_i}{\mathrm{\partial }{x}_j}=\beta {\delta }_{1i}-\frac{1}{\rho }\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\left(\nu +{\nu }_t\right)\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }{x}_j}\right),$(20) where the first term on the right side is the driving pressure gradient in the streamwise direction, which is used in the fully-developed channel flow simulations and for the hill flow.

4. Results

The new formulation of $f_k$ is evaluated and compared with the old formulation given in Equation (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ). Herafter, they will be referred to as ‘new PANS’ and ‘old PANS’. The comparison is performed in three test cases: fully developed channel flow, the hump flow (see Figure 1), and the hill flow (see Figure 2).

Figure 1. The geometry of the hump.

Figure 2. The geometry of the hill.

Ensemble-averaged quantities are plotted hereafter, i.e. they are averaged in time and over statistically homogeneous directions (which differ with flow cases).

4.1. Channel flow

The Reynolds number is defined as $R{e}_{\tau }=u_{\tau }\delta /\nu =5200$ where δ denotes half channel height and $u_{\tau }$ is the friction velocity. The streamwise, wall-normal and spanwise directions are denoted by x, y and z, respectively. The size of the domain is $x_{max}=3.2$, $y_{max}=2$ and $z_{max}=1.6$. The mesh has $32\times 96\times 32$ cells in the x−y−z directions. Periodic boundary conditions are used in the x and z directions. Therefore, these two directions are considered statistically homogeneous. A precursor DES computation is used as initial condition. The driving pressure gradient (first term on the right-hand side in Equation (20(20) $\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\overline{v}}_j{\overline{v}}_i}{\mathrm{\partial }{x}_j}=\beta {\delta }_{1i}-\frac{1}{\rho }\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\left(\nu +{\nu }_t\right)\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }{x}_j}\right),$(20) ) is used with $\beta =1$. A lower limit of 0.05 is used when computing $f_k$ from Equation (16(16) $f_k=max\left[0,min\left(1,1-\frac{\psi -1}{C_{\epsilon 2}-C_{\epsilon 1}}\right)\right],$(16) ).

Figure 3 presents the mean velocity and the turbulent kinetic energy profiles.As can be seen, the new PANS exhibits a small bump near $y^{+}=1000$; otherwise it agrees well with the DNS profile. The predicted turbulent kinetic energies with both old and new PANS agree well with DNS for $y^{+}⪆500$. However, there is a big difference in how much turbulence is resolved and how much is modelled. For $y^{+}⪆500$, all turbulence is resolved by the new PANS model. In the old PANS model, all turbulence is modelled. The reason is found in Figure 4(b)): the old PANS predicts $f_k\simeq 1$ in the entire region which means that the model is working in RANS mode. This is also seen in Figure 4(a)) in which it can be seen that the turbulent viscosity predicted by the old PANS is one-order magnitude larger than the new PANS model far from the wall. Furthermore, Figure 4(b)) shows that the new PANS exhibits a real plateau of $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}=1$ in the near-wall region, leading to a sharp switch from RANS to LES (at $y^{+}\simeq 200$, see insets). On the contrary, the old PANS does not switch to LES at all but stays in RANS mode in the entire channel. For both old and new PANS, the observed $f_k$ is significantly different from the prescribed $f_k$, see Figure 4(b)). However, for the new PANS, $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ is smaller than $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ everywhere, but the profiles have a quite similar shape.

Figure 3. Channel flow, $R{e}_{\tau }=5200$. Blue lines: new PANS; red lines: old PANS, see Equation (2(2) $f_k=C_{\mu }^{-1/2}{\left(\frac{\mathrm{\Delta }}{L_t}\right)}^{2/3},L_t=\frac{k_{\mathrm{t}\mathrm{o}\mathrm{t}}^{3/2}}{\epsilon }$(2) ) $R{e}_{\tau }=5200$. Markers: DNS [14].

4.2. Hump flow

The Reynolds number of the hump flow is $R{e}_c=936,000$, based on the hump length, c=1, and the inlet mean velocity at the centreline, $U_{\mathrm{i}\mathrm{n},c}$. In the present simulations, the value of ρ, c and $U_{\mathrm{i}\mathrm{n},c}$ have been set to unity. The configuration is given in Figure 1. Experiments were conducted by Greenblatt et al. [15, 16]. The maximum height of the hump, h, and the channel height, H, are given by $h/c=0.128$ and $H/c=0.91$, respectively. The mesh has $304\times 108\times 32$ cells and is taken from the NASA workshop.² The spanwise extent is set to $Z_{max}/c=0.2$. The inlet is located at $x/c=0.5$ and the outlet at $x/c=4.0$.

Figure 4. Channel flow. Viscosity and $f_k$. $R{e}_{\tau }=5200$.

A periodic boundary condition is applied in the spanwise direction z. Therefore, this direction is considered statistically homogeneous. The inlet conditions (U, V, k and ϵ) are taken from a 2D RANS simulation using the AKN $k-\epsilon$ turbulence model [17] coupled to the EARSM model [18]. Synthetic isotropic fluctuations are superimposed on the 2D RANS velocity field. The synthetic fluctuations are scaled with the RANS shear stress profile. To reduce the inlet k, prescribed from 2D RANS, a commutation term $\mathrm{\partial }{f}_k/\mathrm{\partial }x$ is used. For more detail on inlet synthetic fluctuations and the commutation term, see [7]. A lower limit of 0.2 is used when computing $f_k$ from Equation (16(16) $f_k=max\left[0,min\left(1,1-\frac{\psi -1}{C_{\epsilon 2}-C_{\epsilon 1}}\right)\right],$(16) ).

The simulations are initialised as follows [19]: first the 2D RANS equations are solved. Anisotropic synthetic fluctuations, $({\mathcal{V}}_i^{\prime }{)}_m$, are then superimposed to the 2D RANS field which gives the initial LES velocity field. In order to compute $({\mathcal{V}}_i^{\prime }{)}_m$, synthetic fluctuations, $v_{i,\mathrm{s}\mathrm{y}\mathrm{n}\mathrm{t}}^{\prime }$, are computed plane-by-plane (y−z) in the same way as prescribing inlet boundary conditions. The synthetic fluctuations in the y−z planes are coupled with an asymmetric space filter (21) $({\mathcal{V}}_i^{\prime }{)}_m=a({\mathcal{V}}_i^{\prime }{)}_{m-1}+b(v_{\mathrm{s}\mathrm{y}\mathrm{n}\mathrm{t},i}^{\prime }{)}_m,$(21) where m denotes the index of the $x_1$ location and $a=\exp (-\mathrm{\Delta }{x}_1/L_{\mathrm{i}\mathrm{n}\mathrm{t}})$ and $\mathrm{\Delta }{x}_1$ and $L_{\mathrm{i}\mathrm{n}\mathrm{t}}$ denote the grid size and the integral length scale, respectively ($L_{\mathrm{i}\mathrm{n}\mathrm{t}}=0.2$ ).

Figures 5 and 6 compare predictions with experiments and, as can be seen, the new PANS is in a very good agreement with experiments, while the old PANS exhibits a significant discrepancy at $x/c>1$. The predicted skinfrictions show a small bump near the inlet, and the reason is (at least partly) that a different RANS turbulence model (EARSM) was used in the 2D RANS simulations than the underlying RANS model in the PANS simulations. The backflow region is well predicted by both PANS models, but globally, the new PANS predicts the streamwise velocity far better than the old PANS.

Figure 5. Hump flow. Pressure coefficient and skinfriction.

: new PANS;

: old PANS; markers: Experiments [15, 16].

Figure 6. Hump flow. Streamwise velocities.

: new PANS;

: old PANS; markers: Experiments [15, 16].

Figure 7 shows the turbulent viscosity for both models (note that the turbulent viscosities given by the old PANS have been multiplied by 0.1, for the sake of presentation). As in the channel flow, the turbulent viscosities predicted by the old PANS are an order of magnitude larger than those predicted by the new PANS. This means that the solution predicted by the old PANS is closer to RANS than to LES.

Figure 7. Hump flow. Turbulent viscosity.

: new PANS;

: $0.1{\nu }_t/\nu$, old PANS.

Figure 8 presents $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ and $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ for both old and new PANS at three locations, $x/c=0.65,1.1$ and 1.3. The old PANS gives a much larger $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ (more RANS) than the new PANS as was noted above in Figure 7. The difference between the two model is largest in the attached flow region. Furthermore, for both approaches, $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ is significantly different from $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$, $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ being really low, although the flow is attached to the wall at all three locations. This may be due to the fact that the overall instability of the flow generates a big amount of fluctuations, thus increasing $k_{\mathrm{r}\mathrm{e}\mathrm{s}}$, leading to a reduction of $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$.

Figure 8. Hump flow.

:$f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$: new PANS;

: $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$: old PANS;

: $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$: new PANS;

: $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$: old PANS.

4.3. Hill flow

The domain is shown in Figure 2. The size of the domain is $9H\times 3.035H\times 4.5H$ in the streamwise (x), wall-normal (y) and span-wise direction (z), respectively. The grid has $160\times 80\times 32$ cells in the x, y and z direction. Periodic boundary conditions are used in the x and z directions. The z direction is considered statistically homogeneous. Slip conditions are prescribed at the upper wall. The Reynolds number is Re=10,600 based on the hill height and the bulk velocity at the top of the hill. An initial velocity field is prescribed from a 2D RANS solution with the correct bulk Reynolds number. Furthermore, the same technique for synthetic turbulence as for the hump flow, is used to add initial fluctuations. The bulk velocity is then kept constant by adjusting β in Equation (20(20) $\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\overline{v}}_j{\overline{v}}_i}{\mathrm{\partial }{x}_j}=\beta {\delta }_{1i}-\frac{1}{\rho }\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\left(\nu +{\nu }_t\right)\frac{\mathrm{\partial }{\overline{v}}_i}{\mathrm{\partial }{x}_j}\right),$(20) ) at each time step by ensuring that the sum of the forces at the wall (wall shear stress and pressure on the lower wall) balances the driving pressure gradient [20, 21, Section 4.5]. A lower limit of 0.2 is used when computing f_k from Equation (16).

Figure 9 compares the velocity field with LES of [22] and the agreement with the new PANS model is excellent.However, the old PANS model performs very poorly. The reason for the poor performance of the old PANS is seen in Figure 10 (note that the viscosities of the old PANS have been multiplied by 0.1); the turbulent viscosity is one order of magnitude larger with the old PANS compared to the new PANS.

Figure 9. Hill flow. Velocities.

: new PANS;

: old PANS; markers: LES [22].

Figure 11 compares profiles of $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ and $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$ for both PANS approaches. It is interesting to notice that at the top of the hill ($x/H=0.05$), the two $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$ are close to each other, while their observed counterparts are significantly different: the new PANS is very well resolved, while the old one is closer to the RANS level. Further downstream, this global tendency is conserved, except that old and new PANS stray from each other, even after reattachment ($x/H=6$).

Figure 10. Hill flow. Turbulent viscosities.

: new PANS;

: $0.1{\nu }_t/\nu$, old PANS.

Figure 11. Hill flow.

:$f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$, new PANS;

: $f_{k,\mathrm{t}\mathrm{a}\mathrm{r}}$: old PANS;

: $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$: new PANS;

: $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$: old PANS.

5. Concluding remarks

A new formulation for prescribing $f_k$ has been presented. It has been found to perform much better than the standard form of $f_k$ of [4, 23]. It should, however, be mentioned that the standard form of $f_k$ works much better when used in a four-equation turbulence model [23] than in a two-equation model as in the present work, presumably because the four-equation closure contains more accurate near-wall physics.

The new formulation presented here, behaves very much like ‘classic DES’, as one can see in [9], even though they use different mechanisms to separate modelled and resolved scales. In particular, they show very similar performance in separated flows, as well as in attached boundary layers.

Another interesting feature of the new PANS, is the robustness inherited from DES. As explained in Section 2.2, it behaves like a passive control device, in the sense that it does not require the explicit computation of the observed energy ratio $f_{k,\mathrm{o}\mathrm{b}\mathrm{s}}$.

So what is the advantage of using PANS instead of DES? One advantage is that PANS has a much stronger theoretical foundation than DES. PANS is rigorously derived whereas DES is an ad-hoc (but very successful) modification of a RANS model. Another advantage of PANS is that the modified partition between modelled and resolved turbulence due to non-uniform grids can be accounted for by adding a term in the k and momentum equations based on the gradient of $f_k$ [7, 24]. Future work will focus on a thorough theoretical derivation of the relationship between $f_k$ and the grid step, by accounting for inhomogeneity. Another test will consist of prescribing $f_k$ with some more elaborate forms of DES, for instance Improved Delayed Detached Eddy Simulation (IDDES, see e.g. [25, 26]).

Disclosure statement

No potential conflict of interest was reported by the authors.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Notes

1. ϵ is independent of $f_k$ provided that no dissipation is resolved, which corresponds to $f_{\epsilon }=1$

2. https://turbmodels.larc.nasa.gov/nasahump_val.html