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 problem in the classic physics." Historically, some renowned physicists and mathematicians, such as Werner Heisenberg, Richard Feynman, and Andrei Kolmogorov, dedicated their wonderful 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{\partial }{\partial t}}$, ${(:)}_{,i}={\displaystyle \frac{\partial }{\partial 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{v}_{i,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.

So what is the cause of turbulence? This is a core issue in fluid mechanics. Although we are not yet able to predict the occurrence of turbulence with complete accuracy, we do know under what circumstances turbulence will definitely not occur, and that is in the flow of inviscid fluids, which cannot be turbulent. Currently, our understanding of turbulent flow can be macroscopically stated as follows: the interaction between fluid viscosity and velocity gradients is the main cause of turbulence. Viscosity causes the flow with non-zero velocity gradients to produce vorticity, leading to fluid rotation. Due to mass conservation, the fluid must replenish the mass that flows out. Additionally, the nonlinear modulation of the convection term makes the motion pattern 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.

To obtain the governing equations for turbulent flow, let’s derive the momentum conservation equation in fluid mechanics from the "first principles." The term "first principles" originates from ancient Greek philosophy, where Aristotle first introduced this concept in his work "Metaphysics." He believed that first principles are "the most basic propositions or assumptions that cannot be further derived," serving as the ultimate foundation for all knowledge and reasoning. That is, by returning to the essence of things, analyzing the most fundamental principles and assumptions, we can construct new cognitive frameworks or solutions. Elon Musk [32] is a key figure who popularized this term. In an interview, he expressed high praise for the "first principles" thinking approach: "By using first principles, I distill the essence of things and derive from the most fundamental level..." "It is crucial to think from the perspective of first principles rather than relying on comparative thinking. We often fall into the habit of comparing and doing what others have done or are doing, which only leads to incremental improvements. The first principles thinking is to examine the world from a physicist’s perspective, stripping away the external appearance of things layer by layer to insights into their inner essence, and then building from the essence step by step." This is how Musk understands the "first principles thinking model" - tracing the roots of things and rethinking the way of action. The Navier-Stokes equations have been established for over 200 years since 1822, if we want to make breakthroughs, we must start thinking from scratch with the mindset of first principles.

2. General Equations of Motion for an Incompressible Viscous Fluid

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{\partial \mathit{v}}{\partial t}}+\mathit{v}·(\nabla \mathit{v})\right]=\mathit{\sigma }·\nabla ,$$

(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 $\nabla ={\mathit{e}}_i{\partial }_i$ is the gradient operator and $\mathit{\sigma }·\nabla$ is right gradient of the stress tensor $\mathit{\sigma }$, whose computational rule is $\mathit{\sigma }·\nabla =({\sigma }_{ij}{\mathit{e}}_i\otimes {\mathit{e}}_j)·\left({\mathit{e}}_k{\partial }_k\right)={\sigma }_{ij,k}{\mathit{e}}_i({\mathit{e}}_j·{\mathit{e}}_k)={\sigma }_{ij,k}{\mathit{e}}_i{\delta }_{jk}={\sigma }_{ij,j}{\mathit{e}}_i$, in general $\mathit{\sigma }·\nabla \ne \nabla ·\mathit{\sigma }$. Equation (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}$, namely $\mathit{\sigma }=\mathit{\sigma }\left(\mathit{S}\right)$, where $\mathit{S}={\displaystyle \frac{1}{2}}(\nabla \mathit{v}+\mathit{v}\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 [33,34,35,36,37], 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\hfill \\ \hfill & =-p\mathit{I}+\mu (\nabla \mathit{v}+\mathit{v}\nabla )+\lambda {(\nabla \mathit{v}+\mathit{v}\nabla )}^2\hfill \\ \hfill & =\left[-p{\delta }_{ij}+\mu (v_{i,j}+v_{j,i})+\lambda (v_{i,k}+v_{k,i})(v_{k,j}+v_{j,k})\right]{\mathit{e}}_i\otimes {\mathit{e}}_j\hfill \end{array}$$

(2)

where ${\mathit{S}}^2=\mathit{S}·\mathit{S}=(S_{ik}{\mathit{e}}_i\otimes {\mathit{e}}_k)·(S_{lj}{\mathit{e}}_l\otimes {\mathit{e}}_j)=S_{ik}{S}_{lj}({\mathit{e}}_k·{\mathit{e}}_l){\mathit{e}}_i\otimes {\mathit{e}}_j=S_{ik}{S}_{lj}{\delta }_{kl}{\mathit{e}}_i\otimes {\mathit{e}}_j=S_{il}{S}_{lj}{\mathit{e}}_i\otimes {\mathit{e}}_j=[(v_{i,k}+v_{k,i}){\mathit{e}}_i\otimes {\mathit{e}}_k]·[(v_{l,j}+v_{j,l}){\mathit{e}}_l\otimes {\mathit{e}}_j]=(v_{i,k}+v_{k,i})(v_{k,j}+v_{j,k}){\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 [33,34,35,36,37]: $\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 given by

$$\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, $\mathit{A}:\mathit{B}$ is double dot product of two tensor $\mathit{A}=A_{ij}{\mathit{e}}_i\otimes {\mathit{e}}_j$ and $\mathit{B}=B_{ij}{\mathit{e}}_i\otimes {\mathit{e}}_j$ with computational rule: $\mathit{A}:\mathit{B}=(A_{ij}{\mathit{e}}_i\otimes {\mathit{e}}_j):(B_{kl}{\mathit{e}}_k\otimes {\mathit{e}}_l)=A_{ij}{B}_{kl}({\mathit{e}}_i·{\mathit{e}}_k)\otimes ({\mathit{e}}_j·{\mathit{e}}_l)=A_{ij}{B}_{kl}{\delta }_{ik}{\delta }_{jl}=A_{ij}{B}_{ij}$

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 [38], which was published after formulation of the Naver-Stokes equations. The second reason is that even if one could obtain the expression Equation (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 [39] 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 Equation (2) into Equation (1) allows us to derive the fluid motion equation that considers second-order viscosity:

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

(4)

wherein, the kinematic viscosity is $\nu =\mu /\rho$, and the kinematic second-order viscosity is $\kappa =\lambda /\rho$. Equation (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 $\nabla ·\mathit{v}=0$, Equation (4) can be simplified to:

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

(5)

where the Laplacian ${\nabla }^2=\nabla ·\nabla$. In this way, we have improved the Navier-Stokes equation to Equation (5), which is the most general equation of motion for incompressible viscous fluids. The term $\kappa \left[{(\nabla \mathit{v}+\mathit{v}\nabla )}^2\right]·\nabla$ includes the nonlinear combination of the fluid’s second-order viscosity and velocity gradients, which is the main cause of complex flow phenomena. Since Equation (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.

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

$$\begin{array}{cc}\hfill v_{i,t}+v_j{v}_{i,j}& =-{\displaystyle \frac{1}{\rho }}{p}_{,i}+\nu v_{i,kk}+\kappa {\left[(v_{i,k}+v_{k,i})(v_{k,j}+v_{j,k})\right]}_{,j},\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 Equation (5) not only yields the Reynolds number $Re={\displaystyle \frac{LV}{\nu }}$, but Equation (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 Equation (6) can be further expanded in conventional form, in the Cartesian rectangular coordinate system, denoting $v_1=u,v_2=v,v_3=w$, $x_1=x,x_2=y,x_3=z$, ${(:)}_{,x}={\displaystyle \frac{\partial (:)}{\partial x}}$, ${(:)}_{,yy}={\displaystyle \frac{{\partial }^2(:)}{\partial y^2}}$, ${\nabla }^2={\displaystyle \frac{{\partial }^2}{\partial x^2}}+{\displaystyle \frac{{\partial }^2}{\partial y^2}}+{\displaystyle \frac{{\partial }^2}{\partial z^2}}$, 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{1}{\rho }}{p}_{,x}+\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{1}{\rho }}{p}_{,y}+\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{1}{\rho }}{p}_{,z}+\nu {\nabla }^{2}w+\kappa M_z,\hfill \end{array}$$

(9)

where

$$\begin{array}{cc}\hfill M_x& ={[4{\left(u_{,x}\right)}^2+{(v_{,x}+u_{,y})}^2+{(w_{,x}+u_{,z})}^2]}_{,x}\hfill \\ \hfill & +{[2(u_{,x}+v_{,y})(u_{,y}+v_{,x})+(w_{,y}+v_{,z})(w_{,x}+u_{,z})]}_{,y}\hfill \\ \hfill & +{[2(u_{,x}+w_{,z})(u_{,z}+w_{,x})+(v_{,z}+w_{,y})(v_{,x}+u_{,y})]}_{,z},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill M_y& ={[2(u_{,x}+v_{,y})(u_{,y}+v_{,x})+(w_{,x}+u_{,z})(w_{,y}+v_{,z})]}_{,x}\hfill \\ \hfill & +{[{(u_{,y}+v_{,x})}^2+4{\left(v_{,y}\right)}^2+{(w_{,y}+v_{,z})}^2]}_{,y}\hfill \\ \hfill & +{[(u_{,z}+w_{,x})(u_{,y}+v_{,x})+2(v_{,y}+w_{,z})(v_{,z}+w_{,y})]}_{,z},\hfill \end{array}$$

(11)

and

$$\begin{array}{cc}\hfill M_z& ={[2(u_{,x}+w_{,z})(u_{,z}+w_{,x})+(v_{,x}+u_{,y})(v_{,z}+w_{,y})]}_{,x}\hfill \\ \hfill & +{[(u_{,y}+v_{,x})(u_{,z}+w_{,x})+2(v_{,y}+w_{,z})(v_{,z}+w_{,y})]}_{,y}\hfill \\ \hfill & +{[{(u_{,z}+w_{,x})}^2+{(v_{,z}+w_{,y})}^2+4{\left(w_{,z}\right)}^2]}_{,z}.\hfill \end{array}$$

(12)

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}$$

(13)

$$\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}$$

(14)

For heat conduction problems, it is assumed that there is a heat flux density $\mathit{q}=-\chi \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{\partial s}{\partial t}}+\mathit{v}·\nabla s)=& \nabla ·(\chi \nabla T)+2\mu \mathit{S}:\mathit{S}+4\lambda {\mathit{S}}^2:\mathit{S},\hfill \end{array}$$

(15)

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

3. Boundary Layer Theory and Application to Wedge Flow

In 1904, Ludwig Prandtl [39] 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 [40]. In particular, Prandtl discovered that the flow within the boundary layer is also turbulent [41,42], and in 1925, Prandtl proposed the concept of the mixing length to characterize the turbulent viscosity [43], 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 [40] studied the problem of steady flat-plate boundary layers. A thin flat plate is immersed at zero incidence in a uniform stream with speed $U\left(x\right)=const.$ 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. Blasius [40] 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 [44] 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}$$

(16)

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

(17)

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. The exact solution of partial differential equation systems in Equation (16) and Equation (17) can be easily solved by using Maple, the solutions of velocity components $u,$ and $v$ are mainly expressed in terms of Lambert function [46]: LambertW=$W\left(z\right)$, which satisfies: $W\left(z\right)exp\left(W\right(z\left)\right)=z$. As the equation $W\left(z\right)exp\left(W\right(z\left)\right)=z$has an infinite number of solutions y for each (non-zero) value of z, LambertW has an infinite number of branches. Exactly one of these branches is analytic at 0.

Note $U\left(x\right)=const.$ and $U_{,x}=0$, the Maple code of solving above partial differential equations is provided in Appendix A. Hence, denoting $F=ln[c_1(-x+y)+c_{2}t+c_3]$, we have the general solutions of Equation (16) and Equation (17) for $U\left(x\right)=U_0=const.$ as below:

$$\begin{array}{c}\hfill u=c_4+{\int }^{F}exp\left[\xi -W\left(-{\displaystyle \frac{2c_1\kappa }{\nu }}exp\left[{\displaystyle \frac{c_2}{\nu }}(c_3-{\mathrm{e}}^{\xi })\right]\right)+{\displaystyle \frac{c_2}{\nu }}(c_3-{\mathrm{e}}^{\xi })\right]d\xi ,\end{array}$$

(18)

and

$$\begin{array}{cc}\hfill v=& {\int }^{F}exp\left[\xi -W\left(-{\displaystyle \frac{2c_1\kappa }{\nu }}exp\left[{\displaystyle \frac{c_2}{\nu }}(c_3-{\mathrm{e}}^{\xi })\right]\right)+{\displaystyle \frac{c_2}{\nu }}(c_3-{\mathrm{e}}^{\xi })\right]d\xi \hfill \\ & +c_4-{\displaystyle \frac{c_2}{c_1}}-{\displaystyle \frac{c_1{c}_2}{1+W(-\frac{2\kappa c_1}{\nu }exp(-\frac{c_2\left(c_1\left(-x+y\right)+c_{2}t\right)}{\nu }))}}\hfill \\ \hfill & +{\displaystyle \frac{2\kappa c_1^2{c}_2}{\nu }}{\displaystyle \frac{exp[W(-\frac{2\kappa c_1}{\nu }exp(-\frac{c_2(c_1(-x+y)+c_{2}t)}{\nu }))+\frac{-c_2[c_{2}t-(x-y)c_1]}{\nu }]}{1+W(-\frac{2\kappa c_1}{\nu }exp[-\frac{c_2{c}_1(-x+y)+c_{2}t)}{\nu }])}}.\hfill \end{array}$$

(19)

Although we have obtained the solution, it is regrettable that the integral in Equation (18) cannot be calculated explicitly at the moment, making this solution less applicable. Let us try other methods below.

Introducing a stream function $\psi (x,y,t)$ and express the velocity components as follows: $u={\displaystyle \frac{\partial \psi }{\partial y}}={\psi }_{,y},v=-{\displaystyle \frac{\partial \psi }{\partial x}}=-{\psi }_{,x}$, the mass conservation Equation (17) is satisfied, and the momentum conservation Equation (16) becomes

$${\psi }_{,ty}+{\psi }_{,y}{\psi }_{,xy}-{\psi }_{,x}{\psi }_{yy}=U{U}_{,x}+\nu {\psi }_{,yyy}+2\kappa {\psi }_{,yy}{\psi }_{,xyy}$$

(20)

and corresponding boundary conditions.

Applying Sun’s similarity transformations [44]: $\psi =U\left(x\right)\delta \left(x\right)f(\eta ,\tau ),\eta ={\displaystyle \frac{y}{\delta \left(x\right)}},\tau ={\displaystyle \frac{\nu }{{\delta }^2}}t,$. Denoting ${\psi }_{,x}={\displaystyle \frac{\partial \psi }{\partial x}}$ and $f_{,\eta }={\displaystyle \frac{\partial f}{\partial \eta }}$, and 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: ${\psi }_{,x}={\left(U\delta \right)}_{,x}f-U\eta {\delta }_{,x}{f}_{,\eta }-2U\tau {\delta }_{,x}{f}_{,\tau },{\psi }_{,y}=U{f}_{,\eta },{\psi }_{,xy}=U_{,x}{f}_{,\eta }-\eta U{\delta }^{-1}{\delta }_{,x}{f}_{,\eta \eta }-2U\tau {\delta }^{-1}{\delta }_{,x}{f}_{,\eta \tau },{\psi }_{,yy}=U{\delta }^{-1}{f}_{,\eta \eta },{\psi }_{,yyy}=U{\delta }^{-2}{f}_{,\eta \eta \eta },{\psi }_{,ty}=U\nu {\delta }^{-2}{f}_{,\tau \eta }$. Thus the velocity components become

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

(21)

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

(22)

Substituting Equation (21) and Equation (22) into Equation (20), we have a single partial differential equation as follows

$$\begin{array}{cc}\hfill f_{,\eta \eta \eta }& +\alpha f{f}_{,\eta \eta }+\beta [1-{\left(f_{,\eta }\right)}^2]\hfill \\ \hfill & +{\displaystyle \frac{2\kappa }{\nu }}{f}_{,\eta \eta }\left[U_{,x}{f}_{,\eta \eta }-{\displaystyle \frac{U{\delta }_{,x}}{\delta }}{f}_{,\eta \eta }-{\displaystyle \frac{U{\delta }_{,x}}{\delta }}{f}_{,\eta \eta \eta }\right]\hfill \\ \hfill & =f_{,\tau \eta }+\gamma \tau (f_{,\tau }{f}_{,\eta \eta }-f_{,\eta }{f}_{,\tau \eta }),\hfill \end{array}$$

(23)

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

If the coefficient $\alpha$, $\beta$ and $\gamma$ were constants, then we have: $U\left(x\right)=C{x}^m,\delta \left(x\right)={\left[\nu \left|{\displaystyle \frac{(2\alpha -\beta )x}{U\left(x\right)}}\right|\right]}^{1/2}$, where the exponent $m={\displaystyle \frac{\beta }{2\alpha -\beta }}$. It is clear, for arbitrary $U\left(x\right)=C{x}^m$, that the Equation (23) explicitly contains the coordinate x, which cann’t be solved easily for arbitrary exponent m.

In the case of wedge flow with velocity field: $U\left(x\right)=Cx$, so that $m=1$, then $\alpha =\beta$ and $\gamma =0$, and leads to a constant boundary thickness $\delta ={[\nu \alpha x/\left(Cx\right)]}^{1/2}={(\alpha \nu /C)}^{1/2}$. If we set $\alpha =1$, the Equation (23) becomes

$$\begin{array}{cc}\hfill f_{,\eta \eta \eta }& +f{f}_{,\eta \eta }+1-{\left(f_{,\eta }\right)}^2+\Lambda {\left(f_{,\eta \eta }\right)}^2=f_{,\tau \eta },\hfill \end{array}$$

(24)

and boundary conditions: $\eta =0:f=f_{,\eta }=0$, $\eta \to \infty :f_{,\eta }=1$ and $f(\tau ,0)=0$, where $\Lambda ={\displaystyle \frac{2\kappa C}{\nu }}$ is a constant for given a problem. Equation (24) has get rid of coordinate x, can be solved numerically. A full Maple code of solving Equation (24) is provided in Appendix B.

For its corresponding steady problem, the function f is only a function of $\eta$, the Equation (24) can be simplified into:

$$\begin{array}{cc}\hfill f_{,\eta \eta \eta }& +f{f}_{,\eta \eta }+1-{\left(f_{,\eta }\right)}^2+\Lambda {\left(f_{,\eta \eta }\right)}^2=0,\hfill \end{array}$$

(25)

and boundary conditions: $\eta =0:f=f_{,\eta }=0$ and $\eta \to \infty :f_{,\eta }=1$. Equation (25) not only facilitates solving the equation but also allows us to express the obtained results as a single profile in terms of the similar variable $\eta$. A full Maple code of solving Equation (25) is provided in Appendix C.

Although we do not yet have specific data for $\kappa$, we can still assume $\Lambda =0$ for the laminar boundary layer and $\Lambda =0.95$ for the turbulent boundary layer. Using the program provided above, it is easy to calculate the corresponding $f,f_{,\eta },f_{,\eta \eta }$ for both laminar and turbulent flow. The specific results are shown in Figure 1. From the profile of $f_{,\eta }$, it can be seen that the velocity profile in turbulent flow is blunter than that in laminar flow; from the profile of $f_{,\eta \eta }$, it can be observed that the friction in turbulent flow is much greater than that in laminar flow.

To demonstrate the difference in frictional resistance between turbulent and laminar flow, let’s calculate the shear stress ${\tau }_{xy}=\mu (u_{,y}+v_{,x})+2\lambda \left({\underbrace{u_{,x}+v_{,y}}}_{=0}\right)(u_{,y}+v_{,x})=\mu (u_{,y}+v_{,x})$. Considering the order of magnitude approximation estimates of the Prandtl boundary layer theory in Table 1, the shear stress can be approximated as: ${\tau }_{xy}\approx \mu u_{,y}$. Therefore, the local wall shear stress ${\tau }_w$ is given by ${\tau }_{\mathrm{w}}={\tau }_{xy}{|}_{y=0}=\mu u_{,y}{{|}_{y=0}=\mu {\displaystyle \frac{U}{\delta }}{f}_{,\eta \eta }\left(\eta \right)|}_{\eta =0}$:

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{w}}& =x\rho \sqrt{\nu C^3}{f}_{,\eta \eta }\left(0\right)=\left\{\begin{array}{cc}1.233x\rho \sqrt{\nu C^3},& \mathrm{Laminar}\mathrm{flow}:\Lambda =0,\\ 1.923x\rho \sqrt{\nu C^3},& \mathrm{Turbulent}\mathrm{flow}:\Lambda =0.95,\end{array}\right.\hfill \end{array}$$

(26)

This expression reveals that the resistance of turbulent flow is much greater than that of laminar flow. Table 2 lists the different $\Lambda$ corresponding to different $f_{\eta \eta }\left(0\right)$; the larger the $\Lambda$, the greater the $f_{\eta \eta }\left(0\right)$.

4. Discussions and Conclusions

Based on the above theoretical derivation and application examples, the quadratic term $\kappa {\mathit{S}}^2$ plays a fundamental role in characterizing turbulent motion. It can be said that the Navier-Stokes equations without the quadratic term $\kappa {\mathit{S}}^2$ are actually equations that primarily describe laminar flow and cannot be used to describe turbulent motion. Due to the high velocity gradients in turbulent motion, the quadratic term $\kappa {\mathit{S}}^2$ cannot be omitted and must be retained; only then can the derived equations be used to describe turbulent motion.

Although the context of this paper is directed towards turbulence, the equations Equation (5) or Equation (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, such as gas aerodynamics and sound generated by turbulence. Using the constitutive equation Equation (2), the motion equations for compressible viscous fluids can be established by the same method. 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.