A Self-Closed Turbulence Model for the Reynolds-Averaged Navier-Stokes Equations

1. Introduction

Solving the Reynolds-averaged Navier-Stokes (RANS) equations [1,2] to predict turbulent flows has been a central subject of turbulent modelling, in which a major challenge is turbulence model aiming to construct a closure for the Reynolds stresses [3–13], In respect of the turbulence modelling, a myriad of tentative theories, such as zero-equation models, one- or two- equations models (so called N transport equations), on the the Reynolds stresses, have been proposed and each with its own doctrines and beliefs [3–6,12,13]. Unfortunately those models include some parameters from data fitting and even functions appear in the N transport equations due to lack of universal principle [3–6,12,13]. There is a great need to develop a turbulence model that does not contain any tunable parameters.

2. The Reynolds-averaged Navier-Stokes (RANS) equations

The Navier–Stokes equations of incompressible flows can be expressed as follows:

$$\begin{array}{cc}\hfill \frac{\partial \mathit{u}}{\partial t}+\mathit{u}·(\mathsf{\nabla }\mathit{u})& =-\frac{1}{\rho }\mathsf{\nabla }p+\nu {\mathsf{\nabla }}^2\mathit{u},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \mathsf{\nabla }·\mathit{u}& =0.\hfill \end{array}$$

(2)

where $\mathit{u}(\mathit{x},t)$ is the flow velocity field, $\mathit{x}$ is the vector of spatial coordinates and t is time, $\rho$ is constant mass density, $p(\mathit{x},t)$ is flow pressure, $\mu$ is dynamical viscosity, $\mathit{x}=x^k{\mathit{e}}_k$ are position coordinates, ${\mathit{e}}_k$ is a base vector, $\mathit{u}$ is flow velocity, $\mathsf{\nabla }={\mathit{e}}_k\frac{\partial }{\partial x^k}$ is a gradient operator, and ${\mathsf{\nabla }}^2(:)=\mathsf{\nabla }·\mathsf{\nabla }(:)$.

Applying the divergence operation to both sides of the momentum equation in Eq.1 and using the mass conservation $\mathsf{\nabla }·\mathit{u}=0$ leads to:

$$\begin{array}{cc}\hfill {\mathsf{\nabla }}^{2}p& =-\rho \mathsf{\nabla }·(\mathit{u}·\mathsf{\nabla }\mathit{u}).\hfill \end{array}$$

(3)

Reynolds [1], assuming that turbulent motion already exists, sought to establish a criterion, which decides whether the turbulent character will increase or diminish, or remain stationary. Reynolds [1] proposed decomposing the flow velocity $\mathit{u}$ and pressure p into their respective time-averaged quantities $\overline{\mathit{u}}$ and $\overline{p}$ and their respective fluctuating quantities ${\mathit{u}}^{\prime }$ and $p^{\prime }$. As such, the Reynolds decompositions are

$$\begin{array}{ccc}\hfill \mathit{u}(\mathit{x},t)& =& \overline{\mathit{u}}\left(\mathit{x}\right)+{\mathit{u}}^{\prime }(\mathit{x},t),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill p(\mathit{x},t)& =& \overline{p}\left(\mathit{x}\right)+p^{\prime }(\mathit{x},t),\hfill \end{array}$$

(5)

which convert the four independent unknowns, $\mathit{u}$ and p, into the eight independent unknowns, $\overline{\mathit{u}}$, $\overline{p}$, ${\mathit{u}}^{\prime }$, and $p^{\prime }$.

According to Reynolds [1], the time-averaged velocity and pressure are defined as integration transformation as follows [2,12,13]:

$$\begin{array}{ccc}\hfill \overline{\mathit{u}}\left(\mathit{x}\right)& =& \frac{1}{T}{\int }_t^{t+T}\mathit{u}(\mathit{x},t)dt,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \overline{p}\left(\mathit{x}\right)& =& \frac{1}{T}{\int }_t^{t+T}p(\mathit{x},t)dt,\hfill \end{array}$$

(7)

where T is the time period over which the averaging takes place and must be sufficiently large to give meaningful averages to measure mean values. Naturally, the time-fluctuating velocity is ${\mathit{u}}^{\prime }=\mathit{u}-\overline{\mathit{u}}\left(\mathit{x}\right)$ and the time-fluctuating pressure is $p^{\prime }=p-\overline{p}\left(\mathit{x}\right)$, both with vanishing time averages, namely, $\overline{{\mathit{u}}^{\prime }}=\mathbf{0}$ and $\overline{p^{\prime }}=0$, respectively.

The Reynolds decomposition transforms the Navier-Stokes equations into equations for the mean velocity $\overline{\mathit{u}}$, mean pressure $\overline{p}$, velocity-fluctuation ${\mathit{u}}^{\prime }$ and fluctuation pressure $p^{\prime }$ as follows:

$$\begin{array}{c}\overline{\mathit{u}}·(\mathsf{\nabla }\otimes \overline{\mathit{u}})=-\frac{1}{\rho }\mathsf{\nabla }·\left(\overline{p}\mathit{I}\right)+\nu {\mathsf{\nabla }}^2\overline{\mathit{u}}+\frac{1}{\rho }\mathsf{\nabla }·\mathit{\tau },\end{array}$$

(8)

$$\begin{array}{c}\mathsf{\nabla }·\overline{\mathit{u}}=0,\end{array}$$

(9)

$$\begin{array}{c}{\mathit{u}}_{,t}^{\prime }+\mathsf{\nabla }·(\overline{\mathit{u}}\otimes {\mathit{u}}^{\prime }+{\mathit{u}}^{\prime }\otimes \overline{\mathit{u}}+{\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime })\hfill \\ =-\frac{1}{\rho }\mathsf{\nabla }·\left(p^{\prime }\mathit{I}\right)+\nu {\mathsf{\nabla }}^2{\mathit{u}}^{\prime }-\frac{1}{\rho }\mathsf{\nabla }·\mathit{\tau },\hfill \end{array}$$

(10)

$$\begin{array}{c}\mathsf{\nabla }·{\mathit{u}}^{\prime }=0,\end{array}$$

(11)

where $\mathit{I}={\mathit{e}}_k\otimes {\mathit{e}}_k$ is an identity tensor and the Reynolds stress tensor is defined by

$$\begin{array}{cc}\hfill \mathit{\tau }\left(\mathit{x}\right)& =-\rho {\displaystyle \overline{{\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime }}}=-\frac{\rho }{T}{\int }_t^{t+T}({\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime })dt,\hfill \end{array}$$

(12)

and

$$\begin{array}{cc}\hfill \nabla ·\mathit{\tau }& =-\rho \nabla ·\left(\overline{{\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime }}\right)=-\rho \left(\overline{{\mathit{u}}^{\prime }·(\mathsf{\nabla }\otimes {\mathit{u}}^{\prime })}\right)\hfill \\ \hfill & =-\frac{\rho }{T}{\int }_t^{t+T}{\mathit{u}}^{\prime }·(\mathsf{\nabla }\otimes {\mathit{u}}^{\prime })dt.\hfill \end{array}$$

(13)

and $\overline{{\mathit{u}}^{\prime }·(\mathsf{\nabla }\otimes {\mathit{u}}^{\prime })}=\frac{1}{T}{\int }_t^{t+T}{\mathit{u}}^{\prime }·(\mathsf{\nabla }\otimes {\mathit{u}}^{\prime })dt$.

If you look at the system of equations in Eqs. 8, 9, 10 and 11, you can easily identify that the system has eight equations for the eight unknowns that are the mean fields ${\overline{u}}_1,{\overline{u}}_2,{\overline{u}}_3,\overline{p}$ and fluctuation fields $u_1^{\prime },u_2^{\prime },u_3^{\prime },p^{\prime }$. However, the summation of Eq. 8 and 10 will cancel out the contribution of the Reynolds stress tensor, namely $\nabla ·\mathit{\tau }$. It means that the systems of equations in Eqs. 8, 9, 10 and 11 can not be directly solved without extra information of the Reynolds stress tensor $\mathit{\tau }$. We need to do some operations on the system of the equations to get rid of the terms containing the Reynolds stress.

3. A self-closed turbulence model

Since the velocity fluctuations in the Reynolds decompositions are rapid varying function respect to the time t, it 1st derivative to the time may not be able to capture its feature, it is natural to think about its 2nd derivatives to the time, because differential operation can make the function more smooth.

Since $\overline{\mathit{u}}\left(\mathit{x}\right)$, $\overline{p}\left(\mathit{x}\right)$, $\mathit{\tau }\left(\mathit{x}\right)$ and $\overline{{\mathit{u}}^{\prime }·\mathsf{\nabla }{\mathit{u}}^{\prime }}$ are the function of space coordinates $\mathit{x}$, hence their derivatives with respect to the time equal to zero, namely ${\overline{\mathit{u}}}_{,t}=0$, ${\overline{p}}_{,t}=0$ and

$$\frac{\partial \mathit{\tau }}{\partial t}={\mathit{\tau }}_{,t}\equiv 0$$

(14)

From Eqs. (10) and (11), we have following equations for the fluctuation fields

$$\begin{array}{cc}\hfill & {\mathit{u}}_{,tt}^{\prime }-\nu {\mathsf{\nabla }}^2{\mathit{u}}_{,t}^{\prime }+\frac{1}{\rho }\mathsf{\nabla }·\left(\mathit{I}{p}_{,t}^{\prime }\right)\hfill \\ \hfill & =-\mathsf{\nabla }·(\overline{\mathit{u}}\otimes {\mathit{u}}_{,t}^{\prime }+{\mathit{u}}_{,t}^{\prime }\otimes \overline{\mathit{u}}+{\mathit{u}}_{,t}^{\prime }\otimes {\mathit{u}}^{\prime }+{\mathit{u}}^{\prime }\otimes {\mathit{u}}_{,t}^{\prime }),\hfill \end{array}$$

(15)

The beauty of Eq.15 is that there is no any terms containing the Reynolds stress.

After having Eq.15, let’s us expand the Reynolds stress tensor in Taylor series. Denoting $\mathit{\sigma }={\mathit{u}}^{\prime }(\mathit{x},t)\otimes {\mathit{u}}^{\prime }(\mathit{x},t)$, and applying the product rule of integration, namely $\int udv=uv-\int vdu$, to integration ${\int }_t^{t+T}\mathit{\sigma }(\mathit{x},\xi )d\xi$, leads to

$$\begin{array}{cc}\hfill & {\int }_{t_0}^{t_0+T}\mathit{\sigma }(\mathit{x},t)dt={\int }_{t_0}^{t_0+T}\mathit{\sigma }(\mathit{x},t)d(t-t_0-T)\hfill \\ \hfill & =T\mathit{\sigma }(\mathit{x},t_0)+\frac{1}{2!}{T}^2\frac{d\mathit{\sigma }}{dt}{{|}_{t=t_0}+\frac{1}{3!}{T}^3\frac{d^2\mathit{\sigma }}{d{t}^2}|}_{t=t_0}\dots ,\hfill \\ \hfill & =\sum _{n=0}^{\infty }\frac{T^{n+1}}{(n+1)!}\left[\frac{{\partial }^n\mathit{\sigma }(\mathit{x},t)}{\partial t^n}\right]{|}_{t=t_0}.\hfill \end{array}$$

(16)

According to Reynolds [1], the turbulent motion is assumed already exists at $t_0$, the Reynolds stress tensor is hence expressed as follows

$$\begin{array}{cc}\hfill \mathit{\tau }\left(\mathit{x}\right)& =-\frac{\rho }{T}{\int }_{t_0}^{t_0+T}\mathit{\sigma }(\mathit{x},t)dt\hfill \\ \hfill & =-\rho \sum _{n=0}^{\infty }\frac{T^n}{(n+1)!}{\left[\frac{{\partial }^n({\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime })}{\partial t^n}\right]}_{t=t_0}.\hfill \end{array}$$

(17)

Since the period T must be sufficiently large to give meaningful averages to measure mean values, hence the series in Eq.17 is divergent. For a sufficient large T, the series summation can’t go to infinite and must truncated. How to fix the divergence problem ?

The period T can be considered as observable time scale, for a boundary thickness $\delta$, mass density $\rho$ and viscosity $\nu$, by dimensional analysis, then period T can be proposed as follows

$$T=\frac{\rho {\delta }^2}{\mu }=\frac{{\delta }^2}{\nu }.$$

(18)

To have some ideas about this periodic scale, for example, the kinematic viscosity of water at 20^o is about $\nu ={10}^6{\mathrm{m}}^2{\mathrm{s}}^{-1}$, and if the length scale of the fluid is $\delta =0.001$, Eq.18 gives the time scale $T={10}^3$ s. As can be seen, the time scale of such a definition is considerable. The time scale in Eq.18 indicates that, for given viscosity, the length scale $\delta$ will effect the turbulence.

If the length scale is taken as the longest ${\delta }_{max}$ of the problem, there is the maximum time $T_{max}={\delta }_{max}^2/\nu$. Chen and Sun proposed a new summation of divergent series by Möbius inversion [15].According to Chen and Sun [15], for a divergent series in Eq.17 can be expressed as

$$\begin{array}{cc}\hfill \mathit{\tau }\left(\mathit{x}\right)& \approx -\rho \sum _{n=0}^{[{\delta }_{max}/\delta ]}\frac{T^n}{(n+1)!}{\left[\frac{{\partial }^n({\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime })}{\partial t^n}\right]}_{t=t_0}.\hfill \end{array}$$

(19)

It is easy to see that there is an essential difference between Eq. 17 and Eq. 19. The conventional summation in Eq. 17 is made up to infinity, however, the modified summation in Eq. 19 is made up to a finite number $[T_{max}/T]$. The infinite terms summation leads to divergent, however, the controllable finite terms of summation bring a convergent, since the number of summation is a function of the ratio of $[T_{max}/T]$. This strategy of dealing with the divergent series in Eq. 19 has been successfully used for studying hypersonic compressible fluid by Renard et al. [16], whose proposed model has benn shown to be an efficient, accurate, and robust alternative to classic Navier–Stokes methods for the simulation of compressible flows.

In practical computations, it is very often to keep the divergent series only up to its first two terms [15]. It seems surprising that in this case the approximation solution is quite accurate with a reasonably large range of T. As an application, if the $T_{max}$ is defined as $T\le T_{max}<2T$, or in ratio $1\le T_{max}/T<2$, which leads to $[T_{max}/T]=1$, hence the Reynolds stress tensor can be approximated as follows

$$\begin{array}{cc}\hfill \mathit{\tau }\left(\mathit{x}\right)\approx & -\rho {\left[{\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime }\right]}_{t=t_0}-\frac{\rho }{2!}\frac{{\delta }^2}{\nu }{\left[\frac{\partial ({\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime })}{\partial t}\right]}_{t=t_0}.\hfill \end{array}$$

(20)

From Eq.20, we have the divergence of the Reynolds stress

$$\begin{array}{cc}\hfill \mathsf{\nabla }·\mathit{\tau }\left(\mathit{x}\right)=& -\rho {\left[{\mathit{u}}^{\prime }·(\mathsf{\nabla }\otimes {\mathit{u}}^{\prime })\right]}_{t=t_0}\hfill \\ \hfill & -\frac{\rho }{2}\frac{{\delta }^2}{\nu }{\left\{\frac{\partial [{\mathit{u}}^{\prime }·(\mathsf{\nabla }\otimes {\mathit{u}}^{\prime })]}{\partial t}\right\}}_{t=t_0}.\hfill \end{array}$$

(21)

Using alternative expression $\mathsf{\nabla }\otimes \mathit{u}=\mathsf{\nabla }\mathit{u}$, and since ${\mathit{u}}^{\prime }=0$ at $t=t_0$, hence, the self-closed turbulence model for the Reynolds-averaged Navier-Stokes equations are proposed as follows

$$\begin{array}{cc}\hfill \overline{\mathit{u}}·(\mathsf{\nabla }\overline{\mathit{u}})& =-\frac{1}{\rho }\mathsf{\nabla }·\left(\overline{p}\mathit{I}\right)+\nu {\mathsf{\nabla }}^2\overline{\mathit{u}}\hfill \\ & -\frac{1}{2}\frac{{\delta }^2}{\nu }{\left\{\frac{\partial [{\mathit{u}}^{\prime }·(\mathsf{\nabla }{\mathit{u}}^{\prime })]}{\partial t}\right\}}_{t=t_0},\hfill \end{array}$$

(22)

$$\begin{array}{c}\mathsf{\nabla }·\overline{\mathit{u}}=0,\end{array}$$

(23)

$$\begin{array}{cc}\hfill {\mathit{u}}_{,tt}^{\prime }& +\overline{\mathit{u}}·\mathsf{\nabla }{\mathit{u}}_{,t}^{\prime }+{\mathit{u}}_{,t}^{\prime }·\mathsf{\nabla }\overline{\mathit{u}}+{\mathit{u}}_{,t}^{\prime }·\mathsf{\nabla }{\mathit{u}}^{\prime }+{\mathit{u}}^{\prime }·\mathsf{\nabla }{\mathit{u}}_{,t}^{\prime }\hfill \\ & =-\frac{1}{\rho }\mathsf{\nabla }·\left(\mathit{I}{p}_{,t}^{\prime }\right)+\nu {\mathsf{\nabla }}^2{\mathit{u}}_{,t}^{\prime }-\frac{1}{\rho }\mathsf{\nabla }·\mathit{\tau }\hfill \end{array}$$

(24)

$$\begin{array}{c}\mathsf{\nabla }·{\mathit{u}}^{\prime }=0.\end{array}$$

(25)

The above system of equations has eight equations for the eight unknowns, namely the mean velocity $\overline{\mathit{u}}$, mean pressure $\overline{p}$, velocity-fluctuation ${\mathit{u}}^{\prime }$ and fluctuation pressure $p^{\prime }$. It is obvious that there is no any adjustable constant in the system of the equations. Therefore, we call the system as self-closed turbulence model of the RANS.

4. A simple model

One simple model can be made by directly omitting $\mathit{\tau }$ in Eq. (10), we have simpler governing equations

$$\begin{array}{cc}\hfill \overline{\mathit{u}}·(\mathsf{\nabla }\overline{\mathit{u}})& =-\frac{1}{\rho }\mathsf{\nabla }·\left(\overline{p}\mathit{I}\right)+\nu {\mathsf{\nabla }}^2\overline{\mathit{u}}\hfill \\ & -\frac{1}{2}\frac{{\delta }^2}{\nu }{\left\{\frac{\partial [{\mathit{u}}^{\prime }·(\mathsf{\nabla }{\mathit{u}}^{\prime })]}{\partial t}\right\}}_{t=t_0},\hfill \end{array}$$

(26)

$$\begin{array}{c}\mathsf{\nabla }·\overline{\mathit{u}}=0,\end{array}$$

(27)

$$\begin{array}{c}{\mathit{u}}_{,t}^{\prime }+\mathsf{\nabla }·(\overline{\mathit{u}}\otimes {\mathit{u}}^{\prime }+{\mathit{u}}^{\prime }\otimes \overline{\mathit{u}}+{\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime })\hfill \\ =-\frac{1}{\rho }\mathsf{\nabla }·\left(p^{\prime }\mathit{I}\right)+\nu {\mathsf{\nabla }}^2{\mathit{u}}^{\prime }\hfill \end{array}$$

(28)

$$\begin{array}{c}\mathsf{\nabla }·{\mathit{u}}^{\prime }=0.\end{array}$$

(29)

5. Validation rule of turbulence model

Now the question is how to validate the the self-closed turbulence model. Since there is no exact solution for any turbulent flows has even been found, for isotropic homogenous turbulence, Kolmogorov [14] used dimensional method and obtained $E\left(\kappa \right)=C_{\kappa }{\epsilon }^{2/3}{\kappa }^{-5/3}$, where $C_{\kappa }$ is the Kolmogorov constant, $E\left(\kappa \right)$ is turbulence energy density, $\epsilon$ is dissipation energy rate and $\kappa$ wave number. The Kolmogorov $-5/3$ law is so well established that, as noted by Rogallo and Moin [4], theoretical or numerical predictions are regarded with skepticism if they fail to reproduce the Kolmogorov $-5/3$ law. Its standing is as important as the law of wall [4,5,13]. Let’s us to see whether the self-closed turbulence model will lead to a self-validation rule.

Although the Kolmogorov $-5/3$ law was not derived from the the Reynolds-averaged Navier-Stokes (RANS) but rather from the dimensional arguments, the Kolmogorov $-5/3$ law is still the only validation law of turbulence models. The validation laws for real flow turbulence models have not been developed from the Reynolds-averaged Navier-Stokes (RANS) equations [4,5,13].

Since $\overline{\mathit{u}}={\mathit{u}}^{\prime }=\mathbf{0}$ on the boundary surfaces of the region $\partial V$ or at infinity, namely ${\int }_V\mathsf{\nabla }·(\overline{\mathit{u}}\otimes {\mathit{u}}_{,t}^{\prime }+{\mathit{u}}_{,t}^{\prime }\otimes \overline{\mathit{u}}+{\mathit{u}}_{,t}^{\prime }\otimes {\mathit{u}}^{\prime }+{\mathit{u}}^{\prime }\otimes {\mathit{u}}_{,t}^{\prime })d^3\mathit{x}={\int }_{\partial V}(\overline{\mathit{u}}\otimes {\mathit{u}}_{,t}^{\prime }+{\mathit{u}}_{,t}^{\prime }\otimes \overline{\mathit{u}}+{\mathit{u}}_{,t}^{\prime }\otimes {\mathit{u}}^{\prime }+{\mathit{u}}^{\prime }\otimes {\mathit{u}}_{,t}^{\prime })·d^2\mathit{x}=0$. Hence we have a validation rule for our turbulence model as follows:

$${\int }_V\left({\mathit{u}}_{,tt}^{\prime }-\nu {\mathsf{\nabla }}^2{\mathit{u}}_{,t}^{\prime }+\frac{1}{\rho }\mathsf{\nabla }·\left(p_{,t}^{\prime }\mathit{I}\right)\right)d^3\mathit{x}\equiv 0.$$

(30)

Obviously, the proposed validation rule in Eq.30 have an advantage, namely the validation laws do not include the Reynolds stress tensor and successfully bypassed the closure issue for the Reynolds stresses, which reveals that the validation laws in Eqs.30 will not be effected by turbulence modelling on the Reynolds stresses, while the Kolmogorov law $E\left(\kappa \right)=C_{\kappa }{\epsilon }^{2/3}{\kappa }^{-5/3}$ is linked with the diagonal of the Reynolds stresses and will be effected by the stresses models due to the fact that ${\int }_0^{\kappa }E\left(\kappa \right)d\kappa =\frac{1}{2}\overline{{\mathit{u}}^{\prime }\otimes {\mathit{u}}^{\prime }}=\frac{1}{2}(\overline{u_1^{\prime 2}}+\overline{u_1^{\prime 2}}+\overline{u_1^{\prime 2}})$.

It is worth pointing out that our validation laws are rigorously derived from the Reynolds-averaged Navier-Stokes equations and applicable to real fluids, while Kolmogorov "-5/3" law was not derived from the Reynolds-averaged Navier-Stokes equations and supposed to be only valid for isotropic homogenous turbulence.

6. Two dimensional turbulent boundary layers

As an application, let’s us consider a two dimensional turbulent boundary layers. A thin flat plate is immersed at zero incidence in a uniform stream as shown in Figure 1, which flows with speed $U\left(x\right)$ and is assumed not to be affected by the presence of the plate, except in the boundary layer. The fluid is supposed unlimited in extent, and the origin of coordinates is taken at the leading edge, with x measured downstream along the plate and y perpendicular to it. Assuming that the turbulent flow is steady with pressure gradient along the x axis. The pressure gradient drives the against the shear stresses at the wall.

The unsteady Navier-Stokes equations of the two dimensional boundary layers flow under gradient are given by

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}& =0,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}& =U\frac{dU}{dx}+\nu \frac{d^{2}u}{d{y}^2},\hfill \end{array}$$

(32)

Introducing the Reynolds velocity decomposition $u=\overline{u}+u^{\prime }$ and $v=\overline{v}+v^{\prime }$, and substituting them into Eq.31 and Eq.32, then applying the Reynolds average operation and notice ${\overline{u}}_{,t}\equiv 0$, we can get the mean field equation as follows

$$\overline{u}{\overline{u}}_{,x}+\overline{v}{\overline{u}}_{,y}=U{U}_{,x}+\nu {\overline{u}}_{,yy}-(\overline{u^{\prime }{u}_{,x}^{\prime }}+\overline{v^{\prime }{u}_{,y}^{\prime }}),$$

(33)

and mean mass conservation

$${\overline{u}}_{,x}+{\overline{v}}_{,y}=0,$$

(34)

and fluctuation field equation

$$\begin{array}{cc}\hfill u_{,t}^{\prime }& +\overline{u}{u}_{,x}^{\prime }+u^{\prime }(\overline{u}+u^{\prime }){,}_{,x}+\overline{v}{u}_{,y}^{\prime }+v^{\prime }(\overline{u}+u^{\prime }){,}_{,y}\hfill \\ \hfill & =\nu u_{,yy}^{\prime }+(\overline{u^{\prime }{u}_{,x}^{\prime }}+\overline{v^{\prime }{u}_{,y}^{\prime }})\hfill \end{array}$$

(35)

as well as fluctuation mass conservation

$$u_{,x}^{\prime }+v_{,y}^{\prime }=0,$$

(36)

The boundary conditions:

$$\begin{array}{cc}\hfill y& =0:\overline{u}=\overline{v}=0,u^{\prime }=0,v^{\prime }=0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill y& =\delta \left(x\right):\overline{u}=U\left(x\right),u^{\prime }=0,v^{\prime }=0,\hfill \end{array}$$

(38)

Taking into account of the fluctuation mass conservation in Eq.36, the Reynolds turbulent stress, $\overline{u^{\prime }{u}_{,x}^{\prime }}+\overline{v^{\prime }{u}_{,y}^{\prime }}$, can be rewritten as follows: $\overline{u^{\prime }{u}_{,x}^{\prime }}+\overline{v^{\prime }{u}_{,y}^{\prime }}=\overline{{\left(u^{\prime }{u}^{\prime }\right)}_{,x}}+\overline{{\left(v^{\prime }{u}^{\prime }\right)}_{,y}}-\overline{u^{\prime }(u_{,x}^{\prime }+v_{,y}^{\prime }})=\overline{{\left(u^{\prime }{u}^{\prime }\right)}_{,x}}+\overline{{\left(v^{\prime }{u}^{\prime }\right)}_{,y}}={\overline{\left(u^{\prime }{u}^{\prime }\right)}}_{,x}+{\overline{\left(v^{\prime }{u}^{\prime }\right)}}_{,y}$.

Introducing a mean stream function $\psi (x,y)$ and express the mean velocity components as follows

$$\overline{u}={\psi }_{,y},\overline{v}=-{\psi }_{,x},$$

(39)

with the relation in Eq.(41), the mean mass conservation Eq.(34) is satisfied, and the mean momentum conservation Eq.(33) becomes

$${\psi }_{,y}{\psi }_{,xy}-{\psi }_{,x}{\psi }_{,yy}=U{U}_{,x}+\nu {\psi }_{,yyy}-(\overline{u^{\prime }{u}_{,x}^{\prime }}+\overline{v^{\prime }{u}_{,y}^{\prime }})$$

(40)

Similarly, introducing a mean stream function $\varphi (x,y)$ and express the fluctuation velocity components as follows

$$u^{\prime }={\varphi }_{,y},v^{\prime }=-{\varphi }_{,x},$$

(41)

the fluctuation mass conservation Eq.(36) is satisfied, and the fluctuation momentum conservation Eq.(35) becomes

$$\begin{array}{cc}\hfill & {\varphi }_{,t}+{\varphi }_{,y}{\varphi }_{,xy}-{\varphi }_{,x}{\varphi }_{,yy}\hfill \\ \hfill & +{\psi }_{,y}{\varphi }_{,xy}-{\varphi }_{,x}{\psi }_{xy}-{\psi }_{,x}{\varphi }_{,yy}-{\varphi }_{,x}{\psi }_{,yy}\hfill \\ \hfill & =\nu {\varphi }_{,yyy}+(\overline{u^{\prime }{u}_{,x}^{\prime }}+\overline{v^{\prime }{u}_{,y}^{\prime }})\hfill \end{array}$$

(42)

By using similar transformations proposed by Sun [10,11], we have

$$\begin{array}{cc}\hfill \psi & =U\left(x\right)\delta \left(x\right)f\left(\eta \right),\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \varphi & =U\left(x\right)\delta \left(x\right)g(\eta ,\tau ),\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill \eta & =\frac{y}{\delta \left(x\right)},\tau =\frac{\nu }{{\delta }^2}t,\hfill \end{array}$$

(45)

Some useful derivatives of function $\varphi$ respect to both x and y, for simplification, we denote ${\varphi }_{,x}=\frac{\partial \varphi }{\partial x}$ and $g_{,\eta }=\frac{\partial g}{\partial \eta }$ and so on. Noting the ${\eta }_{,x}=-\eta {\delta }^{-1}{\delta }_{,x}$ and ${\tau }_{,x}=-2\tau {\delta }^{-1}{\delta }_{,x}$, and by the chain rule for derivatives, we can obtain some useful relations as follows:

$$\begin{array}{cc}\hfill {\psi }_{,x}& ={\left(U\delta \right)}_{,x}f-U\eta {\delta }_{,x}{f}_{,\eta },\hfill \\ \hfill {\psi }_{,y}& =U{f}_{,\eta },\hfill \\ \hfill {\psi }_{,xy}& =U_{,x}{f}_{,\eta }-\eta U{\delta }^{-1}{\delta }_{,x}{f}_{,\eta \eta },\hfill \\ \hfill {\psi }_{,yy}& =U{\delta }^{-1}{f}_{,\eta \eta },\hfill \\ \hfill {\psi }_{,yyy}& =U{\delta }^{-2}{f}_{,\eta \eta \eta },\hfill \end{array}$$

(46)

and

$$\begin{array}{cc}\hfill {\varphi }_{,x}& ={\left(U\delta \right)}_{,x}g-U\eta {\delta }_{,x}{g}_{,\eta }-2U\tau {\delta }_{,x}{g}_{,\tau },\hfill \\ \hfill {\varphi }_{,y}& =U{g}_{,\eta },\hfill \\ \hfill {\varphi }_{,xy}& =U_{,x}{g}_{,\eta }-\eta U{\delta }^{-1}{\delta }_{,x}{g}_{,\eta \eta }-2U\tau {\delta }^{-1}{\delta }_{,x}{g}_{,\eta \tau },\hfill \\ \hfill {\varphi }_{,yy}& =U{\delta }^{-1}{g}_{,\eta \eta },\hfill \\ \hfill {\varphi }_{,yyy}& =U{\delta }^{-2}{g}_{,\eta \eta \eta },\hfill \\ \hfill {\varphi }_{,ty}& =U\nu {\delta }^{-2}{g}_{,\tau \eta }.\hfill \end{array}$$

(47)

Thus the mean velocity components become

$$\begin{array}{cc}\hfill \overline{u}& =U{f}_{,\eta },\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill \overline{v}& =-\left[{\left(U\delta \right)}_{,x}f-U\eta {\delta }_{,x}{f}_{,\eta }\right].\hfill \end{array}$$

(49)

and the fluctuation velocity components become

$$\begin{array}{cc}\hfill u^{\prime }& =U{g}_{,\eta },\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill v^{\prime }& =-\left[{\left(U\delta \right)}_{,x}g-U\eta {\delta }_{,x}{g}_{,\eta }-2U\tau {\delta }_{,x}{g}_{,\tau }\right].\hfill \end{array}$$

(51)

Substituting Eqs.46 into Eq.40, we have

$$\begin{array}{cc}\hfill f_{,\eta \eta \eta }& +\alpha f{f}_{,\eta \eta }+\beta [1-{\left(f_{,\eta }\right)}^2]\hfill \\ \hfill & =\beta \overline{{\left(g_{,\eta }\right)}^2}-\alpha \overline{g{g}_{,\eta \eta }}+\gamma \overline{\tau (g_{,\tau }{g}_{,\eta \eta }-g_{,\eta }{g}_{,\tau \eta })},\hfill \end{array}$$

(52)

Substituting Eqs.47 into Eq.42, we have

$$\begin{array}{cc}\hfill g_{,\eta \eta \eta }& +\alpha (g{g}_{,\eta \eta }+f{g}_{,\eta \eta }+g{f}_{,\eta \eta })\hfill \\ \hfill & -\beta [{\left(g_{,\eta }\right)}^2+2f_{,\eta }{g}_{,\eta }]\hfill \\ \hfill & =g_{,\tau \eta }+\gamma \tau [g_{,\tau }{(f+g)}_{,\eta \eta }-g_{,\tau \eta }{(f+g)}_{,\eta }]\hfill \\ \hfill & -\beta \overline{{\left(g_{,\eta }\right)}^2}+\alpha \overline{g{g}_{,\eta \eta }}-\gamma \overline{\tau (g_{,\tau }{g}_{,\eta \eta }-g_{,\eta }{g}_{,\tau \eta })},\hfill \end{array}$$

(53)

where the coefficients are $\alpha =\frac{\delta }{\nu }\frac{dU\delta }{dx}$, $\beta =\frac{{\delta }^2}{\nu }\frac{dU}{dx}$, and $\gamma =\frac{U}{\nu }\frac{d{\delta }^2}{dx}$.

If the coefficient $\alpha$, $\beta$ and $\gamma$ were constants, the Eq.52 is solvable because it will get rid of the variable x. Following from the relation $\frac{d{\delta }^{2}U}{dx}=(2\alpha -\beta )\nu$, its integration leads to ${\delta }^{2}U=(2\alpha -\beta )\nu x$, canceling out ${\delta }^2$ by $\beta =\frac{{\delta }^2}{\nu }\frac{dU}{dx}$, we have $\frac{dU}{U}=\frac{\beta }{(2\alpha -\beta )}\frac{dx}{x}$, hence when $\beta \ne 2\alpha$, the general solution of this equation is of the form:

$$\begin{array}{cc}\hfill U\left(x\right)& =C{x}^m,\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill \delta \left(x\right)& ={\left[\nu \left|\frac{(2\alpha -\beta )x}{U\left(x\right)}\right|\right]}^{1/2},\hfill \end{array}$$

(55)

where the exponent $m=\frac{\beta }{2\alpha -\beta }$. In the same way, we have $\frac{d{\delta }^{2}U}{dx}=(\beta +\gamma )\nu$, leads to a relation $2\alpha -\beta =\beta +\gamma$, namely $\gamma =2(\alpha -\beta )$. It implies that $\gamma$ is constant if both $\alpha$ and $\beta$ were constants.

With the boundary layer thickness $\delta \left(x\right)$, we can estimate the period $T=\frac{{\delta }^2}{\nu }=\frac{(2\alpha -\beta )x}{U\left(x\right)}=\frac{2\alpha -\beta }{C}{x}^{1-m}$, general speaking, the period is the function of time. For example, plate boundary layer, $U=U_0,\alpha =1,\beta =0$, hence $T=\frac{2x}{U_0}$.

Considering that, the Reynolds turbulent stress, $-\beta \overline{{\left(g_{,\eta }\right)}^2}+\alpha \overline{g{g}_{,\eta \eta }}-\gamma \overline{\tau (g_{,\tau }{g}_{,\eta \eta }-g_{,\eta }{g}_{\tau \eta })}$, is independent of time, in order to eliminate it in Eq.57, we can take the partial derivative of both sides of Eq.57 with respect to dimensionless time $\tau$, noting $f{\left(\eta \right)}_{,\tau }=0$, and obtain the equation as follows:

$$\begin{array}{cc}\hfill g_{,\eta \eta \eta \tau }& +\alpha {\left(g{g}_{,\eta \eta }\right)}_{,\tau }+\alpha (f{g}_{,\eta \eta \tau }+g_{,\tau }{f}_{,\eta \eta })\hfill \\ \hfill & -2\beta (g_{,\eta }+f_{,\eta })g_{,\eta \tau }\hfill \\ \hfill & =g_{,\tau \eta \tau }+\gamma [g_{,\tau }{(f+g)}_{,\eta \eta }-g_{,\tau \eta }{(f+g)}_{,\eta }]\hfill \\ \hfill & +\gamma \tau {\left[g_{,\tau }{(f+g)}_{,\eta \eta }-g_{,\tau \eta }{(f+g)}_{,\eta }\right]}_{,\tau }.\hfill \end{array}$$

(56)

The equations Eq.52 and Eq.56 together proper inial boundary conditions form a governing equation of the turbulent boundary layers.

Regarding the computational strategy of Eq.52 and Eq.56, the algorithm is proposed in the Table 1

7. Somple model and plate turbulent boundary layers

A simpler model can be proposed in Eq.57, omitting $-\beta \overline{{\left(g_{,\eta }\right)}^2}+\alpha \overline{g{g}_{,\eta \eta }}-\gamma \overline{\tau (g_{,\tau }{g}_{,\eta \eta }-g_{,\eta }{g}_{\tau \eta })}$, hence we have

$$\begin{array}{cc}\hfill g_{,\eta \eta \eta }& +\alpha (g{g}_{,\eta \eta }+f{g}_{,\eta \eta }+g{f}_{,\eta \eta })\hfill \\ \hfill & -\beta [{\left(g_{,\eta }\right)}^2+2f_{,\eta }{g}_{,\eta }]\hfill \\ \hfill & =g_{,\tau \eta }+\gamma \tau [g_{,\tau }{(f+g)}_{,\eta \eta }-g_{,\tau \eta }{(f+g)}_{,\eta }].\hfill \end{array}$$

(57)

For plate boundary layers in an uniform flow with velocity $U=U_0$ constant. We have $\alpha =1,\beta =0,\gamma =2$, and $T=\frac{2x}{U_0}$. The governing equations of the plate boundary layers are

$$\begin{array}{cc}\hfill f_{,\eta \eta \eta }& +f{f}_{,\eta \eta }=-\overline{g{g}_{,\eta \eta }}+2\overline{\tau (g_{,\tau }{g}_{,\eta \eta }-g_{,\eta }{g}_{,\tau \eta })},\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill g_{,\eta \eta \eta }& +(f+g)g_{,\eta \eta }+g{f}_{,\eta \eta }\hfill \\ \hfill & =g_{,\tau \eta }+2\tau [g_{,\tau }{(f+g)}_{,\eta \eta }-g_{,\tau \eta }{(f+g)}_{,\eta }].\hfill \end{array}$$

(59)

The system of equations Eq.58 and Eq.59 can be solved by iterative approach:

• Step 0: set $g_0=0$
• Step 1: to find Blasius solution $f_{0,\eta \eta \eta }+f_0{f}_{0,\eta \eta }=0$
• Step 2: to find $g_1$ by equation $g_{1,\eta \eta \eta }+(f_0+g_1)g_{1,\eta \eta }+g_1{f}_{0,\eta \eta }=g_{1,\tau \eta }+2\tau [g_{1,\tau }{(f_0+g_1)}_{,\eta \eta }-g_{1,\tau \eta }{(f_0+g_1)}_{,\eta }]$
• Step 3: to find $f_1$ by equation $f_{1,\eta \eta \eta }+f_1{f}_{1,\eta \eta }=-\overline{g_1{g}_{1,\eta \eta }}+2\overline{\tau (g_{1,\tau }{g}_{1,\eta \eta }-g_{1,\eta }{g}_{1,\tau \eta })}$

8. Conclusions

By introducing a period of time scale, and applying the product rule of integration and divergent series truncated summation, the Reynolds stresses have been expressed in an explicit form, a self-closed turbulence model for the Reynolds-averaged Navier-Stokes equations has been successfully formulated. The proposed self-closed turbulence model does not contain any adjustable parameter. The validation rule for the self-closed model is derived rigorously.

Although we theoretically propose a self-closed turbulence model, its correctness are still to be validated.