Fundamental Equations for Turbulent Motion of an Incompressible Viscous Fluid

1. Introduction

The mention of turbulence immediately brings to mind a saying: turbulence is one of the most challenging problems in physics and is known as an "unsolved classic physics problem." Historically, some renowned physicists and mathematicians, such as Werner Heisenberg, Richard Feynman, and Andrei Kolmogorov, dedicated their entire lives to the problem of turbulence without solving it. It is even said that the famous physicist Horace Lamb hoped to ask God about turbulence after reaching heaven [1,2]. From the perspective of quantitative analysis, the daunting scientific problem of turbulence began in 1895 with Reynolds [3] attempting to study the Navier-Stokes equations for viscous fluids using average statistical methods, namely Reynolds averaged Navier-Stokes (RANS) equations, which has been ongoing for 130 years. Over these 130 years, it can be said that all available mathematical tools and computational resources have been employed to study the Navier-Stokes equations for viscous fluids, yet no fundamental progress or breakthrough has been made [4,5,6,7,8]. Since the establishment of the Reynolds averaged Navier-Stokes (RANS) equations, due to the needs of numerous engineering applications, a variety of approximation methods have been employed to close the Reynolds averaged Navier-Stokes (RANS) equation set. Various models have been proposed for the Reynolds stresses [9,10,11,12,13,14,15,16,17,18]. Is it possible to establish a universal closure turbulence model theory that can predict turbulent motion across various scales? Unfortunately, after 130 years of collective effort, such a theory has yet to be established. On the contrary, the numerous models proposed to date have all relied on non-rigorous hypotheses usually based on experimental observations and have their own limitations [17,18,19].

Although direct numerical simulation (DNS) has made gratifying progress [20,21,22,23,24], and in recent years, the use of big data and AI to study turbulence has emerged, the essence of turbulent motion has not yet been revealed. Standing at the dawn of the AI era and looking to the future, we want to loudly ask, where is the path to solving turbulence? Is it only possible to unravel the mysteries of turbulence through large-scale numerical simulations when quantum computers are developed? Over the past 130 years of numerical simulations, regardless of the computational strategy, a fundamental fact is that all simulations are based on the Navier-Stokes equations, as it is believed that they encompass all the information of turbulent motion [2,8,15,17,18,25,26]. Now, we are faced with two lines of thought: one is to continue on the path based on the Navier-Stokes equations, and the other is to return to the origin, to re-examine the equations of fluid motion from the starting point of their establishment, that is, to reflect on whether the Navier-Stokes equations need to be further refined and improved to encompass both turbulent and laminar motion without using the Reynolds statistical averaging methods [3].

To cast aside the clouds and see the sun, let us take a general look back at the establishment of the equations of fluid motion. Euler played a crucial role in conceptualizing the mathematical description of fluid flow. He described the flow using a three-dimensional pressure and velocity field that varies in space, modeling the flow as a collection of continuous, infinitesimally small fluid elements. By applying the basic principles of conservation of mass and Newton’s second law, Euler arrived at two coupled nonlinear partial differential equations involving the flow fields of pressure and velocity. Although these Euler equations represent a significant intellectual breakthrough in theoretical fluid dynamics, obtaining their general solutions is a difficult task. Euler did not consider the effect of frictional forces on the movement of fluid elements; in other words, he ignored viscosity, and the Euler equations are not suitable for solving real flow problems. It was only a century after Euler that his equations were modified to account for the influence of frictional forces within the flow field. The resulting set of equations is a more complex system of nonlinear partial differential equations now known as the Navier-Stokes equations, initially derived by Navier in 1822 [27] and later independently by Stokes in 1845 [28]. That is, the Navier-Stokes Equations (system): the momentum conservation equation: $v_{i,t}+v_j{v}_{i,j}=-{\displaystyle \frac{1}{\rho }}{p}_{,i}+\nu v_{i,kk}$, and the incompressible mass conservation equation: $v_{i,i}=0$, where $v_i$ is the flow velocity field, p is the flow pressure, t is time, $\rho$ is the constant mass density, and $\nu$ is the kinematic viscosity, ${(:)}_{,t}={\displaystyle \frac{𝜕}{𝜕t}}$, ${(:)}_{,i}={\displaystyle \frac{𝜕}{𝜕x_i}}$ and $x_i$ is spatial coordinates. To this day, the Navier-Stokes equations remain the gold standard for the mathematical description of fluid flow. Unfortunately, their general analytical solutions have not yet been found [29,30,31].

In 1895, Reynolds studied turbulent motion using statistical averaging methods, decomposing the fluid velocity into the sum of mean velocity ${\overline{v}}_i$ and fluctuating velocity ${\tilde{v}}_i$, that is, $v_i=\overline{v_i}+\widetilde{v_i}$. By time-averaging the Navier-Stokes equations, he obtained the Reynolds-Averaged Navier-Stokes (RANS) equations, which are expressed as ${\overline{v}}_{i,t}+{\overline{v}}_j{\overline{v}}_{i,j}=-{\displaystyle \frac{1}{\rho }}{\overline{p}}_{,i}+\nu {\overline{v}}_{i,kk}-\overline{{\tilde{v}}_i{\tilde{v}}_j}$, with the hope of obtaining the average values of important parameters in fluid motion. However, due to the appearance of the fluctuating velocity correlation term $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ (termed Reynolds stress but is not a stress at all! [19]) in the Reynolds-Averaged Navier-Stokes (RANS) equations, which are not closed in themselves.

It is not difficult to notice that various turbulence models are all based on the Reynolds-Averaged Navier-Stokes (RANS) equations, with the aim of guessing the relationship between the turbulent viscosity coefficient ${\nu }_t$ and the mean velocity field from different perspectives using the Boussinesq hypothesis $\overline{{\tilde{v}}_i{\tilde{v}}_j}={\displaystyle \frac{2}{3}}k{\delta }_{ij}-{\nu }_t\left({\overline{v}}_{i,j}+{\overline{v}}_{j,i}\right)$, where k is the turbulent kinetic energy density. In other words, the Boussinesq hypothesis essentially treats $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ as a kind of stress and, by mimicking the linear constitutive relationship of Newtonian viscous fluids, artificially constructs a possible relationship between $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ and the mean flow velocity. The reason why people propose various models for $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ is that they believe turbulence is caused by $\overline{{\tilde{v}}_i{\tilde{v}}_j}$, which is actually a misconception.

In fact, $\overline{{\tilde{v}}_i{\tilde{v}}_j}$originates from the convective term in the fluid motion acceleration, which is $v_{j}vi,j$, and is produced by the Reynolds averaging process. That is, $\overline{v_{i,t}+v_j{v}_{i,j}}$ = $\overline{({\overline{v}}_{i,t}+{\tilde{v}}_{i,t})+({\overline{v}}_j+{\tilde{v}}_j)({\overline{v}}_{i,j}+{\tilde{v}}_{i,j})}$ = ${\overline{v}}_{i,t}+\overline{{\tilde{v}}_j{\tilde{v}}_{i,j}}$. Using the mass conservation relationship $v_{i,i}=0$, the Reynolds average of the fluid motion acceleration can be rewritten as $\overline{v_{i,t}+v_j{v}_{i,j}}={\overline{v}}_{i,t}+{\left(\overline{{\tilde{v}}_j{\tilde{v}}_i}\right)}_{,j}$. This means that the generation of$\overline{{\tilde{v}}_i{\tilde{v}}_j}$ has nothing to do with fluid viscosity; it is merely the projection of fluid velocity onto its velocity gradient. In fact, the inviscid Euler equation $v_{i,t}+v_j{v}_{i,j}=-p_{,i}/\rho$ ($\nu v_{i,jj}$ is omitted for inviscid flow) will also produce $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ after Reynolds averaging. Therefore, George [18] stated that $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ is not a stress at all, but simply a re-worked version of the fluctuation contribution to the nonlinear acceleration terms. And we know that inviscid fluids cannot generate turbulent motion, so it can be said that $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ is not the cause of turbulence.

What, then, is the cause of turbulence? This is a fundamental question in fluid mechanics. Macroscopically, it can be said that the interaction between fluid viscosity and velocity gradients is the main cause of turbulence. Viscosity causes flows with non-zero velocity gradients to produce vorticity, leading to fluid rotation. Due to the conservation of mass, the fluid must replenish the mass that flows out, and the nonlinear modulation of the convective terms makes the patterns of motion very complex. When the velocity gradient is very small, the flow is mainly laminar; when the velocity gradient is very large, the flow is mainly turbulent.

Within the current framework of the Navier-Stokes equations, only one term, $\nu v_{i,kk}$, in the Navier-Stokes equation $v_{i,t}+v_j{v}_{i,j}=-p_{,i}/\rho +\nu v_{i,kk}$ contains the viscosity coefficient $\nu$ and the derivative of the velocity gradient $v_{i,k}$, which is the divergence $v_{i,kk}$. Unfortunately, because $\overline{{\tilde{v}}_{i,kk}}=0$, after Reynolds averaging, the divergence $v_{i,kk}$ only retains the mean field ${\overline{v}}_{i,kk}$, without including the fluctuating components. This may imply that the Navier-Stokes equations established based on $\nu v_{i,kk}$, $v_{i,t}+v_j{v}_{i,j}=-{\displaystyle \frac{1}{\rho }}{p}_{,i}+\nu v_{i,kk}$ and $v_{i,i}=0$, may not include the complete information of viscous fluid motion and are imperfect, necessitating a modification. Since $\nu v_{i,kk}$ is derived from the linear constitutive relation of viscous fluids combined with the mass conservation condition $v_{j,j}=0$, that is, $\nu {(v_{i,j}+v_{j,i})}_j=\nu v_{i,kk}$, it is necessary to re-examine the constitutive relation of fluids in order to include the complete information of fluid motion.

2. General Equations of Motion for an Incompressible Viscous Fliud

We know that before introducing the fluid constitutive relation, the mass conservation equation for incompressible fluids is $v_{i,i}=0$, and the momentum conservation equation is: $\rho (v_{i,t}+v_j{v}_{i,j})={\sigma }_{ij,j}$. The momentum conservation equation can also be rewritten in the form of tensor total quantity as follows:

$$\rho \left[{\displaystyle \frac{𝜕\mathit{v}}{𝜕t}}+\mathit{v}·\left(\mathsf{\nabla }\mathit{v}\right)\right]=\mathsf{\nabla }·\mathit{\sigma },$$

(1)

where $\mathit{\sigma }={\sigma }_{ij}{\mathit{e}}_i\otimes {\mathit{e}}_j$ is the stress tensor, ${\sigma }_{ij}$ are the components of the stress tensor, ${\mathit{e}}_i$ are the unit base vectors in the Cartesian coordinate system, and $\mathsf{\nabla }={\mathit{e}}_i{𝜕}_i$ is the gradient operator. Eq.(1) is the most general form of the momentum balance equation for a continuum and can be applied to describe the motion of any continuum.

To obtain a closed set of governing equations, we need to introduce the fluid’s constitutive relation, which is the relationship between the stress tensor $\mathit{\sigma }$ and the rate-of-strain tensor $\mathit{S}$: $\mathit{\sigma }=\mathit{\sigma }\left(\mathit{S}\right)$, where $\mathit{S}={\displaystyle \frac{1}{2}}(\mathsf{\nabla }\mathit{v}+\mathit{v}\mathsf{\nabla })=S_{ij}{\mathit{e}}_i\otimes {\mathit{e}}_j$ is the rate-of-strain tensor, and its components are $S_{ij}={\displaystyle \frac{1}{2}}(v_{i,j}+v_{j,i})$.

It should be noted that in the process of introducing the fluid’s constitutive relation, there is an element of human interpretation involved in understanding the physical properties of the fluid, which is the most subjective part in establishing the equations of fluid motion. From Navier to Stokes, it is assumed that the fluid is isotropic and incompressible, and a linear constitutive relationship between the stress tensor $\mathit{\sigma }$ and the rate-of-strain tensor $\mathit{S}$ is hypothesized: $\mathit{\sigma }=-p\mathit{I}+2\mu \mathit{S}$, where $\mu$ is the viscosity coefficient. This constitutive relation was proposed in an era when there was no understanding of tensors and has been adopted by all subsequent literature. In the modern era, equipped with knowledge of tensors, we must ask whether the constitutive relationship between the stress tensor $\mathit{\sigma }$ and the rate-of-strain tensor $\mathit{S}$ can be rationally expressed, and if so, what form it should take.

Assuming the fluid’s physical properties are isotropic, since both the stress tensor $\mathit{\sigma }$ and the rate-of-strain tensor $\mathit{S}$ are second-order tensors, according to the tensor representation theory [32,33,34,35,36], the most general expression for the constitutive relation of an incompressible fluid $\mathit{\sigma }=\mathit{\sigma }\left(\mathit{S}\right)$ can always be written as follows:

$$\begin{array}{cc}\hfill \mathit{\sigma }& =-p\mathit{I}+2\mu \mathit{S}+4\lambda {\mathit{S}}^2=-p\mathit{I}+\mu (\mathsf{\nabla }\mathit{v}+\mathit{v}\mathsf{\nabla })+\lambda {(\mathsf{\nabla }\mathit{v}+\mathit{v}\mathsf{\nabla })}^2,\hfill \end{array}$$

(2)

where ${\mathit{S}}^2=\mathit{S}·\mathit{S}=(S_{ik}{\mathit{e}}_i\otimes {\mathit{e}}_k)·(S_{\mathcal{l}j}{\mathit{e}}_{\mathcal{l}}\otimes {\mathit{e}}_j)=S_{ik}{S}_{\mathcal{l}j}({\mathit{e}}_k·{\mathit{e}}_{\mathcal{l}}){\mathit{e}}_i\otimes {\mathit{e}}_j=S_{ik}{S}_{\mathcal{l}j}{\delta }_{k\mathcal{l}}{\mathit{e}}_i\otimes {\mathit{e}}_j=S_{i\mathcal{l}}{S}_{\mathcal{l}j}{\mathit{e}}_i\otimes {\mathit{e}}_j$, the term $\lambda$ is not referred to by name in the literature; in this text, it is termed the second-order viscosity coefficient. From the thermodynamic entropy inequality, the inequality can be obtained as cited in [32,33,34,35,36]: $\mu \mathrm{tr}{\mathit{S}}^2+2\lambda \mathrm{tr}{\mathit{S}}^3\ge 0$. By utilizing the Cayley-Hamilton theorem, this inequality can be further simplified to: $\mu \mathrm{tr}{\mathit{S}}^2+6\lambda det\left(\mathit{S}\right)\ge 0$.

The fluid viscosity leads to the dissipation of energy, which is ultimately converted into heat. The energy dissipation per unit time for an incompressible viscous fluid is

$$\begin{array}{cc}\hfill {\dot{E}}_d=& -{\int }_{\Omega }(2\mu \mathit{S}+4\lambda {\mathit{S}}^2):\mathit{S}{d}^3\mathit{x},\hfill \end{array}$$

(3)

where $\Omega$ represents the volume configuration of the fluid control body.

It is particularly worth noting that there is an issue that does not need to be hidden, which is the data problem of the physical parameter $\kappa$. Currently, there are no data available for the parameter $\kappa$, and it is a research topic that needs to be measured in the future.

In the literature of fluid mechanics, from Navier and Stokes to the present, fluids have been regarded as Newtonian fluids, using the linear constitutive relation $\mathit{\sigma }=-p\mathit{I}+2\mu \mathit{S}$, without considering the second-order term $\lambda {\mathit{S}}^2$. As for why the second-order term $\lambda {\mathit{S}}^2$ is not considered, there may be two reasons. The first is that during the time of Newton, Navier, and Stokes, the nonlinear effects were not observable, and there were no tools for tensor analysis, let alone tensor representation theory and the Cayley-Hamilton theorem. The second reason is that even if one could obtain the expression Eq.(2), it was considered that the second-order term $\lambda {\mathit{S}}^2$ is small enough to be negligible.

This paper argues that neglecting the second-order term $\lambda {\mathit{S}}^2$ is a significant oversight in the establishment of fluid mechanics equations. Because the spatial variation of velocity, such as the velocity gradient near a wall, is very large, neither the first nor the second power of the velocity gradient can be omitted. This is analogous to the situation in the boundary layer. Prandtl [37] stated when proposing the boundary layer theory that the Reynolds number $Re$ in ${\displaystyle \frac{1}{Re}}{\nabla }^2\mathit{v}$ cannot be omitted, no matter how large $Re$ (or how small ${\displaystyle \frac{1}{Re}}$) is, because the velocity gradient is so large that the entire term ${\displaystyle \frac{1}{Re}}{\nabla }^2\mathit{v}$ cannot be neglected.

Bringing Eq.2 into Eq.1 allows us to derive the fluid motion equation that considers second-order viscosity:

$${\displaystyle \frac{𝜕\mathit{v}}{𝜕t}}+\mathit{v}·\left(\mathsf{\nabla }\mathit{v}\right)=-{\displaystyle \frac{1}{\rho }}\mathsf{\nabla }p+2\mathsf{\nabla }·\left(\nu \mathit{S}\right)+4\mathsf{\nabla }·\left(\kappa {\mathit{S}}^2\right),$$

(4)

wherein, the kinematic viscosity is $\nu =\mu /\rho$, and the kinematic second-order viscosity is $\kappa =\lambda /\rho$. Eq.4 represents the most general equation of motion for viscous fluids. The viscosity coefficients $\nu ,\kappa$ are generally functions of pressure p and temperature T, and they cannot be factored out of the gradient operator.

In most cases, the viscosity coefficients $\nu ,\kappa$ do not vary significantly within the fluid and can be considered as constants. Under the incompressible condition $\mathsf{\nabla }·\mathit{v}=0$, Eq.4 can be simplified to:

$${\displaystyle \frac{𝜕\mathit{v}}{𝜕t}}+\mathit{v}·\left(\mathsf{\nabla }\mathit{v}\right)=-{\displaystyle \frac{1}{\rho }}\mathsf{\nabla }p+\nu {\nabla }^2\mathit{v}+\kappa \mathsf{\nabla }·\left[{(\mathsf{\nabla }\mathit{v}+\mathit{v}\mathsf{\nabla })}^2\right],$$

(5)

where the Laplacian ${\nabla }^2=\mathsf{\nabla }·\mathsf{\nabla }$. In this way, we have improved the Navier-Stokes equation to Eq.5, which is the most general equation of motion for incompressible viscous fluids. The term $\mathsf{\nabla }·\left[\kappa {(\mathsf{\nabla }\mathit{v}+\mathit{v}\mathsf{\nabla })}^2\right]$ includes the nonlinear combination of the fluid’s second-order viscosity and velocity gradients, which is the main cause of complex flow phenomena. Since Eq. (5) already includes higher-order terms of the velocity gradient, it can be directly used to simulate turbulence without the need to employ Reynolds’ statistical averaging method, which is the statistical decomposition of the fluid velocity.

Eq.(5) can be written in component form as follows:

$$\begin{array}{cc}\hfill v_{j,t}+v_i{v}_{j,i}& =-{\displaystyle \frac{1}{\rho }}{p}_{,j}+\nu v_{j,kk}+\kappa {\left[(v_{k,i}+v_{i,k})(v_{k,j}+v_{j,k})\right]}_{,i}.\hfill \end{array}$$

(6)

It should be particularly noted that for problems with specific characteristic length $L$ and characteristic velocity $V$, the dimensionless process of Eq.5 not only yields the Reynolds number $Re={\displaystyle \frac{LV}{\nu }}$, but Eq.5 also generates a new dimensionless parameter, denoted as $K={\displaystyle \frac{L^2}{\kappa }}$. The parameter $K$ arises due to the consideration of second-order viscosity. For a fluid with a given $\kappa$, this parameter $K$ is independent of the characteristic velocity and only related to the characteristic length $L$, or in other words, for flow problems of a given length scale, the parameter $K$ affects the fluid motion for all characteristic velocity scale. Therefore, it is necessary to retain the second-order viscosity coefficient.

The tensorial Eq.6 can be further expanded in conventional form, in the Cartesian rectangular coordinate system, the three-dimensional fluid momentum equation can be read as follows:

$$\begin{array}{cc}\hfill u_{,t}+u{u}_{,x}+v{u}_{,y}+w{u}_{,z}=& -{\displaystyle \frac{p_{,x}}{\rho }}+\nu {\nabla }^{2}u+\kappa M_x,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill v_{,t}+u{v}_{,x}+v{v}_{,y}+w{v}_{,z}=& -{\displaystyle \frac{p_{,y}}{\rho }}+\nu {\nabla }^{2}v+\kappa M_y,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill w_{,t}+u{w}_{,x}+v{w}_{,y}+w{w}_{,z}=& -{\displaystyle \frac{p_{,z}}{\rho }}+\nu {\nabla }^{2}w+\kappa M_z,\hfill \end{array}$$

(9)

where $v_1=u,v_2=v{v}_3=w$, $x_1=x,x_2=y,x_3=z$, ${(:)}_{,x}={\displaystyle \frac{𝜕(:)}{𝜕x}}$, ${(:)}_{,yy}={\displaystyle \frac{{𝜕}^2(:)}{𝜕y^2}}$, ${\nabla }^2={\displaystyle \frac{{𝜕}^2}{𝜕x^2}}+{\displaystyle \frac{{𝜕}^2}{𝜕y^2}}+{\displaystyle \frac{{𝜕}^2}{𝜕z^2}}$, $M_x={[4{\left(u_{,x}\right)}^2+{(v_{,x}+u_{,y})}^2+{(w_x+u_{,z})}^2]}_{,x}+{[2u_{,x}(u_{,y}+v_{,x})+2v_{,y}(v_{,x}+u_{,y})+(w_{,y}+v_{,z})(w_{,x}+u_{,z})]}_{,y}+{[2u_{,x}(u_{,z}+w_{,x})+(v_{,z}+w_{,y})(v_{,x}+u_{,y})+2w_{,z}(w_{,x}+u_{,z})]}_{,z}$, $M_y={[2u_{,x}(u_{,y}+v_{,x})+2v_{,y}(v_{,x}+u_{,y})+(w_x+u_{,z})(w_{,y}+v_{,z})]}_{,x}+{[{(u_{,y}+v_{,x})}^2+4{\left(v_{,y}\right)}^2+{(w_{,y}+v_{,z})}^2]}_{,y}+{[(u_{,z}+w_{,x})(u_{,y}+v_{,x})+2v_{,y}(v_{,z}+w_{,y})+2w_{,z}(w_{,y}+v_{,z})]}_{,z}$, and $M_z={[2u_{,x}(u_{,z}+w_{,x})+(v_{,y}+u_{,y})(v_{,z}+w_{,y})+2w_{,z}(w_x+u_{,z})]}_{,x}+{[(u_{,y}+v_{,x})(u_{,z}+w_{,x})+2v_{,y}(v_{,z}+w_{,y})+2w_{,z}(w_{,y}+v_{,z})]}_{,y}+{[{(u_{,z}+w_{,x})}^2+{(v_{,z}+w_{,y})}^2+4{\left(w_{,z}\right)}^2]}_{,z}$. The difference from the Navier-Stokes equations can be seen from the above equations, namely, each equation includes higher-order terms of the velocity gradient.

After using the continuity equation $u_{,x}+v_{,y}=0,$ the two-dimensional momentum equation is expressed as follows:

$$\begin{array}{cc}\hfill u_{,t}+u{u}_{,x}+v{u}_{,y}=& -{\displaystyle \frac{1}{\rho }}{p}_{,x}+\nu (u_{,xx}+u_{,yy})+\kappa {[4{\left(u_{,x}\right)}^2+{(u_{,y}+v_{,x})}^2]}_{,x}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill v_{,t}+u{v}_{,x}+v{v}_{,y}=& -{\displaystyle \frac{1}{\rho }}{p}_{,y}+\nu (v_{,xx}+v_{,yy})+\kappa {[4{\left(v_{,y}\right)}^2+{(u_{,y}+v_{,x})}^2]}_{,y},\hfill \end{array}$$

(11)

For heat conduction problems, it is assumed that there is a heat flux density $\mathit{q}=-\chi \mathsf{\nabla }T$, where T is the temperature and $\chi$ is the thermal conductivity. The general equation of heat transfer is given by

$$\begin{array}{cc}\hfill \rho T({\displaystyle \frac{𝜕s}{𝜕t}}+\mathit{v}·\mathsf{\nabla }s)=& \mathsf{\nabla }·\left(\chi \mathsf{\nabla }T\right)+2\mu \mathit{S}:\mathit{S}+4\lambda {\mathit{S}}^2:\mathit{S},\hfill \end{array}$$

(12)

where s is the entropy density per mass [7].

3. Boundary Layer Theory

In 1904, Ludwig Prandtl [37] introduced the concept of the boundary layer at the third International Congress of Mathematicians, suggesting that the influence of frictional forces is confined to a thin boundary layer near the surface of an object. Within the boundary layer, the velocity gradient is large, as is the shear stress, leading to frictional resistance that cannot be neglected. The boundary layer equations are much simpler than the Navier-Stokes equations and can be solved numerically [38]. In particular, Prandtl discovered that the flow within the boundary layer is also turbulent [39,40], and in 1925, Prandtl proposed the concept of the mixing length to characterize the turbulent viscosity [41], which can be used to derive the logarithmic law for the wall in turbulent boundary layers. This was a significant attempt to close the Reynolds-Averaged Navier-Stokes (RANS) equations.

In 1908, Prandtl’s student Blasius studied the problem of steady flat-plate boundary layers. He performed a similarity transformation on Prandtl’s equations $u{u}_{,x}+v{u}_{,y}=U{U}_{,x}+\nu u_{,yy}$ and $u_{,x}+v_{,y}=0$ and successfully obtained a numerical solution for the problem. In 2024, Sun [42] studied unsteady laminar boundary layers and applied a similarity transformation to the equations $u_{,t}+u{u}_{,x}+v{u}_{,y}=U{U}_{,x}+\nu u_{,yy}$ and $u_{,x}+v_{,y}=0$, successfully obtaining an exact solution for the unsteady flat-plate laminar boundary layer. This solution is primarily expressed in terms of Kummer functions, and other flow problems were also investigated.

Based on the Prandtl’s magnitude analysis and approximations as shown in Table 1, we can obtain the boundary layer equations as follows:

$$\begin{array}{cc}\hfill u_{,t}+u{u}_{,x}+v{u}_{,y}& =U{U}_{,x}+\nu u_{,yy}+\underline{2\kappa u_{,y}{u}_{,xy}},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill u_{,x}+v_{,y}& =0,\hfill \end{array}$$

(14)

in which, the underlined segment represents the manifestation of the second-order theory proposed in this paper, which is where it differs from Prandtl’s boundary layer equations. Another difference is that here the fluid velocity u is used instead of the velocity average $\overline{u}$ as in Prandtl’s equations.

4. Discussions and Conclusions

Although the context of this paper is directed towards turbulence, the equations Eq. 5 or Eq. 6 established here are indeed general equations of motion for viscous fluids. The derivation of these equations is based on rationality at every step, without any artificial introduction of approximations. These equations can be applied to various movements of incompressible viscous fluids. The only regret is that the data for the physical parameter $\kappa$ has not been found in the current literature, and it seems that its determination will be a subject for future research.

This article serves as a significant inspiration for understanding the causes of turbulence due to its innovative reevaluation of the fluid motion equations, addressing a long-standing issue in turbulence research. By retaining second-order terms in the fluid’s constitutive equation, it challenges conventional wisdom and provides a more accurate representation of turbulent flows. This work not only advances our understanding of fluid dynamics but also offers a new framework for future research, potentially revolutionizing the derivation of hydrodynamic equations for turbulent motion. Its significance lies in its potential to refine the Navier-Stokes equations, making it a pivotal contribution to the field.