A framework of the transport model for high-order eddy viscosity tensor in two-dimensional turbulent flow

\nameXingguang Zhou, Dalin Zhang*, Wenxi Tian, Guanghui Su and Suizheng Qiu CONTACT Prof. Dalin Zhang Email: [email protected] State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, 710049, Xi’an, Shaanxi, China

Abstract

Motivated by the concept of eddy viscosity tensor in improved Boussinesq hypothesis, a transport model of high-order eddy viscosity tensor in 2D-3C turbulence structure is derived from the second-order moment model by tensorial analysis. Compared with the previous high-order eddy viscosity tensor mathematical model [1, 2, 3] for turbulent flow under special operating conditions, we have obtained a general high-order eddy viscosity tensor transport model, which can show the anisotropy and transport characteristics of Reynolds stress. The transport equation of a high-order eddy viscosity tensor includes the transient term, convection term, production term, diffusion term, positive-definite source term, and dissipation term. In addition, when the components of high-order eddy viscosity tensor are positive or negative, the physical meaning of the source terms is also discussed. Two numerical simulation cases are also interpreted to verify the model. First, we apply the transport model of high-order eddy viscosity to the numerical calculation analysis of two-dimensional straight channel turbulent flow, and the results show that this model can effectively characterize the anisotropy of Reynolds stress. With a combination of analytical and numerical methods, we analyze the flow fields and the components of high-order eddy viscosity tensor to investigate the evolution law of high-order eddy viscosity tensor. The characteristics of the principal components and tensorial invariance of dimensionless Reynolds stress anisotropy tensor of two-dimensional straight channel turbulent flow also have been investigated with the Lumley triangle [4, 5]. The results show that this model has a good agreement with the DNS calculation results for the characteristics of dimensionless Reynolds stress anisotropy tensor. Then, we make the numerical investigation for a two-dimensional planar asymmetric diffuser and compare the results with the experiments of Obi et.al [6]. The transport model can resolve the separation bubble. The detachment point, attachment point, and separation bubble length are also compared with the experimental measurements, and the results of LES, are shown in good agreement. Finally, other discussions between 2D-3C turbulence structure and EV-RSM are also introduced, which may supply the further calibration of this model in the future. Such a transport model of high-order eddy viscosity tensor is expected to be useful in the theoretical and numerical analysis of high-order eddy viscosity tensor with improved Boussinesq hypothesis for turbulent flow.

keywords:

High-order eddy viscosity tensor, improved Boussinesq hypothesis, tensorial analysis, turbulence modeling

^†^†articletype: ORIGINAL MANUSCRIPT

1 Introduction

An incompressible fluid flow is described by the governing equations, as

	$\frac{\partial u_{i}}{\partial x_{i}}=0,$		(1a)
	$\rho\frac{\partial u_{i}}{\partial t}+\rho u_{j}\frac{\partial u_{i}}{\partial x_{j}}=\frac{\partial\tau_{ij}}{\partial x_{j}}+f_{i}.$		(1b)

In which $u_{i}$ is the velocity, $\rho$ is the density of the fluid, $\tau_{ij}$ is a second-order symmetric tensor which represents the forces acting on surfaces, $f_{i}$ is the external body forces such as gravity. For a Newtonian fluid, Navier and Stokes proposed a series of assumptions about modeling the stress tensor $\tau_{ij}$ . The constitutive equation of $\tau_{ij}$ is derived as

\tau_{ij}=-p\delta_{ij}+c_{ijkl}\frac{\partial u_{k}}{\partial x_{l}},

(2)

where $p$ is the pressure, $\tau_{ij}^{V}=c_{ijkl}\frac{\partial u_{k}}{\partial x_{l}}$ is the viscous stress tensor, which is a linear homogeneous tensor function of the local velocity gradient. $c_{ijkl}$ is a fourth-order tensor used to characterize the molecular viscosity of the fluid and the index of $i$ and $j$ are symmetric. For general fluids, $c_{ijkl}$ is usually reduced to an isotropic fourth-order tensor, with symmetric index of $i$ and $j$

c_{ijkl}=\nu\delta_{ij}\delta_{kl}+\mu\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)

(3)

in which $\nu$ and $\mu$ are both scalars, $\delta_{ij}$ is the Kronecker symbol. Substituting Equation (3) into Equation (2), the constitutive relation of incompressible fluid flow can be derived as

\tau_{ij}=-p\delta_{ij}+2\mu S_{ij},

(4)

where $\mu$ is the dynamic viscosity, and $S_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)$ is the rate of the strain tensor. Therefore, the Navier-Stokes equation for incompressible fluid flow can be derived as

\frac{\partial u_{i}}{\partial t}+u_{j}\frac{\partial u_{i}}{\partial x_{j}}=-\frac{1}{\rho}\frac{\partial p}{\partial x_{i}}+\nu\frac{\partial}{\partial x_{j}}\frac{\partial u_{i}}{\partial x_{j}}+g_{i},

(5)

in which $\nu$ is the kinematic viscosity, and $g_{i}$ is the gravitational acceleration.

For turbulent flow, the Reynolds time-averaging method is usually used to get the RANS (Reynolds-Averaged Navier-Stokes) equation to describe the motion of fluids in large-scale turbulent flow

\frac{\partial\overline{u}_{i}}{\partial t}+\overline{u}_{j}\frac{\partial\overline{u}_{i}}{\partial x_{j}}=-\frac{1}{\rho}\frac{\partial\overline{p}}{\partial x_{i}}+\nu\frac{\partial}{\partial x_{j}}\frac{\partial\overline{u}_{i}}{\partial x_{j}}-\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial x_{j}}+g_{i},

(6)

where $\overline{u}_{i}$ is the time-averaging velocity, $u^{\prime}_{i}=u_{i}-\overline{u}_{i}$ is the fluctuating velocity.

In RANS, an unclosed term called Reynolds stress $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ has occurred. Analogous to the constitutive relationship between stress tensor and rate of strain tensor in laminar flow as Equation (4), the eddy viscosity model (EVM) is proposed by Boussinesq hypothesis [7] to model the Reynolds stress

-\overline{u^{\prime}_{i}u^{\prime}_{j}}=-\frac{2}{3}k\delta_{ij}+\nu_{t}\left(\frac{\partial\overline{u}_{i}}{\partial x_{j}}+\frac{\partial\overline{u}_{j}}{\partial x_{i}}\right),

(7)

in which $k=\frac{1}{2}\overline{u^{\prime}_{i}u^{\prime}_{i}}$ is the turbulent kinetic energy, $\nu_{t}$ is a scalar called the eddy kinematic viscosity. Equation (7) relates Reynolds stress to time-averaging quantities. According to the dimensional analysis, many two-equations models, such as standard $k-\varepsilon$ [8] and standard $k-\omega$ [9], have been proposed to solve $\nu_{t}$ . However, in these two-equations models, only partly transportation property of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ can be considered through $k$ equation and $\varepsilon$ equation or $\omega$ equation [10]. The anisotropic, relaxation, and historical effect of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ are not well characterized.

For anisotropic of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ , the NLEVM (Non-Linear Eddy Viscosity Model) [11] is proposed. In NLEVM, $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ is modeled with a series of tensor basis [12, 13]

\overline{u^{\prime}_{i}u^{\prime}_{j}}=\frac{2}{3}k\delta_{ij}+\sum_{n=1}^{N}\alpha^{\prime}_{n}T_{ij}^{\left(n\right)},

(8)

where $\alpha^{\prime}_{n}$ is the expansion coefficients, and $T_{ij}^{\left(n\right)}\left(n=1,\cdots,N\right)$ is a given tensor basis. According to the definition of dimensionless Reynolds stress anisotropy tensor $b_{ij}=\overline{u^{\prime}_{i}u^{\prime}_{j}}/\overline{u^{\prime}_{k}u^{\prime}_{k}}-\frac{1}{3}\delta_{ij}$ , Equation (8) can be simplified to

b_{ij}=\frac{\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\overline{u^{\prime}_{k}u^{\prime}_{k}}}=\sum_{n=1}^{N}\alpha^{\prime}_{n}T_{ij}^{\left(n\right)}.

(9)

With the assumption of the functional dependency of $b_{ij}=b_{ij}\left(S_{kl},W_{kl},k/\varepsilon\right)$ , in which $W_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial u_{j}}{\partial x_{i}}\right)$ is the vorticity tensor, $T_{ij}^{\left(n\right)}$ can be solved by isotropic tensor function transformation and Cayley-Hamilton theorem. The expansion coefficients $\alpha^{\prime}_{n}$ need to be calibrated by DNS data or experimental data [11]. NLEVM has been interpreted in many research works, such as Shih quadratic model [14] and Lien cubic model [15]. However, for the tensor basis $T_{ij}^{\left(n\right)}$ , some assumptions are made and cannot guarantee the choice of the best representation basis [11]. Regarding the property of NLEVM, the more terms of tensor series expansion, the more accurate the calculation result is. However, more terms of tensor series expansion will make $\alpha^{\prime}_{n}$ more difficult to be calibrated by DNS or experimental measurements, and the oscillation of non-physical phenomena may occur [16].

In the aspect of RANS, Reynolds stress is a physical quantity with large-scale behavior that is characterized by the flow field. And the stress tensor in laminar flow is a molecular-scale behavior. Therefore, Reynolds stress does not have the same instantaneous property as the stress tensor in laminar flow. To characterize the anisotropy, historical effect, and relaxation property of Reynolds stress and to avoid the assumptions and higher-order tensor basis expansion in NLEVM. The concept of eddy viscosity tensor [1, 2, 17, 3, 18] is proposed.

To date, many researches on high-order eddy viscosity tensor have been proposed. Dubrulle and Frisch [1] use the multiscale formalism theory to discuss the eddy viscosity tensor for the incompressible flow of arbitrary dimensionality subject to forcing periodic in space and time, and the explicit expressions of eddy viscosity tensor are given for basic flow with low Reynolds numbers, and when the basic flow is layered. Gama et al. [2] make the theoretical and numerical research of negative eddy viscosity in isotropic forced two-dimensional flow. But the flow is also assumed that the basic flow in space-time periodic with parity-invariance. Wirth et al. [3] give the detailed theoretical and numerical results for eddy viscosity tensor of three-dimensional forced spatially periodic incompressible flow with weak large-scale perturbation theory. Wang and Jiang [18] consider the high-order eddy viscosity tensor for the modeling of shear turbulence. The high-order eddy viscosity tensor is inversely solved by using Moore-Penrose generalized inverse matrix theory in abstract algebra for Reynolds stress. However, this method is not good for the closure of turbulence model. And the high-order eddy viscosity tensor directly acts on the rate of strain tensor rather than the velocity gradient tensor, which is a simplified form of the high-order eddy viscosity tensor.

Most of the previous theoretical works are aimed at the simple and basic turbulent flow. Both the high-order eddy viscosity tensor and the eddy viscosity coefficient are the modeling methods to characterize Reynolds stress. Moreover, Reynolds stress is a conserved quantity, which has a conservation equation and is known as the second-order moment transport equation. Hence, we consider that the high-order eddy viscosity tensor can also be regarded as a conserved quantity, which means that there should be a conservation equation for high-order eddy viscosity tensor. This conservation equation should be an explicit fourth-order tensor transport equation and contains the transient term, convection term, generation term, dissipation term, and other source terms to reflect the physical evolution of high-order eddy viscosity tensor. In this article, we carry out the study of this transport equation for high-order eddy viscosity tensor and try to explain the evolution of high-order eddy viscosity in the aspect of physics.

The rest of this article is organized as follows. In Section 2, we build the constitutive relation between Reynolds stress and velocity gradient with the fourth-order eddy viscosity tensor in 2D-3C turbulence structure. Then we investigate the second-order moment transport equation of Reynolds stress and make the closure by LRR-IP [19] method. The constitutive relation is regarded as the tensor function. With the mathematical transformation and tensorial analysis of the modeled second-order moment transport equation, the transport equation of high-order eddy viscosity tensor in 2D-3C turbulence structure is obtained, and each term of the transport equation is analyzed in the aspect of physics. In Section 3, numerical analysis and model validation are carried out. First, we study the evolution of high-order eddy viscosity tensor in the two-dimensional straight channel turbulence. The characteristics of the principal components and tensorial invariance of dimensionless Reynolds stress anisotropy tensor are also investigated with the Lumley triangle [4, 5]. For complex turbulent flow, the two-dimensional turbulent flow in an asymmetric planar diffuser also be calculated, and the results are validated by the experimental measurements. With the transport model of high-order eddy viscosity tensor, the distribution of each independent component of high-order eddy viscosity tensor in an asymmetric planar diffuser is also obtained. In Section 4, we review the assumptions that are made in the constitutive relation between Reynolds stress and velocity gradient with high-order eddy viscosity tensor in 2D-3C turbulence structure. The transport equation in a strict form of eddy viscosity coefficient is derived, and the relation of the Spalart-Allmaras one-equation model [20] is also introduced. Finally, we summarize our results and derivations of this framework of the transport model for high-order eddy viscosity tensor. The outlook and potential extensions are also discussed in this study.

2 Mathematical method and tensorial analysis

2.1 Constitutive relation and assumptions

Analogous to the constitutive of laminar flow as Equation (2), $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ can be modeled by

-\overline{u^{\prime}_{i}u^{\prime}_{j}}=-\frac{2}{3}k\delta_{ij}+\zeta_{ijkl}\frac{\partial u_{l}}{\partial x_{k}},

(10)

in which $\zeta_{ijkl}$ is a fourth-order eddy viscosity tensor.

Considering a statistically two-dimensional flow in which statistics are independent of $x_{3}$ , and which is statistically invariant under reflections of the $x_{3}$ coordinate axis. The PDF (Probability Density Function) of velocity can be implied as

	$\frac{\partial f}{\partial x_{3}}=0,$		(11a)
	$f\left(u_{1},u_{2},u_{3};x_{1},x_{2},x_{3},t\right)=f\left(u_{1},u_{2},-u_{3};x_{1},x_{2},x_{3},t\right).$		(11b)

At $x_{3}=0$ , we can get $\overline{u}_{3}=-\overline{u}_{3}$ and $\overline{u}_{3}=0$ . According to the symmetry, we also can get $\overline{u^{\prime}_{1}u^{\prime}_{3}}=\overline{u^{\prime}_{2}u^{\prime}_{3}}=0$ . Hence for a statistically two-dimensional flow, $\overline{u}_{3}=0$ and Reynolds stress can be derived as

R_{ij}=\overline{u^{\prime}_{i}u^{\prime}_{j}}=\begin{pmatrix}\overline{u^{\prime}_{1}u^{\prime}_{1}}&\overline{u^{\prime}_{1}u^{\prime}_{2}}&0\\ \overline{u^{\prime}_{2}u^{\prime}_{1}}&\overline{u^{\prime}_{2}u^{\prime}_{2}}&0\\ 0&0&\overline{u^{\prime}_{3}u^{\prime}_{3}}\end{pmatrix}.

(12)

Considering the linear EVM constitutive relation in Equation (7), for a two-dimensional turbulent flow, the components of normal Reynolds viscous stress of the $x_{3}$ coordinate axis can be obtained

R_{33}^{V}=-\overline{u^{\prime}_{3}u^{\prime}_{3}}+\frac{2}{3}k=2\nu_{t}\frac{\partial\overline{u}_{3}}{\partial x_{3}}=0.

(13)

Therefore, considering the linear EVM constitutive relation, an assumption is made of the transport model of $\zeta_{ijkl}$ , that

R_{33}^{V}=-\overline{u^{\prime}_{3}u^{\prime}_{3}}+\frac{2}{3}k=\zeta_{33kl}\frac{\partial\overline{u}_{l}}{\partial x_{k}}=0.

(14)

Wang and Jiang [18] also adopted the same assumption as shown in Equation (14) when using the simplified high-order eddy viscosity tensor to calculate the two-dimensional turbulent flow. Hence for two-dimensional turbulent flow, the Reynolds viscous stress in the transport model of $\zeta_{ijkl}$ can be obtained

R_{ij}^{V}=-\overline{u^{\prime}_{i}u^{\prime}_{j}}+\frac{2}{3}k\delta_{ij}=\begin{pmatrix}\zeta_{11kl}\frac{\partial\overline{u}_{l}}{\partial x_{k}}&\zeta_{12kl}\frac{\partial\overline{u}_{l}}{\partial x_{k}}&\\ \zeta_{21kl}\frac{\partial\overline{u}_{l}}{\partial x_{k}}&\zeta_{22kl}\frac{\partial\overline{u}_{l}}{\partial x_{k}}&\\ &&0\end{pmatrix}.

(15)

From Equation (10), by contracting the index of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ , it is obvious that $R_{ij}^{V}$ is traceless. Therefore, we can get an identity relation

\zeta_{11kl}\frac{\partial\overline{u}_{l}}{\partial x_{k}}+\zeta_{22kl}\frac{\partial\overline{u}_{l}}{\partial x_{k}}=0.

(16)

Based on the above derivation, $\zeta_{ijkl}$ has 12 independent components because of the symmetric of index $i$ and $j$ , that $\zeta_{12kl}=\zeta_{21kl}$ .

Therefore, the constitutive relation and mathematical form of high-order eddy viscosity tensor are determined in 2D-3C turbulence structure. Based on the derivation in this section, the transport model of high-order eddy viscosity tensor is studied. And the assumptions which are made for constitutive relation will be also discussed in the following sections.

2.2 Closure of second-order moment model and derivation of transport equation of $\zeta_{ijkl}$

From Equation (10), only anisotropy of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ can be characterized. Spalart and Allmaras [20] develop a one-equation model, which builds a transport equation of $\nu_{t}$ in an empirical way. Spalart-Allmaras one-equation model still uses the linear EVM constitutive relation, which characterizes the isotropic turbulence. However, the transport equation of $\nu_{t}$ reveals the relaxation and history effect of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ .

In the following sections, we need to point out the following two points to explain the derivation concisely and clearly:

1.

If there are no special instructions, time-averaging velocity $\overline{u}_{i}$ will be simplified to mark as $u_{i}$ , time-averaging pressure $\overline{p}$ will be simplified to mark as $p$ .
2.

The derivation of this article is carried out in the Cartesian coordinate system, which does not need to make the clear specification of contravariant or covariant to the tensor bases and components.

Researching Equation (10) in the aspects of continuum mechanics and tensorial analysis. From the definition of finite differential of tensor and the derivative of tensor function, since the $\frac{\partial u_{l}}{\partial x_{k}}$ is not a symmetric tensor, that Equation (10) can be written in

\frac{\partial R_{ij}^{V}}{\partial\left(\frac{\partial u_{l}}{\partial x_{k}}\right)}\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}=\frac{\partial\left(-\overline{u^{\prime}_{i}u^{\prime}_{j}}+\frac{2}{3}k\delta_{ij}\right)}{\partial\left(\frac{\partial u_{l}}{\partial x_{k}}\right)}\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}=\zeta_{ijkl}\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l},

(17)

in which $\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}$ is the base of the fourth-order tensor in the Cartesian coordinate system.

Equation (17) shows that we regard the constitutive relation of Equation (10) as the form of tensor function

\mathbf{R}^{V}=\mathbf{R}^{V}\left(\nabla\mathbf{u}\right)=\boldsymbol{\zeta}:\nabla\mathbf{u}.

(18)

$\mathbf{R}^{V}$ is the function of $\nabla\mathbf{u}$ , where $\nabla\mathbf{u}$ is the argument. In continuum mechanics, such as solid mechanics and the theory of elasticity, the constitutive relation between elastic stress tensor and rate of strain tensor is similar to Equation (10). However, the explicit correlation of ordinary material between elastic stress tensor and rate of strain tensor can be measured by experiments [21, 22]. In turbulent fluid flow, the explicit correlation of Reynolds viscous stress field $R_{ij}^{V}$ and velocity gradient fields $\frac{\partial u_{l}}{\partial x_{k}}$ is difficult to measure because of the nonlinear relationship of RANS. For complex turbulent flow, the analytic solution is basically impossible to be obtained. Hence, motivated by Spalart-Allmaras one-equation model to build the transport equation of $\nu_{t}$ . We are also able to build a transport model of high-order eddy viscosity tensor $\zeta_{ijkl}$ to preliminary reveal its evolution law. Different from Spalart-Allmaras one-equation model which is built in an empirical way. We research the second-order moment transport model, to get the transport model of $\zeta_{ijkl}$ in a theoretical way.

The strict form of the second-order moment transport model can be derived from the turbulent fluctuating momentum equation by Reynolds time-averaging

$\displaystyle\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial t}+u_{k}\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial x_{k}}$	$\displaystyle=\underbrace{-\left(\overline{u^{\prime}_{i}u^{\prime}_{k}}\frac{\partial u_{j}}{\partial x_{k}}+\overline{u^{\prime}_{j}u^{\prime}_{k}}\frac{\partial u_{i}}{\partial x_{k}}\right)}_{P_{ij}}$	(19)
	$\displaystyle\underbrace{+\overline{\frac{p^{\prime}}{\rho}\left(\frac{\partial u^{\prime}_{j}}{\partial x_{i}}+\frac{\partial u^{\prime}_{i}}{\partial x_{j}}\right)}}_{\Phi_{ij}}$
	$\displaystyle\underbrace{-\frac{\partial}{\partial x_{k}}\left(\overline{u^{\prime}_{i}u^{\prime}_{j}u^{\prime}_{k}}-\nu\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial x_{k}}+\frac{\delta_{ik}}{\rho}\overline{u^{\prime}_{j}p^{\prime}}+\frac{\delta_{jk}}{\rho}\overline{u^{\prime}_{i}p^{\prime}}\right)}_{D_{ij}}$
	$\displaystyle-\underbrace{2\nu\overline{\left(\frac{\partial u^{\prime}_{i}}{\partial x_{k}}\right)\left(\frac{\partial u^{\prime}_{j}}{\partial x_{k}}\right)}}_{\varepsilon_{ij}},$

where $P_{ij}$ is the production term, $\Phi_{ij}$ is the pressure rate of strain term, $D_{ij}$ is the diffusion term, $\varepsilon_{ij}$ is the dissipation term. By Green-Gauss theorem, Equation (19) can be written in the integral form with vector notation

\frac{\partial}{\partial t}\int_{V}\mathbf{R}dV+\oint_{A}\mathbf{n}\cdot\mathbf{u}\otimes\mathbf{R}dA=\int_{V}\mathbf{P}dV+\int_{V}\mathbf{\Phi}dV+\oint_{A}\mathbf{n}\cdot\mathbf{D}dA-\int_{V}\boldsymbol{\varepsilon}dV,

(20)

where $\mathbf{R}=\overline{u^{\prime}_{i}u^{\prime}_{j}}\mathbf{e}_{i}\mathbf{e}_{j}$ is the Reynolds stress, $\mathbf{n}$ is the unit normal vector of the surface. Substituting Equation (10) into Equation (20), we can get the transport equation of Reynolds viscous stress tensor $\mathbf{R}^{V}=-\overline{u^{\prime}_{i}u^{\prime}_{j}}\mathbf{e}_{i}\mathbf{e}_{j}+\frac{2}{3}k\delta_{ij}\mathbf{e}_{i}\mathbf{e}_{j}$

	$\displaystyle\frac{\partial}{\partial t}\int_{V}\mathbf{R}^{V}dV-\frac{\partial}{\partial t}\int_{V}\frac{2}{3}k\mathbf{I}dV$	$\displaystyle+\oint_{A}\mathbf{n}\cdot\mathbf{u}\otimes\mathbf{R}^{V}dA-\oint_{A}\mathbf{n}\cdot\mathbf{u}\otimes\frac{2}{3}k\mathbf{I}dA$		(21)
		$\displaystyle=-\int_{V}\mathbf{P}dV-\int_{V}\mathbf{\Phi}dV-\oint_{A}\mathbf{n}\cdot\mathbf{D}dA+\int_{V}\boldsymbol{\varepsilon}dV,$		(21)

in which $\mathbf{I}=\delta_{ij}\mathbf{e}_{i}\mathbf{e}_{j}$ . In Equation (21), it is obvious that turbulent kinetic energy transport equation is included. By contracting the index of $i$ and $j$ for Equation (19), we can obtain the transport equation of $k$

\frac{\partial k}{\partial t}+u_{m}\frac{\partial k}{\partial x_{m}}=\underbrace{-\overline{u^{\prime}_{i}u^{\prime}_{m}}\frac{\partial u_{i}}{\partial x_{m}}}_{P_{k}}\underbrace{-\frac{\partial}{\partial x_{m}}\left(\frac{\overline{p^{\prime}u^{\prime}_{m}}}{\rho}+\overline{u^{\prime}_{m}u^{\prime}_{i}u^{\prime}_{i}}-\nu\frac{\partial k}{\partial x_{m}}\right)}_{D_{k}}-\underbrace{\nu\overline{\left(\frac{\partial u^{\prime}_{i}}{\partial x_{m}}\right)\left(\frac{\partial u^{\prime}_{i}}{\partial x_{m}}\right)}}_{\varepsilon},

(22)

in which $P_{k}$ is the production term, $D_{k}$ is the diffusion term, $\varepsilon$ is the turbulent kinetic energy dissipation term. Equation (22) can also be written in the integral form with vector notation

\frac{\partial}{\partial t}\int_{V}kdV+\oint_{A}\mathbf{n}\cdot\mathbf{u}kdA=\int_{V}P_{k}dV+\oint_{A}\mathbf{n}\cdot D_{k}dA-\int_{V}\varepsilon dV.

(23)

Tensor multiplying $\frac{2}{3}\mathbf{I}$ on both sides to Equation (23), we can obtain

\frac{\partial}{\partial t}\int_{V}\frac{2}{3}k\mathbf{I}dV+\oint_{A}\mathbf{n}\cdot\mathbf{u}\otimes\frac{2}{3}k\mathbf{I}dA=\int_{V}P_{k}\otimes\frac{2}{3}\mathbf{I}dV+\oint_{A}\mathbf{n}\cdot D_{k}\otimes\frac{2}{3}\mathbf{I}dA-\int_{V}\varepsilon\otimes\frac{2}{3}\mathbf{I}dV.

(24)

Add Equation (21) and Equation (24), the transport equation of Reynolds viscous stress $\mathbf{R}^{V}$ in strict form is obtained

$\displaystyle\frac{\partial}{\partial t}\int_{V}\mathbf{R}^{V}dV+\oint_{A}\mathbf{n}\cdot\mathbf{u}\otimes\mathbf{R}^{V}dA$	$\displaystyle=$	(25)
	$\displaystyle\underbrace{-\int_{V}\mathbf{P}dV-\int_{V}\mathbf{\Phi}dV-\oint_{A}\mathbf{n}\cdot\mathbf{D}dA+\int_{V}\boldsymbol{\varepsilon}dV}_{Source\ term\ by\ second-order\ moment\ equation}$
	$\displaystyle\underbrace{+\int_{V}P_{k}\otimes\frac{2}{3}\mathbf{I}dV+\oint_{A}\mathbf{n}\cdot D_{k}\otimes\frac{2}{3}\mathbf{I}dA-\int_{V}\varepsilon\otimes\frac{2}{3}\mathbf{I}dV}_{Source\ term\ by\ turbulent\ kinetic\ energy\ equation}.$

Equation (25) is the strict form that contains unclosed terms. Hence, we need to model these unclosed terms before making the tensorial analysis of $\zeta_{ijkl}$ .

The pressure rate of strain term $\Phi_{ij}$ modeling is crucial and the subject of extensive research [23]. LRR-IP model [19] is used in this article, that $\Phi_{ij}$ can be divided into three terms

\Phi_{ij}=\overline{\frac{p^{\prime}}{\rho}\left(\frac{\partial u^{\prime}_{j}}{\partial x_{i}}+\frac{\partial u^{\prime}_{i}}{\partial x_{j}}\right)}=\Phi_{ij}^{r}+\Phi_{ij}^{s}+\Phi_{ij}^{w},

(26)

where $\Phi_{ij}^{r}$ is the rapid pressure strain term, $\Phi_{ij}^{s}$ is the slow pressure strain term, and $\Phi_{ij}^{w}$ is the wall reflection term. With rapid-distortion theory and Rotta’s model [19, 23], the basic model for $\Phi_{ij}$ can be obtained

\Phi_{ij}=\underbrace{-2C_{R}\varepsilon b_{ij}}_{\Phi_{ij}^{r}}\underbrace{-C_{2}\left(P_{ij}-\frac{2}{3}P_{k}\delta_{ij}\right)}_{\Phi_{ij}^{s}},

(27)

where $C_{R}=1.8,C_{2}=0.6$ are empirical constants. The wall reflection term $\Phi_{ij}^{w}$ is not considered in this basic model.

For the diffusion tensor of Reynolds stress, Shir model [24] is adopted in this article

\overline{u^{\prime}_{i}u^{\prime}_{j}u^{\prime}_{m}}=-C_{S}\frac{k^{2}}{\varepsilon}\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial x_{m}},

(28)

in which $C_{S}=0.09$ is an empirical constant. The fluctuating pressure-velocity correlation term $\frac{\delta_{ik}}{\rho}\overline{u^{\prime}_{j}p^{\prime}}+\frac{\delta_{jk}}{\rho}\overline{u^{\prime}_{i}p^{\prime}}$ is ignored [19].

In LRR-IP model, the dissipation tensor of Reynolds stress adopts the homogeneous approximation

\varepsilon_{ij}=2\nu\overline{\frac{\partial u^{\prime}_{i}}{\partial x_{k}}\frac{\partial u^{\prime}_{j}}{\partial x_{k}}}=\frac{2}{3}\varepsilon\delta_{ij}.

(29)

For the diffusion term $D_{k}$ in turbulent kinetic energy equation, linear eddy diffusion model [25] is adopted, and ignore the fluctuating pressure-velocity correlation term $\overline{p^{\prime}u^{\prime}_{m}}/\rho$

\overline{u^{\prime}_{m}u^{\prime}_{i}u^{\prime}_{i}}=-C_{S}\frac{k^{2}}{\varepsilon}\frac{\partial k}{\partial x_{m}}.

(30)

Substitute Equation (26) to (30) into Equation (25), the modeled transport equation of Reynolds viscous stress tensor $R_{ij}^{V}$ can be obtained

\frac{\partial}{\partial t}\int_{V}R_{ij}^{V}dV+\oint_{A}n_{m}u_{m}R_{ij}^{V}dA=P_{ij}^{M}+D_{ij}^{M}+\Phi_{R,ij}^{M}+\Phi_{S,ij}^{M}+P_{k,ij}^{M}+D_{k,ij}^{M},

(31)

where $P_{ij}^{M}$ , $D_{ij}^{M}$ , $\Phi_{R,ij}^{M}$ , $\Phi_{S,ij}^{M}$ , $P_{k,ij}^{M}$ , and $D_{k,ij}^{M}$ are the modeled terms of Equation (25), which can be shown as

	$P_{ij}^{M}=-\int_{V}\left(R_{im}^{V}\frac{\partial u_{j}}{\partial x_{m}}+R_{jm}^{V}\frac{\partial u_{i}}{\partial x_{m}}\right)dV+\int_{V}\left(\frac{2}{3}k\delta_{im}\frac{\partial u_{j}}{\partial x_{m}}+\frac{2}{3}k\delta_{jm}\frac{\partial u_{i}}{\partial x_{m}}\right)dV,$		(32a)
	$D_{ij}^{M}=\oint_{A}n_{m}\left(C_{S}\frac{k^{2}}{\varepsilon}+\nu\right)\frac{\partial R_{ij}^{V}}{\partial x_{m}}dA-\oint_{A}n_{m}\left(C_{S}\frac{k^{2}}{\varepsilon}+\nu\right)\frac{\partial}{\partial x_{m}}\left(\frac{2}{3}k\delta_{ij}\right)dA,$		(32b)
	$\Phi_{ij}^{R}=-\int_{V}C_{R}\frac{\varepsilon}{k}R_{ij}^{V}dV,$		(32c)
	$\Phi_{ij}^{S}=\int_{V}C_{2}\left[\left(R_{im}^{V}\frac{\partial u_{j}}{\partial x_{m}}+R_{jm}^{V}\frac{\partial u_{i}}{\partial x_{m}}-\frac{2}{3}k\delta_{im}\frac{\partial u_{j}}{\partial x_{m}}-\frac{2}{3}k\delta_{jm}\frac{\partial u_{i}}{\partial x_{m}}\right)-\frac{2}{3}\delta_{ij}R_{rs}^{V}\frac{\partial u_{r}}{\partial x_{s}}\right]dV,$		(32d)
	$P_{k,ij}^{M}=\int_{V}\frac{2}{3}\delta_{ij}R_{rs}^{V}\frac{\partial u_{r}}{\partial x_{s}}dV,$		(32e)
	$D_{k,ij}^{M}=\oint_{A}n_{m}\frac{2}{3}\delta_{ij}\left(C_{k}\frac{k^{2}}{\varepsilon}+\nu\right)\frac{\partial k}{\partial x_{m}}dA,$		(32f)

where $C_{k}=C_{S}=0.09$ [19, 25]. After the modeling of unclosed terms, we can find that the dissipation term does not appear in Equation (31). Then, we will make the tensorial derivation of each tensor function term in Equation (31) with the tensor function relation in Equation (18).

For the transient term in Equation (31), taking its derivative to $\nabla\mathbf{u}$ , we can get

\frac{d}{d\left(\nabla\mathbf{u}\right)}\frac{\partial}{\partial t}\int_{V}\mathbf{R}^{V}dV=\frac{\partial}{\partial t}\int_{V}\frac{d\mathbf{R}^{V}}{d\left(\nabla\mathbf{u}\right)}dV=\frac{\partial}{\partial t}\int_{V}\zeta_{ijkl}dV\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}.

(33)

For the convection term $\oint_{A}\mathbf{n}\cdot\mathbf{u}\otimes\mathbf{R}^{V}dA$ in Equation (31). Let $\mathbf{H}\left(\nabla\mathbf{u}\right)=\mathbf{n}\cdot\mathbf{u}\otimes\mathbf{R}^{V}$ and make the tensorial derivative to $\nabla\mathbf{u}$ , that we can get

\mathbf{H}^{\prime}\left(\nabla\mathbf{u}\right):\mathbf{C}=\mathbf{n}\cdot\left(\mathbf{u}^{\prime}:\mathbf{C}\right)\otimes\mathbf{R}^{V}+\mathbf{n}\cdot\mathbf{u}\otimes\left(\mathbf{R}^{V}\right)^{\prime}:\mathbf{C},

(34)

where $\mathbf{C}$ is an arbitrary second-order tensor from the definition of finite differential of tensor. In Equation (34), the unit normal vector of the surface $\mathbf{n}$ is a constant vector function field, therefore $d\mathbf{n}/d\left(\nabla\mathbf{u}\right)=\mathbf{0}$ . Since $\mathbf{u}=\mathbf{u}\left(\nabla\mathbf{u}\right)$ and $\mathbf{R}^{V}=\mathbf{R}^{V}\left(\nabla\mathbf{u}\right)$ are both the functions of $\nabla\mathbf{u}$ , according to the definition of the finite differential of tensor, an implicit term has occurred. However, the transport equation of $\zeta_{ijkl}$ need to be solved explicitly. We can know that $\mathbf{u}^{\prime}=d\mathbf{u}/d\left(\nabla\mathbf{u}\right)$ is a three-rank tensor with the quotient law of tensor, and the unit of dimension is in meter. It means that the velocity gradient can only work at a certain distance to make the velocity in change. If this distance is very small, we can regard $\mathbf{u}^{\prime}\left(\nabla\mathbf{u}\right)\approx\mathbf{0}$ through a very dense control volume and mesh. Hence, the convection term after the tensorial derivative can be obtained

\frac{d}{d\left(\nabla\mathbf{u}\right)}\oint_{A}\mathbf{n}\cdot\mathbf{u}\otimes\mathbf{R}^{V}dA=\oint_{A}n_{m}u_{m}\frac{\partial R_{ij}^{V}}{\partial\left(\frac{\partial u_{l}}{\partial x_{k}}\right)}dA\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}=\oint_{A}n_{m}u_{m}\zeta_{ijkl}dA\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}.

(35)

For $P_{ij}^{M}$ in Equation (31), we divide it into two parts

P_{ij}^{M}=-\underbrace{\int_{V}\left(R_{im}^{V}\frac{\partial u_{j}}{\partial x_{m}}+R_{jm}^{V}\frac{\partial u_{i}}{\partial x_{m}}\right)dV}_{\mathbf{H}_{1}}+\underbrace{\int_{V}\left(\frac{2}{3}k\delta_{im}\frac{\partial u_{j}}{\partial x_{m}}+\frac{2}{3}k\delta_{jm}\frac{\partial u_{i}}{\partial x_{m}}\right)dV}_{\mathbf{H}_{2}}.

(36)

For $\mathbf{H}_{1}=\mathbf{R}^{V}\cdot\nabla\mathbf{u}+\left(\nabla\mathbf{u}\right)^{T}\cdot\mathbf{R}^{V}$ , before we make the tensorial derivative of $\nabla\mathbf{u}$ , the mathematical transformation needs to be carried out. For two-dimensional incompressible turbulent flow, we can get,

	$H_{1,11}=R_{11}^{V}\frac{\partial u_{1}}{\partial x_{1}}+R_{12}^{V}\frac{\partial u_{1}}{\partial x_{2}}+R_{11}^{V}\frac{\partial u_{1}}{\partial x_{1}}+R_{12}^{V}\frac{\partial u_{1}}{\partial x_{2}}=2\left(R_{11}^{V}\frac{\partial u_{1}}{\partial x_{1}}+R_{12}^{V}\frac{\partial u_{1}}{\partial x_{2}}\right),$		(37a)
	$H_{1,12}=R_{11}^{V}\frac{\partial u_{2}}{\partial x_{1}}+\underbrace{R_{12}^{V}\frac{\partial u_{2}}{\partial x_{2}}+R_{21}^{V}\frac{\partial u_{1}}{\partial x_{1}}}_{\nabla\cdot\mathbf{u}=0}+R_{22}^{V}\frac{\partial u_{1}}{\partial x_{2}}=R_{11}^{V}\frac{\partial u_{2}}{\partial x_{1}}+R_{22}^{V}\frac{\partial u_{1}}{\partial x_{2}},$		(37b)
	$H_{1,21}=R_{22}^{V}\frac{\partial u_{1}}{\partial x_{2}}+\underbrace{R_{21}^{V}\frac{\partial u_{1}}{\partial x_{1}}+R_{12}^{V}\frac{\partial u_{2}}{\partial x_{2}}}_{\nabla\cdot\mathbf{u}=0}+R_{11}^{V}\frac{\partial u_{2}}{\partial x_{1}}=R_{22}^{V}\frac{\partial u_{1}}{\partial x_{2}}+R_{11}^{V}\frac{\partial u_{2}}{\partial x_{1}},$		(37c)
	$H_{1,22}=R_{21}^{V}\frac{\partial u_{2}}{\partial x_{1}}+R_{22}^{V}\frac{\partial u_{2}}{\partial x_{2}}+R_{21}^{V}\frac{\partial u_{2}}{\partial x_{1}}+R_{22}^{V}\frac{\partial u_{2}}{\partial x_{2}}=2\left(R_{21}^{V}\frac{\partial u_{2}}{\partial x_{1}}+R_{22}^{V}\frac{\partial u_{2}}{\partial x_{2}}\right).$		(37d)

Equation (37a) to (37d) are identity transformations with the continuity equation of incompressible flow and the symmetric property of $\mathbf{R}^{V}$ . For Equation (37b) and (37c), consider $\mathbf{R}^{V}$ is traceless, which is interpreted in Equation (16). Hence Equation (37b) and (37c) can be rewritten in

	$H_{1,12}=R_{11}^{V}\frac{\partial u_{2}}{\partial x_{1}}+R_{22}^{V}\frac{\partial u_{1}}{\partial x_{2}}=R_{11}^{V}\left(\frac{\partial u_{2}}{\partial x_{1}}-\frac{\partial u_{1}}{\partial x_{2}}\right),$		(38a)
	$H_{1,21}=R_{22}^{V}\frac{\partial u_{1}}{\partial x_{2}}+R_{11}^{V}\frac{\partial u_{2}}{\partial x_{1}}=R_{22}^{V}\left(\frac{\partial u_{1}}{\partial x_{2}}-\frac{\partial u_{2}}{\partial x_{1}}\right).$		(38b)

Therefore, by re-summarizing $\mathbf{H}_{1}=\mathbf{R}^{V}\cdot\nabla\mathbf{u}+\left(\nabla\mathbf{u}\right)^{T}\cdot\mathbf{R}^{V}$ , we can get

H_{1,ij}=R_{im}^{V}\frac{\partial u_{j}}{\partial x_{m}}+R_{jm}^{V}\frac{\partial u_{i}}{\partial x_{m}}=2\delta_{ij}R_{\underline{i}m}^{V}\frac{\partial u_{\underline{i}}}{\partial x_{m}}+\epsilon_{ij}\omega_{Z}R_{\underline{i}\underline{i}}^{V},

(39)

in which $\epsilon_{ij}$ is a two-dimensional Eddington tensor, that $\epsilon_{11}=\epsilon_{22}=0$ , $\epsilon_{12}=1$ , and $\epsilon_{21}=-1$ . $\omega_{Z}=\frac{\partial u_{2}}{\partial x_{1}}-\frac{\partial u_{1}}{\partial x_{2}}$ is the vorticity of the $x_{3}$ coordinate axis. It should be noted that, in Equation (39), a horizontal line has been added below the index $i$ . This means that index $i$ is neither free nor dummy, only plays the following role, and the relevant definitions can be referred to tensor algebra and analysis [26, 27, 28].

After the above identical mathematical transformation, the tensorial derivative of $\mathbf{H}_{1}$ to $\nabla\mathbf{u}$ can be expressed uniquely and explicitly, as

\mathbf{H}^{\prime}_{1}=\int_{V}\left[\left(2\delta_{ij}R_{\underline{i}m}^{V}\delta_{\underline{i}l}\delta_{mk}+2\delta_{ij}\zeta_{\underline{i}mkl}\frac{\partial u_{\underline{i}}}{\partial x_{m}}\right)+\left(\epsilon_{ij}R_{\underline{i}\underline{i}}^{V}\Omega_{kl}+\epsilon_{ij}\omega_{Z}\zeta_{\underline{i}\underline{i}kl}\right)\right]dV\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l},

(40)

where $\mathbf{\Omega}$ is a second-rank tensor, expressed as

\omega_{Z}=\mathbf{\Omega}:\nabla\mathbf{u},\nabla\mathbf{u}=\begin{pmatrix}\frac{\partial u_{1}}{\partial x_{1}}&\frac{\partial u_{2}}{\partial x_{1}}\\ \frac{\partial u_{1}}{\partial x_{2}}&\frac{\partial u_{2}}{\partial x_{2}}\end{pmatrix}\Rightarrow\frac{d\omega_{Z}}{d\left(\nabla\mathbf{u}\right)}=\mathbf{\Omega}=\begin{pmatrix}0&1\\ -1&0\end{pmatrix}.

(41)

For the $\mathbf{H}_{2}$ part in Equation (36), we can find that the turbulent kinetic energy $k$ appears. Since turbulent kinetic energy $k$ and turbulent kinetic dissipation rate $\varepsilon$ are not explicit tensor functions of $\nabla\mathbf{u}$ . We can consider such a condition, for the steady turbulent kinetic energy transport equation of incompressible fluid flow

\nabla\cdot\left(\mathbf{u}k\right)=P_{k}+\nabla\cdot\left(\Gamma_{k}\nabla k\right)-\varepsilon,

(42)

where $\Gamma_{k}=\nu+\nu_{t}^{k}$ is considered as the general diffusion coefficient of $k$ equation by linear eddy diffusion model. Add a small disturbance which is induced by velocity gradient to Equation (42), we can obtain

\nabla\cdot\left[\left(\mathbf{u}+\delta\mathbf{u}\right)\left(k+\delta k\right)\right]=\left(P_{k}+\delta P_{k}\right)+\nabla\cdot\left[\Gamma_{k}\nabla\left(k+\delta k\right)\right]-\left(\varepsilon+\delta\varepsilon\right),

(43)

in which according to the definition of the finite differential of tensor. We can get $\delta\phi=\frac{d\phi}{d\left(\nabla\mathbf{u}\right)}:d\left(\nabla\mathbf{u}\right)$ , in which $\phi=\left[\mathbf{u},k,P_{k},\varepsilon\right]^{T}$ . Substitute Equation (42) into Equation (43). Then, volume integration of control volume is performed according to Equation (43), and omitting the second-order small quantities, we can get

\oint_{A}\mathbf{n}\cdot\left[\left(\mathbf{u}\delta k\right)-\left(\Gamma_{k}\nabla\delta k\right)\right]dA=\int_{V}\delta P_{k}dV-\int_{V}\delta\varepsilon dV.

(44)

An assumption is made here for the sake of the simplicity of the model, that

\int_{V}\delta P_{k}dV\approx\int_{V}\delta\varepsilon dV.

(45)

This assumption shows that the net flux of $\delta k$ , which is induced by the velocity gradient $\nabla\mathbf{u}$ , is nearly zero. It means that no matter how the $\nabla\mathbf{u}$ changes, $\delta k$ changes little. Comparing with $\delta\mathbf{R}^{V}=\frac{d\mathbf{R}^{V}}{d\left(\nabla\mathbf{u}\right)}:d\left(\nabla\mathbf{u}\right)$ , initial fields and boundary conditions are both implemented with $\delta k\approx 0$ and $\delta\mathbf{R}^{V}\approx\mathbf{0}$ . With the same $d\left(\nabla\mathbf{u}\right)$ , the change of $\delta k$ is smaller than that of $\delta\mathbf{R}^{V}$ . Therefore, as a zero-order approximation, $dk/d\left(\nabla\mathbf{u}\right)$ will be ignored.

In fact, in the original Boussinesq hypothesis [7], the constitutive relation should be written in

-\overline{u^{\prime}_{i}u^{\prime}_{j}}=-\frac{1}{\rho}p_{t}\delta_{ij}+\zeta_{ijkl}\frac{\partial u_{l}}{\partial x_{k}},

(46)

where $p_{t}$ is the pressure that is induced by fluctuating velocity, and it is corresponding to the pressure $p$ in Equation (2). Then, $p_{t}=\frac{2}{3}\rho k$ is adopted to the following EVM methods [10, 13, 25]. Physically, the effect of velocity gradient on $p_{t}=\frac{2}{3}\rho k$ should be similar to that on $p$ . According to the Helmholtz velocity decomposing theorem, the velocity gradient $\nabla\mathbf{u}$ can be decomposed into the rate of strain tensor $\mathbf{S}=\frac{1}{2}\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)$ for elastic deformation and the rotation tensor $\mathbf{W}=\frac{1}{2}\left(\nabla\mathbf{u}-\nabla\mathbf{u}^{T}\right)$ for rigid rotation. There is no direct and explicit relation between the isotropic pressure and the velocity gradient. Therefore, in the boundary layer theory, the pressure in the fully developed region of turbulence will directly act on the near-wall region under the first-order approximation [29]. The assumption in Equation (45) is derived from this physical consideration. For turbulent kinetic energy dissipation rate $\varepsilon$ , which is only directly acting on $k$ or $p_{t}$ , the same assumption and consideration are adopted.

After the above derivation, the tensorial derivative of $P_{ij}^{M}$ to $\nabla\mathbf{u}$ can be obtained

	$\displaystyle{\mathbf{P}^{M}}^{\prime}=$	$\displaystyle-\int_{V}\left[\left(2\delta_{ij}R_{\underline{i}m}^{V}\delta_{\underline{i}l}\delta_{mk}+2\delta_{ij}\zeta_{\underline{i}mkl}\frac{\partial u_{\underline{i}}}{\partial x_{m}}\right)+\left(\epsilon_{ij}R_{\underline{i}\underline{i}}^{V}\Omega_{kl}+\epsilon_{ij}\omega_{Z}\zeta_{\underline{i}\underline{i}kl}\right)\right]dV$		(47)
		$\displaystyle+\int_{V}\frac{2}{3}k\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)dV.$		(47)

From Equation (31), it is clearly seen that the tensorial derivative of the remaining terms to $\nabla\mathbf{u}$ can be referred to the procedure and operation in Equation (47). Therefore, the tensorial derivative of the remaining terms can refer in Appendix A.

Based on the transport equation of $\mathbf{R}^{V}$ as Equation (47), and tensorial derivative to $\nabla\mathbf{u}$ as Equation (33) to (47) and Appendix A. A transport equation of $\zeta_{ijkl}$ is obtained

\frac{\partial\zeta_{ijkl}}{\partial t}+u_{m}\frac{\partial\zeta_{ijkl}}{\partial x_{m}}=P_{ijkl}^{Strain}+P_{ijkl}^{Vortivity}+D_{ijkl}^{\zeta}+S_{ijkl}^{Positive}+E_{ijkl}^{Dissipation},

(48)

in which,

	$P_{ijkl}^{Strain}=\left(C_{2}-1\right)\left(2\delta_{ij}R_{\underline{i}k}^{V}\delta_{\underline{i}l}+2\delta_{ij}\zeta_{\underline{i}mkl}\frac{\partial u_{\underline{i}}}{\partial x_{m}}\right)+\frac{2}{3}\left(1-C_{2}\right)\delta_{ij}\left(\zeta_{rskl}\frac{\partial u_{r}}{\partial x_{s}}+R_{kl}^{V}\right),$		(49a)
	$P_{ijkl}^{Vorticity}=\left(C_{2}-1\right)\left(\epsilon_{ij}R_{\underline{i}\underline{i}}^{V}\Omega_{kl}+\epsilon_{ij}\omega_{Z}\zeta_{\underline{i}\underline{i}kl}\right),$		(49b)
	$D_{ijkl}^{\zeta}=\frac{\partial}{\partial x_{m}}\left[\left(C_{S}\frac{k^{2}}{\varepsilon}+\nu\right)\frac{\partial\zeta_{ijkl}}{\partial x_{m}}\right],$		(49c)
	$S_{ijkl}^{Positive}=\frac{2}{3}\left(1-C_{2}\right)k\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right),$		(49d)
	$E_{ijkl}^{Dissipation}=-C_{R}\frac{\varepsilon}{k}\zeta_{ijkl}.$		(49e)

The tensor base of Equation (48) is $\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}$ , and the dimension is $m^{2}/s^{2}$ .

The transport model of $\zeta_{ijkl}$ contains five source terms. $P_{ijkl}^{Strain}$ is the strain production term, which indicates that $\zeta_{ijkl}$ is generated by the normal elastic deformation. $P_{ijkl}^{Vorticity}$ is the vorticity production term, that $\zeta_{ijkl}$ is generated by the change of time-averaging vorticity, which is used to generate small-scale turbulence. $D_{ijkl}^{\zeta}$ is the diffusion term, which contains molecular diffusion and turbulence-induced diffusion. $\zeta_{ijkl}$ diffuses from the high turbulent intensity region to the low turbulent intensity region, which makes the turbulence isotropic. Source terms contain a positive-definite source $S_{ijkl}^{Positive}$ and a dissipation term $E_{ijkl}^{Dissipation}$ . From the research of Gama et al. and Wirth et al.[2, 3], we can know that $\zeta_{ijkl}$ can be positive or negative. When $\zeta_{ijkl}$ is positive, from constitutive Equation (10), it is obvious that the direction of the component of velocity gradient $\frac{\partial u_{l}}{\partial x_{k}}$ is opposite to Reynolds stress, which means the momentum transportation by fluctuating velocity is from the high region to the low region. In this situation, $S_{ijkl}^{Positive}$ performs as a source, and $E_{ijkl}^{Dissipation}<0$ performs as a dissipation. Therefore, for positive $\zeta_{ijkl}$ , flow fields with high turbulent kinetic energy $k$ will increase $\zeta_{ijkl}$ . At the same time, flow fields with high turbulent specific kinetic energy dissipation rate $\omega\varpropto\varepsilon/k$ will decrease $\zeta_{ijkl}$ . When $\zeta_{ijkl}$ is negative, from constitutive Equation (10), we can know that the direction of the component of velocity gradient $\frac{\partial u_{l}}{\partial x_{k}}$ is the same to $\mathbf{R}^{V}$ , which means that the momentum transportation by fluctuating velocity is from the low region to the high region. This is not obeying Fick’s law [30], which is called by CGT (Converse Gradient Transport) phenomenon [31, 32, 33]. CGT is first observed by experiments in fluid mechanics [34] and is characterized by the negative turbulent kinetic energy generation term. It is obvious that linear EVM constitutive relation cannot characterize the CGT, because the eddy viscosity $\nu_{t}\varpropto k^{2}/\varepsilon$ is always positive. When $\zeta_{ijkl}$ is negative, $S_{ijkl}^{Positive}$ is still positive, and acts as a sink; $E_{ijkl}^{Dissipation}>0$ also acts as a dissipation to prevent the $\zeta_{ijkl}$ more negative. This means the CGT of momentum transportation by fluctuating velocity, which is caused by $\zeta_{ijkl}<0$ , will attenuate faster than the general transport phenomenon. Equation (48) shows a preliminary transport model of high-order eddy viscosity tensor $\zeta_{ijkl}$ from modeled second-order moment transport equation with tensorial analysis.

From the perspective of properties of Reynolds stress, Equation (48) can consider the anisotropy of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ which compares to Spalart-Allmaras one-equation model. Also, this model can consider the relaxation and history effect of $\overline{u^{\prime}_{i}u^{\prime}_{j}}$ which compares to non-linear EVM. In the following sections, we will make the numerical analysis and model validation to investigate the characteristics of Equation (48) in the complex two-dimensional turbulent flow.

3 Numerical analysis and model validation

3.1 Numerical methodology and approximation

As discussed in the previous sections, Equation (48) is the high-order eddy viscosity tensor transport model to describe the two-dimensional turbulent flow. Since the model cannot compute the Reynolds stress directly like the second-order moment model. Therefore, $k$ equation and $\varepsilon$ equation still need to be supplemented to close the turbulence model. Contracting the index $i$ and $j$ of the second-order moment transport model will get the transport equation of turbulent kinetic energy. In Section 2.2, Equation (30) is modeled with the diffusion term of $k$ equation. Hence the modeled $k$ equation will be re-written here

\frac{\partial k}{\partial t}+u_{j}\frac{\partial k}{\partial x_{j}}=-\overline{u^{\prime}_{i}u^{\prime}_{j}}\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left[\left(C_{k}\frac{k^{2}}{\varepsilon}+\nu\right)\frac{\partial k}{\partial x_{j}}\right]-\varepsilon.

(50)

For ordinary fluid (Prandtl number near 1), $\varepsilon$ equation can be modeled [35, 19] as

\frac{\partial\varepsilon}{\partial t}+u_{j}\frac{\partial\varepsilon}{\partial x_{j}}=-C_{\varepsilon 1}\frac{\varepsilon}{k}\overline{u^{\prime}_{i}u^{\prime}_{j}}\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left[\left(C_{\varepsilon}\frac{k^{2}}{\varepsilon}+\nu\right)\frac{\partial\varepsilon}{\partial x_{j}}\right]-C_{\varepsilon 2}\frac{\varepsilon^{2}}{k},

(51)

in which $C_{\varepsilon}=0.07,C_{\varepsilon 1}=1.44$ , and $C_{\varepsilon 2}=1.92$ are empirical constants. Thus the turbulence model is closed and can be solved by Equations (6), (10), (48), (50), and (51). In this article, a turbulence model for two-dimensional turbulent flow in 2D-3C turbulence structure which named as $k-\varepsilon-\zeta$ is developed.

The finite volume method (FVM) is adopted in this article, and the model is implemented in OpenFOAM. OpenFOAM adopts the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm [36] to solve the coupling of pressure and velocity during the calculation process. The SIMPLE algorithm solves the inner and outer iterations, the inner iteration solves the large sparse matrix, and the outer iteration updates the coefficients of the equation with the results of the inner iteration. The SIMPLE algorithm does not consider the influence of the adjacent grid velocity correction on the current grid velocity in the discretization. Hence, the field relaxation factors need to be introduced to ensure the stability of the iterative process. However, $k-\varepsilon-\zeta$ model, which is developed in this article, needs to solve 17 equations for two-dimensional turbulent flow in OpenFOAM. The variants contain $u_{1}$ , $u_{2}$ , $p$ , $k$ , $\varepsilon$ , and 12 independent components of $\zeta_{ijkl}$ . The second-order moment model directly performs the transport equations for Reynolds stress, which can characterize all properties of Reynolds stress, but its numerical convergence and computational cost have been criticized [10, 13]. Nonetheless, the second-order moment model only needs to solve 8 equations in a two-dimensional turbulent flow. The variants contain $u_{1}$ , $u_{2}$ , $p$ , $\varepsilon$ , $\overline{u^{\prime}_{1}u^{\prime}_{1}}$ , $\overline{u^{\prime}_{1}u^{\prime}_{2}}$ , $\overline{u^{\prime}_{2}u^{\prime}_{2}}$ , $\overline{u^{\prime}_{3}u^{\prime}_{3}}$ . It is obvious that $k-\varepsilon-\zeta$ model even needs to solve 9 more equations than the second-order moment model. Therefore, how to ensure convergence and stability is an urgent problem when solving this large group of strongly coupled nonlinear equations of $k-\varepsilon-\zeta$ model. Through the numerical experiments on OpenFOAM, we find that only adjusting the field relaxation factor and matrix relaxation factor is hard to guarantee the stability of the solution. Hence an approximation is made for the production term in turbulent kinetic energy transport equation

	$P_{k}=-\overline{u^{\prime}_{i}u^{\prime}_{j}}\frac{\partial u_{i}}{\partial x_{j}}=\left(\zeta_{ijkl}\frac{\partial u_{l}}{\partial x_{k}}\right)\frac{\partial u_{i}}{\partial x_{j}},$		(52a)
	$P_{k}^{apx}=-\overline{u^{\prime}_{i}u^{\prime}_{j}}\frac{\partial u_{i}}{\partial x_{j}}=\nu_{t}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\frac{\partial u_{i}}{\partial x_{j}}.$		(52b)

$\nu_{t}=C_{\mu}\frac{k^{2}}{\varepsilon}$ is the eddy viscosity which is calculated in the form of linear EVM. Equation (52a) is the strict form of $k-\varepsilon-\zeta$ model in this article, and Equation (52b) is the approximation form used in numerical simulation with OpenFAOM to ensure the stability and convergence of the numerical solution. It is obvious that $\nu_{t}>0$ in $P_{k}^{apx}$ which cannot characterize the CGT phenomena. From the analysis in Section 2.2, we can know that $\zeta_{ijkl}$ can be positive or negative. Hence the strict form $P_{k}$ can reveal the CGT phenomena. For the preliminary research of this model, however, $P_{k}^{apx}$ will be adopted. When the goal is to reveal the CGT phenomenon with this model, a more robust numerical algorithm and the strategy of solution will be adopted to $P_{k}$ , which is not the focus of this article.

3.2 Numerical simulation and analyzation

The first case is a two-dimensional straight channel turbulence, and Re is 40000. We calculate the turbulent flow under fully-developed condition with $k-\varepsilon-\zeta$ model and standard $k-\varepsilon$ model [8], For two-dimensional channel turbulence, we focus on the properties of its dimensionless Reynolds stress characteristic.

Lumley [4] investigated the eigenvalues and characteristics of dimensionless Reynolds stress anisotropy tensor $b_{ij}=\overline{u^{\prime}_{i}u^{\prime}_{j}}/\overline{u^{\prime}_{k}u^{\prime}_{k}}-\frac{1}{3}\delta_{ij}$ . According to tensorial analysis, a symmetric second-order tensor can be diagonalized, and the eigenvalues are real numbers. The characteristic polynomial can be obtained

|\mathbf{b}-\lambda\mathbf{I}|=\lambda^{3}-\Psi_{1}^{\mathbf{b}}\lambda^{2}+\Psi_{2}^{\mathbf{b}}\lambda-\Psi_{3}^{\mathbf{b}}=0,

(53)

where $\lambda$ is the eigenvalue of $b_{ij}$ , in continuum mechanics, $\lambda$ also represents the principal stress. $\Psi_{1}^{\mathbf{b}}$ , $\Psi_{2}^{\mathbf{b}}$ , and $\Psi_{3}^{\mathbf{b}}$ are first, second, and third primary invariants of $b_{ij}$ , respectively.

	$\Psi_{1}^{\mathbf{b}}=tr\left(\mathbf{b}\right)=b_{ii}=0,$		(54a)
	$\Psi_{2}^{\mathbf{b}}=\frac{1}{2}\left(b_{ii}b_{jj}-b_{ij}b_{ij}\right)=-\frac{1}{2}b_{ij}b_{ij},$		(54b)
	$\Psi_{3}^{\mathbf{b}}=\|\mathbf{b}\|.$		(54c)

After the diagonalization of $b_{ij}$ , we can get the dimensionless Reynolds stress anisotropy tensor in principle axes, which only contain the normal stress

\widetilde{b_{ij}}=\begin{pmatrix}\lambda_{1}&&\\ &\lambda_{2}&\\ &&-\left(\lambda_{1}+\lambda_{2}\right)\end{pmatrix}.

(55)

Lumley [4] defined the parameter $\eta$ and $\xi$ , according to the properties of $\widetilde{b_{ij}}$ , we can derive the following relations

	$6\eta^{2}=-2I_{2}=b_{ij}b_{ji}\Rightarrow\eta^{2}=\frac{1}{3}\left(\lambda_{1}^{2}+\lambda_{1}\lambda_{2}+\lambda_{2}^{2}\right),$		(56a)
	$6\xi^{3}=3I_{3}=b_{ij}b_{jk}b_{ki}\Rightarrow\xi^{3}=-\frac{1}{2}\lambda_{1}\lambda_{2}\left(\lambda_{1}+\lambda_{2}\right).$		(56b)

By considering a series of simplified turbulence conditions with Equation (56a) and (56b), the Lumley triangle can be obtained. In other words, if the characteristic quantities of dimensionless Reynolds stress anisotropy tensor could not be wrapped by the Lumley triangle, it means that this turbulence model is not realizable. For $k-\varepsilon-\zeta$ model, characteristic quantities can be derived as

	$6\eta^{2}=b_{ij}b_{ji}=\frac{\left(\zeta_{ijkl}\frac{\partial u_{l}}{\partial x_{k}}\right)\cdot\left(\zeta_{jikl}\frac{\partial u_{l}}{\partial x_{k}}\right)}{4k^{2}},$		(57a)
	$6\xi^{3}=b_{ij}b_{jk}b_{ki}=-\frac{\left(\zeta_{ijmn}\frac{\partial u_{n}}{\partial x_{m}}\right)\cdot\left(\zeta_{jkmn}\frac{\partial u_{n}}{\partial x_{m}}\right)\cdot\left(\zeta_{kimn}\frac{\partial u_{n}}{\partial x_{m}}\right)}{8k^{3}}.$		(57b)

Thus, the results of the dimensionless Reynolds stress anisotropy tensor for $k-\varepsilon-\zeta$ model and standard $k-\varepsilon$ model is obtained, which are shown in Figure 1. For characteristic quantities $\eta$ and $\xi$ of $k-\varepsilon-\zeta$ model and standard $k-\varepsilon$ model are both in the Lumley triangle. It means that $k-\varepsilon-\zeta$ model and standard $k-\varepsilon$ model are both realizable. However, for turbulent flow in the two-dimensional straight channel, the trajectory of the characteristic quantities of $k-\varepsilon-\zeta$ model is along the $\eta=\xi$ , which is in good agreement with DNS channel flow [37].

For a fully-developed straight channel turbulence, it is obvious that $\frac{\partial u_{1}}{\partial x_{1}}=0$ and $\frac{\partial u_{2}}{\partial x_{2}}=0$ . Therefore, standard $k-\varepsilon$ model, which adopts the linear EVM constitutive, cannot capture the differences between $\overline{u^{\prime}_{1}u^{\prime}_{1}}$ and $\overline{u^{\prime}_{2}u^{\prime}_{2}}$ because linear EVM indicates that the turbulence is isotrpic, and the trajectory will pass through $\left(0,0\right)$ in Lumley triangle. Hence, from Figure 1 we can know that $k-\varepsilon-\zeta$ model is superior to the standard $k-\varepsilon$ model in describing the characteristics of turbulence.

Refer to caption — Figure 1: Characteristic quantities distribution of $k-\varepsilon-\zeta$ model and standard $k-\varepsilon$ model in Lumley triangle. Red solid squares represent the turbulent mixing layer obtained by Bell and Mehta [38]. Red solid dots represent the DNS channel flow obtained by Kim et al. [37]. Dark blue star-type nodes represent the results of $k-\varepsilon-\zeta$ model. Light blue fork-type nodes represent the results of standard $k-\varepsilon$ model.

By solving the transport equation of $\zeta_{ijkl}$ , we also can get the distribution of each component of high-order eddy viscosity tensor. Figure 2 shows the numerical results of 12 independent components of high-order eddy viscosity tensor. From a numerical point of view, we can find that $\zeta_{1111}$ and $\zeta_{2222}$ play a major role in Reynolds viscous stress tensor, which $\zeta_{1111}$ acts on $R_{11}^{V}$ , and $\zeta_{2222}$ acts on $R_{22}^{V}$ .

We also can see that 12 independent components can be divided into two parts, the axisymmetric part, and the centrosymmetric part about the center of the straight channel. The main components of Reynolds viscous stress can be written in

	$R_{11}^{V}=\underbrace{\zeta_{1111}\frac{\partial u_{1}}{\partial x_{1}}+\zeta_{1122}\frac{\partial u_{2}}{\partial x_{2}}}_{axisymmetric}+\underbrace{\zeta_{1112}\frac{\partial u_{2}}{\partial x_{1}}+\zeta_{1121}\frac{\partial u_{1}}{\partial x_{2}}}_{centrosymmetric},$		(58a)
	$R_{12}^{V}=\underbrace{\zeta_{1212}\frac{\partial u_{2}}{\partial x_{1}}+\zeta_{1221}\frac{\partial u_{1}}{\partial x_{2}}}_{axisymmetric}+\underbrace{\zeta_{1211}\frac{\partial u_{1}}{\partial x_{1}}+\zeta_{1222}\frac{\partial u_{2}}{\partial x_{2}}}_{centrosymmetric},$		(58b)
	$R_{22}^{V}=\underbrace{\zeta_{2211}\frac{\partial u_{1}}{\partial x_{1}}+\zeta_{2222}\frac{\partial u_{2}}{\partial x_{2}}}_{axisymmetric}+\underbrace{\zeta_{2212}\frac{\partial u_{2}}{\partial x_{1}}+\zeta_{2221}\frac{\partial u_{2}}{\partial x_{1}}}_{centrosymmetric},$		(58c)

As mentioned above, for two-dimensional fully-developed turbulence in the straight channel, $\frac{\partial u_{1}}{\partial x_{1}}=0$ and $\frac{\partial u_{2}}{\partial x_{2}}=0$ . Therefore, axisymmetric parts of $R_{11}^{V}$ and $R_{22}^{V}$ are nearly zero, and centrosymmetric parts reflect the anisotropy of Reynolds stress. For Reynolds viscous shear stress $R_{12}^{V}$ , the main part is the axisymmetric part, and the centrosymmetric part is nearly zero. It can be seen from the value of each component that, for Reynolds viscous stress, the axisymmetric part is the main part, and the centrosymmetric part is more likely to play a regulatory role.

For a complex two-dimensional turbulent flow, the flow in an asymmetric planar diffuser will be adopted in this article to make the numerical simulation of different RANS models. Turbulent flow in the asymmetric planar diffuser is a benchmark for Large Eddy Simulation (LES) because of the presence of an adverse pressure gradient, and the formation of an unsteady separation recirculation in the expansion section of the diffuser, which may carry out many challenges to LES prediction [39]. Many simulations based on LES are carried out for the asymmetric planar diffuser [40, 39, 41]. The turbulent flow in the diffuser contains the behaviors of adverse pressure gradient of flow and anisotropy and relaxation of Reynolds stress, which are also critical to the validation of RANS models. Obi et al. [6] made experimental and numerical researches on the diffuser in a two-dimensional turbulent flow. The LDV measurement of the turbulent separating flow in an asymmetric diffuser was used to capture the behavior of the turbulent flow. Obi et al. [6] demonstrated that standard $k-\varepsilon$ model fails to predict the separation of the flow. In this article, $k-\varepsilon-\zeta$ model will make the simulation on two-dimensional turbulent flow in the diffuser, investigation of the transportation of high-order eddy viscosity tensor will also be carried out.

Figure 3 shows the computational domain and mesh of the expansion part of the diffuser, and the geometric dimensions are referred to the Obi et al. [6] experiments. The height of the inlet channel is $H=2cm$ , and the length of the inlet channel is $100H$ to ensure the fully-developed two-dimensional turbulent channel flow at $Re=2.0\times 10^{4}$ . Obi et al. [6] also examined the two-dimensionality of the flow field, and the difference in mean velocity profile in the spanwise direction was verified to be less than 5% over 90% and 60% of the channel span at the inlet and outlet plane. In this article, the inlet velocity $U_{in}$ is $2m/s$ , and the kinematic viscosity $\nu$ is $10^{-6}$ to ensure the channel flow at $Re=2.0\times 10^{4}$ .

Figure 4 shows the simulation results which are calculated by OpenFOAM for $k-\varepsilon-\zeta$ model and standard $k-\varepsilon$ model. From the results, it is obvious that $k-\varepsilon-\zeta$ model can predict the separation bubble, and standard $k-\varepsilon$ model fails to predict the separation bubble. These phenomena are in the agreement with Obi’s numerical results. Therefore, $k-\varepsilon-\zeta$ model can capture the anisotropy and relaxation of Reynolds stress, while standard $k-\varepsilon$ model cannot. The experimentally obtained flow detachment and reattachment points are $11H$ and $26H$ , respectively [6]. For $k-\varepsilon-\zeta$ model, the calculation results obtain the flow detachment and reattachment points are nearly $6.16H$ and $27H$ , respectively. The prediction of the separation bubble is in good agreement with the experimental measurements. Kaltenbach et al. [39] carried out the three-dimensional LES simulation for the planar asymmetric diffuser, the detachment and reattachment points are also investigated. Crawford and Birk [42] calculated the asymmetric planar diffuser based on the $v^{2}-f$ turbulence model [43] of ANSYS Fluent to study the flow characteristics in two-dimensional turbulent flow, which also contained the detachment and reattachment points in the expansion section of the diffuser. All the results are summarized in Table 1.

\tbl

Comparison of Detachment and Reattachment points and bubble length between different turbulence models and experimental data. Turbulence model or experimental data Detachment Reattachment Bubble length Obi Experiment [6] $11H$ $26H$ $15H$ $k-\varepsilon-\zeta$ model $6.16$ $27$ $20.84$ LES [39] $6H$ $27.5H$ $21.5H$ $v^{2}-f$ [42] $6.15H$ $27.6H$ $21.45H$ \tabnoteNote: All parameters are in the width of the inlet channel $H=0.02m$

To better analyze the turbulent flow characteristics, we carry out quantitative analysis, which compares with the results of Obi et al. [6] experimental measurements. Figure 5 shows the results of dimensionless streamwise mean velocity and dimensionless Reynolds shear stress at $13.2H$ and $19.2H$ for $k-\varepsilon-\zeta$ model, standard $k-\varepsilon$ model, and LRR model, respectively. From Figure 3, we can know that $13.2H$ (26.4 cm) and $19.2H$ (38.4 cm) are both at the expansion section of the asymmetric planar diffuser, which is a strict region for turbulent flow. Obi et al. [6] also carried out a numerical simulation using LRR model, hence we also add a simulation with LRR model to calculate the two-dimensional planar asymmetric diffuser in this article, to make the comparison and analysis with $k-\varepsilon-\zeta$ model and standard $k-\varepsilon$ model.

From the calculation results, neither LRR model nor standard $k-\varepsilon$ model can simulate the separation bubble in the expansion section. $k-\varepsilon-\zeta$ model can simulate the separation bubble, and the results of streamwise mean velocity at $13.2H$ and $19.2H$ are in good agreement with the experimental measurements. At $13.2H$ in the expansion section, the prediction accuracy of the streamwise mean velocity of $k-\varepsilon-\zeta$ model is increased by nearly 300% when compares to standard $k-\varepsilon$ model, and nearly 200% when compares to LRR model. And at $19.2H$ in the expansion section, the prediction accuracy of the streamwise mean velocity of $k-\varepsilon-\zeta$ model is increased by nearly 200% and 100% when compared with standard $k-\varepsilon$ model and LRR model, respectively. However, the streamwise mean velocity of model is higher than the experimental measurement at the upper wall of the asymmetric planar diffuser, because of the numerical approximation of the production term of turbulent kinetic energy in Equation (52b). For the lower wall, the velocity gradient near the wall is not large due to the low velocity in the separation bubble region, hence $k-\varepsilon-\zeta$ model can simulate the streamwise mean velocity well. When near the upper wall region, the velocity is much larger than the separation bubble region, hence the velocity gradient is much larger than the separation bubble region, the effect of the numerical approximation in Equation (52b) is obvious and the numerical deviation appears.

For Reynolds shear stress at $13.2H$ and $19.2H$ in the expansion section, we can find that standard $k-\varepsilon$ model always underestimates the Reynolds shear stress, and $k-\varepsilon-\zeta$ model always overestimates the Reynolds shear stress, which are shown in Figure 5.(b) and (d). Only the LRR model accurately estimates the Reynolds shear stress at $19.2H$ . The main reason for the overestimation of $k-\varepsilon-\zeta$ model is the numerical approximation in Equation (52b), which will overestimate the turbulent kinetic energy and Reynolds shear stress. However, from the overall calculation results, we can see that the $k-\varepsilon-\zeta$ model can preliminarily reveal the turbulent flow characteristics of the two-dimensional asymmetric planar diffuser.

By solving the conserving Equation (48), high-order eddy viscosity tensor $\zeta_{ijkl}$ can be calculated, which contains 12 independent components. Figure 6 shows the results of 12 independent components of high-order eddy viscosity tensor of complex turbulent flow in an asymmetric planar diffuser. The variation of the distribution of $\zeta_{ijkl}$ mainly starts near the flow separation point. And the positive and negative signs of the components of $\zeta_{ijkl}$ will change after the flow separation in the expansion section.

For an asymmetric planar diffuser, numerical analysis [42, 43] shows that the low-Re $v^{2}-f$ model, which considers the curvature correction and the near-wall flow characteristics, has a good performance and quite a high degree of accuracy. For low-Re models, such as $v^{2}-f$ model, $k-\omega$ SST model [44, 45], they capture more accurately the changes of geometric curvature and the changes of velocity gradient near the wall. Therefore, a large part of the reason why they capture the separated flow is that they capture that the $\tau_{\shortparallel}=\mu\frac{\partial u_{\shortparallel}}{\partial x_{\shortparallel}}=0$ near the wall due to the adverse pressure gradient is 0, in which $\tau_{\shortparallel}$ is the shear stress parallel to the wall, $\frac{\partial u_{\shortparallel}}{\partial x_{\shortparallel}}$ is the normal gradient of the velocity which is parallel to the wall.

However, for the high-Re models that need to use wall functions to approximate the near-wall flow characteristics, the capture of geometric curvature and the near-wall flow behavior is weaker than that of low-Re models. Due to $k-\varepsilon-\zeta$ model is derived from LRR model, hence $k-\varepsilon-\zeta$ model is also a high-Re model. It means that the capture of the separation bubble in the expansion region of the diffuser is based on the properties of Reynolds stress which are characterized by the turbulence model, rather than the analysis of flow characteristics near the wall like the low-Re model. The same phenomenon also can be observed when using nonlinear EVM to calculate the internal flow in a square duct, that the induced vortices due to the anisotropy property of Reynolds stress occur on the cross-section [46].

4 Discussion and conclusion

4.1 2D-3C turbulence structure and EV-RSM

In this article, the development of $k-\varepsilon-\zeta$ model is carried out under the 2D-3C turbulence structure. To explicitly obtain the transport model of high-order eddy viscosity tensor, however, analogy from the linear EVM model, an assumption for constitutive relation is made in Equation (15). It means that the Reynolds viscous normal stress on the z-axis is $R_{33}^{V}=0$ , and Reynolds normal stress on the z-axis is $R_{33}=\frac{2}{3}k$ . However, for the second-order moment model, when a two-dimensional turbulence flow in a 2D-3C turbulence structure is established, the Reynolds normal stress on the z-axis $R_{33}$ still needs to be solved by a transport equation. This will result in $R_{33}$ not equal to $\frac{2}{3}k$ . The transport equation of $R_{33}$ in 2D-3C turbulence structure is

$\displaystyle\frac{\partial\overline{u^{\prime}_{3}u^{\prime}_{3}}}{\partial t}+u_{1}\frac{\partial\overline{u^{\prime}_{3}u^{\prime}_{3}}}{\partial x_{1}}+u_{2}\frac{\partial\overline{u^{\prime}_{3}u^{\prime}_{3}}}{\partial x_{2}}=$	$\displaystyle-\frac{\partial\overline{u^{\prime}_{3}u^{\prime}_{3}u^{\prime}_{1}}}{\partial x_{1}}-\frac{\partial\overline{u^{\prime}_{3}u^{\prime}_{3}u^{\prime}_{2}}}{\partial x_{2}}$	(59)
	$\displaystyle-\frac{\partial}{\partial x_{1}}\left(\nu\frac{\partial\overline{u^{\prime}_{3}u^{\prime}_{3}}}{\partial x_{1}}\right)-\frac{\partial}{\partial x_{2}}\left(\nu\frac{\partial\overline{u^{\prime}_{3}u^{\prime}_{3}}}{\partial x_{2}}\right)$
	$\displaystyle-2\nu\overline{\frac{\partial u^{\prime}_{3}}{\partial x_{1}}\cdot\frac{\partial u^{\prime}_{3}}{\partial x_{1}}+\frac{\partial u^{\prime}_{3}}{\partial x_{2}}\cdot\frac{\partial u^{\prime}_{3}}{\partial x_{2}}}.$

This may contain the incompatibility between RSM (Reynolds Stress Model) and EV-RSM (Eddy Viscosity – Reynolds Stress Model). EV-RSM is a RANS turbulence model, which uses the Boussinesq hypothesis to maintain the simplicity and considers incorporating the second-order moment transport equation to improve the accuracy in describing the properties of Reynolds stress. It is obvious that $k-\varepsilon-\zeta$ model is an EV-RSM model, because this model adopts the Boussinesq hypothesis and is derived from the second-order moment model. Spalart-Allmaras one-equation model is also an EV-RSM model, the transport equation of eddy viscosity $\nu_{t}$ has a strict form, which can be built with the second-order moment model. In Spalart-Allmaras one-equation model, the constitutive relation is the linear EVM hypothesis, which is shown in Equation (7). For incompressible turbulent flow, square both sides of Equation (7), and simplify with the continuity equation, we can get

\overline{u^{\prime}_{i}u^{\prime}_{j}}\cdot\overline{u^{\prime}_{i}u^{\prime}_{j}}=\frac{4}{3}k^{2}+4\nu_{t}^{2}S_{kl}S_{kl}.

(60)

Make the material differentiation of Equation (60), we can get

\frac{D}{Dt}\left(\overline{u^{\prime}_{i}u^{\prime}_{j}}\cdot\overline{u^{\prime}_{i}u^{\prime}_{j}}\right)=\frac{8}{3}k\frac{Dk}{Dt}+4\nu_{t}^{2}\frac{D}{Dt}\left(S_{kl}S_{kl}\right)+8S_{kl}S_{kl}\nu_{t}\frac{D\nu_{t}}{Dt}.

(61)

Substitute the constitutive relation of linear EVM in Equation (7), and make the simplification with the continuity equation, we can get

-\frac{D}{Dt}\left(2\nu_{t}S_{ij}\overline{u^{\prime}_{i}u^{\prime}_{j}}\right)=4\nu_{t}^{2}\frac{D}{Dt}\left(S_{kl}S_{kl}\right)+8S_{kl}S_{kl}\nu_{t}\frac{D\nu_{t}}{Dt}.

(62)

With the mathematical transformation, the left-hand side term can be transformed into

-\frac{D}{Dt}\left(2\nu_{t}S_{ij}\overline{u^{\prime}_{i}u^{\prime}_{j}}\right)=-\overline{u^{\prime}_{i}u^{\prime}_{j}}\frac{D}{Dt}\left(2\nu_{t}S_{ij}\right)-2\nu_{t}S_{ij}\frac{D\overline{u^{\prime}_{i}u^{\prime}_{j}}}{Dt}.

(63)

Substitute the constitutive relation of linear EVM in Equation (7), we can get

$\displaystyle-\frac{D}{Dt}\left(2\nu_{t}S_{ij}\overline{u^{\prime}_{i}u^{\prime}_{j}}\right)$	$\displaystyle=2\nu_{t}S_{ij}\frac{D}{Dt}\left(2\nu_{t}S_{ij}\right)-2\nu_{t}S_{ij}\frac{D\overline{u^{\prime}_{i}u^{\prime}_{j}}}{Dt}$	(64)
	$\displaystyle=-2\nu_{t}S_{ij}\frac{D\overline{u^{\prime}_{i}u^{\prime}_{j}}}{Dt}-2\nu_{t}S_{ij}\frac{D\overline{u^{\prime}_{i}u^{\prime}_{j}}}{Dt}$
	$\displaystyle=-4\nu_{t}S_{ij}\frac{D\overline{u^{\prime}_{i}u^{\prime}_{j}}}{Dt}.$

Hence Equation (62) can be converted to

-4\nu_{t}S_{ij}\frac{D\overline{u^{\prime}_{i}u^{\prime}_{j}}}{Dt}=4\nu_{t}^{2}\frac{D}{Dt}\left(S_{kl}S_{kl}\right)+8S_{kl}S_{kl}\nu_{t}\frac{D\nu_{t}}{Dt}.

(65)

We can find that $8S_{kl}S_{kl}\nu_{t}$ is a scalar, which is the reason why we square the linear EVM constitutive relation at first. Divide both sides of the Equation (65) by $8S_{kl}S_{kl}\nu_{t}$

\frac{D\nu_{t}}{Dt}=-\frac{S_{ij}}{2S_{kl}S_{kl}}\frac{D\overline{u^{\prime}_{i}u^{\prime}_{j}}}{Dt}-\frac{\nu_{t}}{2S_{kl}S_{kl}}\frac{D}{Dt}\left(S_{ij}S_{ij}\right).

(66)

Hence, we can get the transport model of eddy viscosity $\nu_{t}$ in the strict form

	$\displaystyle\frac{\partial\nu_{t}}{\partial t}+u_{m}\frac{\partial\nu_{t}}{\partial x_{m}}=$	$\displaystyle-\frac{S_{ij}}{2S_{kl}S_{kl}}\underbrace{\left(\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial t}+u_{m}\frac{\partial\overline{u^{\prime}_{i}u^{\prime}_{j}}}{\partial x_{m}}\right)}_{Second-order\ moment\ model}$		(67)
		$\displaystyle-\frac{\nu_{t}}{2S_{kl}S_{kl}}\left[\frac{\partial\left(S_{ij}S_{ij}\right)}{\partial t}+u_{m}\frac{\partial\left(S_{ij}S_{ij}\right)}{\partial x_{m}}\right].$		(67)

It can be seen that the strict form of transport model of eddy viscosity $\nu_{t}$ contains the second-order moment transport model. Therefore, theoretically, Spalart-Allmaras one-equation model also has incompatibility when calculating the two-dimensional turbulent flow in 2D-3C turbulence structure. Because from the constitutive relation of Spalart-Allmaras one-equation, the Reynolds viscous normal stress on the z-axis is $R_{33}^{V}=0$ and the Reynolds normal stress on the z-axis is $R_{33}=\frac{2}{3}k$ , but the strict form transport model of eddy viscosity $\nu_{t}$ contains the second-order moment model, which needs to solve the transport Equation (59) of $\overline{u^{\prime}_{3}u^{\prime}_{3}}$ in two-dimensional turbulent flow.

It is also not a good choice to directly model and solve Equation (67). Because from Equation (67) we can know that if we want to get eddy viscosity $\nu_{t}$ , we need to solve the second-order moment equation first. If we have already solved the second-order moment equation, the Reynolds stress is obtained, we do not need to solve other equations or make another hypothesis further. In 1994, Spalart and Allmaras [20] empirically developed the one-equation model, which is motivated by algebraic models such as the Baldwin-Lomax model [47], and the Johnson-King model [48]. Spalart and Allmaras considered the one-equation model of four nested operating conditions, from the simplest condition which applicable only to free shear flows to the most complete condition which applicable to viscous flows past solid bodies and with laminar regions. For each new operating condition, some semi-theoretical and semi-empirical source terms would be added to the transport equation of eddy viscosity $\nu_{t}$ . By comparing with the experimental measurements, source terms and empirical constants in the transport model of eddy viscosity would be calibrated and updated. Therefore, Spalart-Allmaras one-equation model is obtained, which is successful in the field of aerodynamics with the small amount of calculation, stable convergence, and explanation of historical effect and relaxation of Reynolds stress. Based on the above analysis, we can find that the incompatibility between EV-RSM and the second-order moment model of two-dimensional turbulent flow under the 2D-3C turbulence structure can be reduced by introducing some empirical constants or source terms artificially depending on experimental measurements or DNS data calibration.

Motivated by Spalart-Allmaras one-equation model, an empirical damping scalar function $f\left(Re_{t},\mathbf{u},\nabla\mathbf{u},\cdots\right)$ could be added on $P_{k}^{apx}$ in Equation (52b), as

P_{k}^{apx}=f\left(Re_{t},\mathbf{u},\nabla\mathbf{u},\cdots\right)\cdot\nu_{t}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\frac{\partial u_{i}}{\partial x_{j}},

(68)

in which $Re_{t}=k^{2}/\left(\nu\varepsilon\right)$ is the turbulent Reynolds number. We think that the damping function is the function of $Re_{t}$ , $\mathbf{u}$ , $\nabla\mathbf{u}$ and other physical quantities. Thus, we can fit the empirical damping function by experimental measurements or DNS data. The fitting function of $f\left(Re_{t},\mathbf{u},\nabla\mathbf{u},\cdots\right)$ can be derived with machine learning or neural network in the future, which is not the focus of this article. Or to find a more stable numerical algorithm and iterative strategy to solve the large nonlinear equations of $k-\varepsilon-\zeta$ model without any numerical approximation and simplification. From an engineering perspective, the $k-\varepsilon-\zeta$ model developed in this article may not be practical, because this model contains a large number of nonlinear equations, and the computational complexity, convergence, and stability are urgent problems to be solved. However, $k-\varepsilon-\zeta$ model supplies a new method to study the high-order eddy viscosity tensor in two-dimensional turbulent flow with 2D-3C turbulence structure, which reveals the evolution law and some physical properties of high-order eddy viscosity tensor in the aspect of physics.

4.2 Conclusion and outlook

A framework of the transport model for high-order eddy viscosity tensor in 2D-3C turbulence structure has been developed which is named as $k-\varepsilon-\zeta$ model in this article. Starting from the constitutive relation of the Boussinesq hypothesis in RANS, by analogy with the constitutive of Navier-Stokes hypothesis, we obtain the constitutive relation in the form of high-order eddy viscosity tensor without simplification and approximation. Through the modeled second-order moment transport equation and turbulent kinetic energy transport equation, we can derive the transport equation of Reynolds viscous stress. Regarding the constitutive relation of high-order eddy viscosity tensor as a tensor function, we obtain the explicit transport equation of high-order eddy viscosity tensor through the derivative of the tensor function and mathematical transformation. The physical meaning of the source terms in the transport equation of high-order eddy viscosity tensor is also analyzed. For the transport equation of high-order eddy viscosity tensor, it contains the transient term, convection term, and source terms. We can figure out that high-order eddy viscosity tensor can be generated by the normal elastic deformation in the strain production term, and also can be generated by the change of time-averaging vorticity in the vorticity production term. The diffusion term also has been concluded in the transport equation to make the high-order eddy viscosity tensor isotropic. The dissipation term of high-order eddy viscosity tensor is obtained to characterize the dissipation behavior of the turbulent flow. In addition, a positive-definite source term also appears in the transport equation, when the components of high-order eddy viscosity tensor are positive, this source term acts as a source that will increase the components. When the components are negative, the positive-definite source acts as a sink which will decrease the components to prevent the CGT phenomena occurs.

The preliminary validation by numerical simulation with OpenFOAM is also carried out in this article. Since the $k-\varepsilon-\zeta$ model contains 17 equations in the 2D-3C turbulence structure, which is more difficult to solve and make convergence than the second-order moment model in the aspect of solving nonlinear equations. Hence a numerical approximation of the production term in turbulent kinetic energy equation is made to enhance the stability and convergence during solving this model. We first analyze the characteristics of the $k-\varepsilon-\zeta$ model by simulating the turbulent flow in the two-dimensional straight channel. The results show that $k-\varepsilon-\zeta$ model can predict the anisotropy of Reynolds stress. The Lumley triangle has also been validated for $k-\varepsilon-\zeta$ model, and the eigenvalue distribution of the dimensionless Reynolds stress anisotropy tensor is in good agreement with DNS data. And then, we also calculate and verify the two-dimensional complex turbulent flow in a plane diffuser. The results show that $k-\varepsilon-\zeta$ model can effectively resolve the separation bubble. The detachment point and reattachment point are also quantitatively checked, which are in good agreement with the experimental measurements and LES calculation results. Compared with the standard $k-\varepsilon$ model, the $k-\varepsilon-\zeta$ model has a better analytical ability for complex turbulent flow.

Future research is needed to calibrate the damping function $f\left(Re_{t},\mathbf{u},\nabla\mathbf{u},\cdots\right)$ with experimental measurements or DNS datasets. The form and arguments of the damping function can be predicted by neural network [49, 50, 51] or machine learning [52, 53, 54, 55]. In addition, other iterative algorithms also can be studied to enhance the stability and convergence of during solving $k-\varepsilon-\zeta$ model. When the damping function is calibrated or other more stable iterative algorithms are developed, more severe and more detailed benchmark tests of the $k-\varepsilon-\zeta$ model should be carried out, to provide more useful information on high-order eddy viscosity tensor. These important developments and tests of $k-\varepsilon-\zeta$ model will be the subject of a future study.

Appendix A Tensor function derivative of source terms in Reynolds viscous stress equation.

In Section 2.2, we carry out the derivation of the tensor function derivative of $P_{ij}^{M}$ in the Reynolds viscous stress transport equation to the velocity gradient tensor. With the mathematical transformation and assumption in the derivation process, we can obtain the tensor function derivative of other source terms in the Reynolds viscous stress transport equation to the velocity gradient tensor, as

${\mathbf{D}^{M}}^{\prime}=\oint_{A}n_{m}\left(C_{S}\frac{k^{2}}{\varepsilon}\right)\frac{\partial\zeta_{ijkl}}{\partial x_{m}}dA,$		(69a)

$\displaystyle{\mathbf{\Phi}_{S}^{M}}^{\prime}$	$\displaystyle=\int_{V}C_{2}\left[\left(2\delta_{ij}R_{\underline{i}m}^{V}\delta_{\underline{i}l}\delta_{mk}+2\delta_{ij}\zeta_{\underline{i}mkl}\frac{\partial u_{\underline{i}}}{\partial x_{m}}\right)+\epsilon_{ij}R_{\underline{i}\underline{i}}^{V}\Omega_{kl}+\epsilon_{ij}\omega_{Z}\zeta_{\underline{i}\underline{i}kl}\right]dV$	(69b)
	$\displaystyle-\int_{V}C_{2}\frac{2}{3}k\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)dV$
	$\displaystyle-\int_{V}C_{2}\frac{2}{3}\delta_{ij}\left(\frac{\partial u_{r}}{\partial x_{s}}\zeta_{rskl}+R_{rs}^{V}\delta_{rl}\delta_{sk}\right)dV,$
${\mathbf{P}_{k}^{M}}^{\prime}=\int_{V}\frac{2}{3}\delta_{ij}\left(\frac{\partial u_{r}}{\partial x_{s}}\zeta_{rskl}+R_{rs}^{V}\delta_{rl}\delta_{sk}\right)dV,$		(69c)
${\mathbf{D}_{k}^{M}}^{\prime}=\mathbf{0}$		(69d)

All above derivatives are in the base of $\mathbf{e}_{i}\mathbf{e}_{j}\mathbf{e}_{k}\mathbf{e}_{l}$ .

Acknowledgement(s)

The authors are grateful to Prof. Nan Jiang at Tianjin University for the useful discussions and suggestions.

Funding

This work is supported by the National Key Research and Development Program of China (No. 2020YBF1902000).

References

[1] B. Dubrulle and U. Frisch. Eddy viscosity of parity-invariant flow. Physical Review A, 43:1–10,10, May 1991.
[2] B. Gama, M. Vergassola, and U. Frisch. Negative eddy viscosity in isotropically forced two-dimensional flow: linear and nonlinear dynamics. Journal of Fluid Mechanics, 260:95–126, 1994.
[3] A. Wirth, S. Gama, and U. Frisch. Eddy viscosity of three-dimensional flow. Journal of Fluid Mechanics, 188:249–264, 1995.
[4] J. L. Lumley. Computational modeling of turbulent flows. Number 18 in Physics. Acdemic Press, New York (NY), 1978.
[5] A. J. Simonsen and P. A. Krogstad. Turbulent stress invariant analysis: Clarification of existing terminology. Physics of Fluids, 17:1–4, 2005.
[6] S. Obi, K. Aoki, and S. Masuda. Experimental and computational study of turbulent separating flow in an asymmetric plane diffuser. Paper presented at: Ninth Symposium on ”Turbulent Shear Flows”, 1993 August 16–18; Kyoto, Japan.
[7] J. Boussinesq. Essai sur la theorie des eaux courantes. Impr. Nationale., 1877.
[8] B. E. Luander and D. B. Spalding. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering, 3:269–289, 1974.
[9] D. C. Wilcox. Reassessment of the scale-determining equation for advanced turbulence models. AIAA Journal, 26:1–12, 1988.
[10] K. Hanjalic and B. E. Launder. Reassessment of modeling turbulence via reynolds averaging: A review of second-moment transport strategy. Physics of Fluids, 33:1–20, 2021.
[11] T. B. Gatski and T. Jongen. Nonlinear eddy viscosity and algebraic stress models for solving complex turbulent flows. Progress in Aerospace Sciences, 36:655–682, 2000.
[12] T. S. Lund and E. A. Novikov. Parameterization of subgrid-scale stress by the velocity gradient tensor. Center for Trubulence Research Annual Research Briefs, 94:27–45, 1992.
[13] S. B. Pope. A more general effective-viscosity hypothesis. Journal of Fluid Mechanics, 72:331–340, 1975.
[14] T. H. Shih, J. Zhu, and J. L. Lumley. A new reynolds stress algebraic equation model. Computer Methods in Applied Mechanics and Engineering, 125:287–302, 1995.
[15] F. S. Lien, W. L. Chen, and M. A. Leschziner. Low-reynolds-number eddy-viscosity modelling based on non-linear stress-strain/vorticity relations. Engineering Turbulence Modelling and Experiments, 3:91–101, 1996.
[16] Z. X. Liu, B. J. Song, C. G. Gu, Y. M. Miao, and M. M. Su. A high-order anisotropic $k-\varepsilon$ model and its numerical investigation in two dimensional turbulent flow fields. Journal of Aerospace Power, 11:393–396, 1996.
[17] K. Hanjalic and B. E. Launder. 50 years of CFD in engineering sciences. Physics. Springer Nature Singapore Pte Ltd., Netherlands, 2020.
[18] Z. D. Wang and N. Jiang. On the modelling of eddy viscosity tensor for shear turbulence. Chinese Journal of Theoretical and Applied Mechanics, 27:137–141, 1995.
[19] B. E. Launder, G. J. Reece, and W. Rodi. Progress in the development of a reynolds-stress turbulence closure. Journal of Fluid Mechanics, 68:537–566, 1975.
[20] P. R. Spalart and S. R. Allmaras. A one-equation turbulence model for aerodynamic flows. La Recherche Aerospatiale, 1:5–21, 1994.
[21] G. Backus. A geometrical picture of anisotropic elastic tensors. Reviews of Geophysics and Space Physics, 8:633–671, 1970.
[22] T. Obermayer, C. Krempaszky, and E. Werner. Determination of the anisotropic elasticity tensor by mechanical spectroscopy. Continuum Mechanics and Thermodynamics, 34:165–184, 2022.
[23] S. Pope. Turbulent flows. Physics. Cambridge University Press, Cambridge, 2000.
[24] C. C. Shir. A preliminary numerical study of atmospheric turbulent flows in the idealized planetary boundary layer. Journal of the Atmospheric Sciences, 30:1327–1339, 1973.
[25] W. Q. Tao. Numerical heat transfer. Physics. Xi’an Jiaotong University Press, China, 2000.
[26] R. L. Bishop and S. I. Goldberg. Tensor analysis on manifolds. Mathematics. Dover, America, 1980.
[27] K. Z. Huang, M. D. Xue, and M. W. Lu. Tensor analysis. Mathematics. Tsinghua University Press, China, 2003.
[28] I. S. Sokolnikoff. Tensor analysis: theory and applications to geometry and mechanics of continua. Mathematics. Wiley, England, 1964.
[29] H. T. Schlichting and K. Gersten. Boundary layer theory. Physics. Springer, German, 1996.
[30] P. D. Davis, G. D. Parbrook, and G. N. C. Kenny. Basic physics and measurement in anaesthesia (Fourth edition). Physics. Butterworth-Heinemann, America, 1995.
[31] V. Dixit, L. Nuijens, and K. C. Helfer. Counter-gradient momentum transport through subtropical shallow convection in icon-lem simulations. Journal of Advances in Modeling Earth Systems, 2021.
[32] J. Jiang, Z. Lu, and Y. Liu. Experimental and theoretical studies on the counter gradient transport phenomena in turbulent flows. Advances in Mechanics, 30:425–432, 2000.
[33] X. Qiu, J. Jiang, and Y. Liu. Effects of pressure-gradient on turbulent counter-gradinet transport. Chinese Journal of Theoretical and Applied Mechanics, 36:163–170, 2004.
[34] S. Eskinazi and H. Yeh. An investigation on fully developed turbulent flows in a curved channel. Paper presented at: Aerospace Research Central Proceedings, 1954 August.
[35] C. J. Chen and S. Y. Jaw. Fundamentals of Turbulence Modeling. Physics. Cambridge University Press, England, 1998.
[36] S. V. Patankar and D. B. Spalding. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. International Journal of Heat and Mass Transfer, 15:1787–1806, 1972.
[37] J. Kim, P. Moin, and R. Moser. Turbulence statistics in fully developed channel flow at low reynolds number. Journal of Fluid Mechanics, 177:133–166, 1987.
[38] J. H. Bell and R. D. Mehta. Development of a two-stream mixing layer from tripped and untripped boundary layers. AIAA Journal, 20:2034–2042, 1990.
[39] H. J. Kaltenbach, M. Fatica, R. Mittal, T. S. Lund, and P. Moin. Study of flow in a planar asymmetric diffuser using large-eddy simulation. Journal of Fluid Mechanics, 390:151–185, 1999.
[40] A. Hajaali and T. Stoesser. Flow separation dynamics in three-dimensional asymmetric diffusers. Flow, Turbulence and Combustion, 108:973–999, 2022.
[41] H. Tang, Y. L. Lei, X. Z. Li, and Y. Fu. Large-eddy simulation of an asymmetric plane diffuser: comparison of different subgrid scale models. Symmetry, 11:1337–1354, 2019.
[42] J. Crawford and A. M. Birk. Influence of inlet boundary conditions on simulations of an asymmetric diffuser with the $v^{2}-f$ turbulence model. Paper presented at: Proceedings of ASME Turbo Expo 2015: Turbine Technical Conference and Exposition, 2015.
[43] P. A. Durbin. Separated flow computations with the $k-\varepsilon-v^{2}$ model. AIAA Journal, 33:659–664, 1995.
[44] F. R. Menter. Zonal two equation k-omega turbulence models for aerodynamic flows. Paper presented at: 24th Fluid Dynamics Conference, 1993.
[45] F. R. Menter. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32:1598–1605, 1994.
[46] C. G. Speziale. On nonlinear $k-l$ and $k-\varepsilon$ models of turbulence. Journal of Fluid Mechanics, 178:459–475, 1987.
[47] B. S. Bladwin and H. Lomax. Thin layer approximation and algebraic model for separated turbulent flows. paper presented at: Aiaa 16th Aerospace Sciences Meeting, 1978.
[48] D. A. Johnson and L. S. King. A mathematically simple turbulence closure model for attached and separated turbulent boundary layers. AIAA Journal, 23:1684–1692, 1985.
[49] R. Maulik and O. San. A neural network approach for the blind deconvolution of turbulent flows. Journal of Fluid Mechanics, 831:151–181, 2017.
[50] C. Y. Xie, X. M. Xiong, and J. C. Wang. Artificial neural network approach for turbulence models: a local framework. Physical Review Fluids, 6:084612–084638, 2021.
[51] R. Fang, D. Sondak, P. Protopapas, and S. Succi. Neural network models for the anisotropic reynolds stress tensor in turbulent channel flow. Journal of Turbulence, 21:525–543, 2020.
[52] K. Fukami, K. Fukagata, and K. Taira. Super-resolution reconstruction of turbulent flows with machine learning. Journal of Fluid Mechanics, 870:106–120, 2019.
[53] K. Fukami, K. Fukagata, and K. Taira. Machine-learning-based spatio-temporal super resolution reconstruction of turbulent flows. Journal of Fluid Mechanics, 909:1–24, 2020.
[54] B. Li, Z. Yang, X. Zhang, G. He, B. Q. Deng, and L. Shen. Using machine learning to detect the turbulent region in flow past a circular cylinder. Journal of Fluid Mechanics, 905:1–26, 2020.
[55] S. Pandey, J. Schumacher, and K. R. Sreenivasan. A perspective on machine learning in turbulent flows. Journal of Turbulence, 21:567–584, 2020.