Fundamental Equations for the Turbulent Motion of an Incompressible Viscous Fluid

1. Introduction

Turbulence is a term that swiftly conjures up a well-known adage: it stands as one of the most formidable challenges in physics and is often cited as one of the "unsolved problems in classical physics." The quest to unravel turbulence has consumed the brilliant careers of notable physicists and mathematicians, including Werner Heisenberg, Richard Feynman, and Andrei Kolmogorov, yet it remains unsolved. It is whimsically legacy that the esteemed physicist Horace Lamb intended to inquire about turbulence in the afterlife, should he meet the divine [1,2].

The quantitative quest to understand turbulence dates back to 1895 with Reynolds’ [3] pioneering efforts to analyze the Navier-Stokes equations for viscous fluids through the lens of average statistical methods, giving rise to the Reynolds averaged Navier-Stokes (RANS) equations. This endeavor has spanned over a century and a third, with every imaginable mathematical apparatus and computational resource brought to bear on the problem. Despite this, no seminal advance or breakthrough has emerged [4,5,6,7,8]. Since the formulation of the RANS equations, the demand of numerous engineering applications has led to the deployment of a spectrum of approximation techniques to close the RANS equation set. A multitude of models have been suggested to account for the elusive Reynolds stresses [9,10,11,12,13,14,15,16,17,18]. The question remains: can a universal closure theory for turbulence be crafted that accurately predicts turbulent motion across a range of scales? Regrettably, after 130 years of collective scientific endeavor, such a theory remains elusive. Instead, the plethora of models proposed to date are underpinned by non-rigorous assumptions, typically grounded in experimental observations, and each carries its own set of limitations [17,18,19].

Despite the significant strides made in direct numerical simulation (DNS) of turbulence, as documented in studies by Orszag [20], Rogallo et al.[21,22], Spalart [23], and Moin [24], and the recent advent of big data and artificial intelligence (AI) techniques in turbulence research, the intrinsic nature of turbulent flow remains an enigma. As we stand at the cusp of the AI era, gazing towards the horizon of discovery, we are compelled to pose a resounding question: What is the trajectory to the resolution of turbulence? Must we await the advent of quantum computing to unravel the enigmas through expansive numerical simulations? Across the span of 130 years of numerical simulation history, a pivotal truth endures: all simulation efforts are grounded in the Navier-Stokes equations, which are considered to encapsulate the entirety of turbulent phenomena, as evidenced by the work of [2,8,15,17,18,26]

Confronted with this legacy, we are at a crossroads with two distinct avenues of thought: one that advocates for a continued pursuit along the path defined by the Navier-Stokes equations, and another that suggests a return to the foundations, a reevaluation of the equations of fluid motion from their inception. This reflective approach challenges whether the Navier-Stokes equations should be further refined or enhanced to account for both turbulent and laminar flows without resorting to the Reynolds averaging methods, as proposed by Reynolds in 1895 [3].

In pursuit of clarity amidst the shadows, let us survey the historical landscape that gave rise to the equations governing fluid motion. Euler’s seminal contribution laid the groundwork for the mathematical depiction of fluid flow, envisioning it as a three-dimensional field of pressure and velocity that evolves through space. He conceptualized the flow as a continuum of infinitesimal fluid elements, deriving two interwoven nonlinear partial differential equations from the principles of mass conservation and Newton’s second law. While Euler’s equations mark a pivotal advance in theoretical fluid dynamics, their general solutions prove elusive. Notably, Euler’s framework omitted the frictional forces affecting fluid element’s viscosity-rendering his equations inadequate for real-world flow phenomena.

It took another century for the incorporation of frictional forces into Euler’s equations, culminating in the more intricate system we now recognize as the Navier-Stokes equations. Introduced independently by Navier in 1822 and Stokes in 1845, this system expands upon Euler’s with a set of nonlinear partial differential equations that account for the internal friction within the flow field. The Navier-Stokes equations, encapsulated in 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$ denotes the flow velocity field, p the pressure, t time, $\rho$ the mass density, and $\nu$ the kinematic viscosity, stand as the quintessential model for fluid flow description. Despite their preeminence, the quest for general analytical solutions to the Navier-Stokes equations persists, an endeavor that has eluded resolution to this day, as noted by [29,30,31].

In the pivotal year of 1895, Reynolds embarked on a groundbreaking investigation into turbulent flow by employing statistical averaging techniques. He deconstructed the fluid velocity into its mean component ${\overline{v}}_i$ and its fluctuating component ${\tilde{v}}_i$, such that the total velocity $v_i$ is given by the sum $v_i=\overline{v_i}+\widetilde{v_i}$. Through the process of time-averaging the Navier-Stokes equations, he derived the Reynolds-Averaged Navier-Stokes (RANS) equations, formulated 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 objective of capturing the mean values of key parameters in fluid dynamics. Despite this, the RANS equations are complicated by the inclusion of the fluctuating velocity correlation term $\overline{{\tilde{v}}_i{\tilde{v}}_j}$, referred to as Reynolds stress, despite not truly being a stress [19]-rendering the equations intrinsically unclosed.

The prevalence of turbulence models grounded in the Reynolds-Averaged Navier-Stokes (RANS) equations is immediately apparent, each model striving to infer the connection between the turbulent viscosity coefficient ${\nu }_t$ and the mean velocity field through diverse lenses. Central to these efforts is the Boussinesq hypothesis, which posits $\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 denotes the turbulent kinetic energy density. Essentially, the Boussinesq hypothesis conceptualizes $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ as a stress-like quantity and, by emulating the linear constitutive relationship of Newtonian fluids, fabricates a potential linkage between $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ and the mean flow velocities. The motivation behind the myriad models proposed for $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ stems from a belief that turbulence originates from this term $\overline{{\tilde{v}}_i{\tilde{v}}_j}$, a notion that is, in fact, a misapprehension.

Indeed, the term $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ arises from the convective acceleration within the fluid motion, encapsulated by $v_j{v}_{i,j}$, and is an artifact of the Reynolds averaging procedure. Specifically, the Reynolds averaging of the fluid acceleration yields $\overline{v_{i,t}+v_j{v}_{i,j}}={\overline{v}}_{i,t}+\overline{{\tilde{v}}_j{\tilde{v}}_{i,j}}$, where the mass conservation condition $v_{i,i}=0$ allows us to express the averaged fluid acceleration 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 indicates that the emergence of $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ is independent of fluid viscosity; it merely represents the interaction of fluid velocity with its gradient. Moreover, the inviscid Euler equation, which omits the viscous term $\nu v_{i,jj}$, also yields $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ upon Reynolds averaging. As George has pointed out [18], $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ is not a stress, but rather a transformed representation of the fluctuating contributions to the nonlinear acceleration terms. Since inviscid fluids are incapable of generating turbulent motion, it is clear that $\overline{{\tilde{v}}_i{\tilde{v}}_j}$ is not the root cause of turbulence.

The enigma of turbulence generation lies at the heart of fluid mechanics. While the precise prediction of turbulence remains an elusive goal, our understanding has crystallized around the conditions under which it does not arise: in the realm of inviscid fluid flows, turbulence is absent. The current conceptualization of turbulent flow can be summarized at a macroscopic level as follows: the interplay between fluid viscosity and velocity gradients is the principal catalyst for turbulence. Viscosity gives rise to vorticity in flows characterized by non-zero velocity gradients, precipitating fluid rotation. In response to mass conservation, the fluid must compensate for the mass that is displaced. Furthermore, the nonlinear nature of the convection term introduces a intricate pattern of motion. As the velocity gradient diminishes, the flow tends towards laminarity; conversely, when the velocity gradient intensifies, the flow transitions towards a turbulent state.

In the existing framework of the Navier-Stokes equations, a single term, $\nu v_{i,kk}$, within the equation $v_{i,t}+v_j{v}_{i,j}=-p_{,i}/\rho +\nu v_{i,kk}$, encapsulates both the viscosity coefficient $\nu$ and the gradient of velocity $v_{i,k}$, specifically its divergence $v_{i,kk}$. However, a limitation arises as the Reynolds averaging process yields $\overline{{\tilde{v}}_{i,kk}}=0$, leaving the divergence $v_{i,kk}$ with only the mean component ${\overline{v}}_{i,kk}$, excluding the fluctuating parts. This suggests that the Navier-Stokes equations, as formulated with $\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$, might not capture the full spectrum of viscous fluid dynamics and are, therefore, potentially incomplete, calling for an adjustment. Since $\nu v_{i,kk}$ is derived from the linear constitutive relation of viscous fluids in conjunction with the mass conservation condition $v_{j,j}=0$, or $\nu {(v_{i,j}+v_{j,i})}_j=\nu v_{i,kk}$, it is imperative to revisit the fluid’s constitutive relationship to ensure the inclusion of the comprehensive dynamics of fluid motion.

To derive the equations that govern turbulent flow, we must return to the foundational concepts of fluid mechanics by deriving the momentum conservation equation from a "first principles" approach. The concept of "first principles" dates back to ancient Greek philosophy, where Aristotle introduced it in his seminal work, "Metaphysics." According to Aristotle, first principles are the most elementary propositions or assumptions that are not derived from anything more fundamental; they serve as the bedrock for all knowledge and logical reasoning. This approach involves delving into the core of a subject, examining its most basic principles and assumptions, to forge new cognitive structures or solutions. Elon Musk has been a prominent advocate for this term, expressing his admiration for the "first principles" thinking method in interviews: "By applying first principles thinking, I distill the essence of matters and work up from the most fundamental truths..." He emphasizes the importance of this perspective over comparative thinking, which can lead to merely incremental advancements. First principles thinking, as Musk describes it, involves analyzing the world with a physicist’s lens, peeling back the layers of appearance to understand the intrinsic nature of things, and then constructing from this essence in a systematic manner. This is Musk’s interpretation of the "first principles thinking model" - a method of tracing the origins of concepts and reevaluating conventional approaches. The Navier-Stokes equations have stood as a cornerstone of fluid dynamics since their establishment in 1822. To achieve significant breakthroughs in our understanding, we must adopt a first principles mindset, beginning our inquiry from the ground up.

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

(1)

where $\mathbf{b}$ is body force density, $\mathbf{\sigma }={\sigma }_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ is the stress tensor, ${\sigma }_{ij}$ are the components of the stress tensor, ${\mathbf{e}}_i$ are the unit base vectors in the Cartesian coordinate system, and $\nabla ={\mathbf{e}}_i{\partial }_i$ is the gradient operator and $\mathbf{\sigma }·\nabla$ is right divergence of the stress tensor $\mathbf{\sigma }$, whose computational rule is $\mathbf{\sigma }·\nabla =({\sigma }_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j)·\left({\mathbf{e}}_k{\partial }_k\right)={\sigma }_{ij,k}{\mathbf{e}}_i({\mathbf{e}}_j·{\mathbf{e}}_k)={\sigma }_{ij,k}{\mathbf{e}}_i{\delta }_{jk}={\sigma }_{ij,j}{\mathbf{e}}_i$, in general $\mathbf{\sigma }·\nabla \ne \nabla ·\mathbf{\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 $\mathbf{\sigma }$ and the rate-of-strain tensor $\mathbf{S}$, namely $\mathbf{\sigma }=\mathbf{\sigma }\left(\mathbf{S}\right)$, where $\mathbf{S}={\displaystyle \frac{1}{2}}(\nabla \mathbf{v}+\mathbf{v}\nabla )=S_{ij}{\mathbf{e}}_i\otimes {\mathbf{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 is important to recognize that the introduction of a fluid’s constitutive relation involves a degree of human interpretation in comprehending the physical attributes of the fluid. This interpretive aspect is perhaps the most subjective element in formulating the equations of fluid motion. From Navier to Stokes, the prevailing assumption has been that fluids are isotropic and incompressible, leading to a proposed linear relationship between the stress tensor $\mathbf{\sigma }$ and the rate-of-strain tensor $\mathbf{S}$ as expressed by $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{S}$, where $\mu$ denotes the viscosity coefficient. This constitutive relation was conceived in a time bereft of tensor theory and has since been universally embraced by the scientific literature. In the context of modern science, with our understanding of tensors, it is incumbent upon us to question whether the relationship between the stress tensor $\mathbf{\sigma }$ and the strain rate tensor $\mathbf{S}$ can be articulated more rationally, and if so, to determine the appropriate form for this new expression.

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

$$\begin{array}{cc}\hfill \mathbf{\sigma }& =-p\mathbf{I}+2\mu \mathbf{S}+4\lambda {\mathbf{S}}^2\hfill \\ \hfill & =-p\mathbf{I}+\mu (\nabla \mathbf{v}+\mathbf{v}\nabla )+\lambda {(\nabla \mathbf{v}+\mathbf{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]{\mathbf{e}}_i\otimes {\mathbf{e}}_j\hfill \end{array}$$

(2)

where ${\mathbf{S}}^2=\mathbf{S}·\mathbf{S}=(S_{ik}{\mathbf{e}}_i\otimes {\mathbf{e}}_k)·(S_{lj}{\mathbf{e}}_l\otimes {\mathbf{e}}_j)=S_{ik}{S}_{lj}({\mathbf{e}}_k·{\mathbf{e}}_l){\mathbf{e}}_i\otimes {\mathbf{e}}_j=S_{ik}{S}_{lj}{\delta }_{kl}{\mathbf{e}}_i\otimes {\mathbf{e}}_j=S_{il}{S}_{lj}{\mathbf{e}}_i\otimes {\mathbf{e}}_j=[(v_{i,k}+v_{k,i}){\mathbf{e}}_i\otimes {\mathbf{e}}_k]·[(v_{l,j}+v_{j,l}){\mathbf{e}}_l\otimes {\mathbf{e}}_j]=(v_{i,k}+v_{k,i})(v_{k,j}+v_{j,k}){\mathbf{e}}_i\otimes {\mathbf{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 with dimensions dim $\lambda =M/L$.

To determine the entropy inequality for the given Cauchy stress tensor $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{S}+4\lambda {\mathbf{S}}^2$, we start from the Clausius-Duhem inequality, which enforces the second law of thermodynamics. For simplicity, we consider isothermal conditions and focus on the mechanical dissipation part of the entropy inequality. The Clausius-Duhem inequality [33,34,35,36] under isothermal conditions reduces to: $\mathbf{\sigma }:\mathbf{S}-\rho \dot{\psi }\ge 0$, where $\mathbf{S}$ is the rate of deformation tensor (symmetric part of the velocity gradient), $\rho$ is the density, and $\psi$ is the Helmholtz free energy, in which $\mathbf{\sigma }:\mathbf{B}$ is double dot product of two tensor $\mathbf{\sigma }={\sigma }_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ and $\mathbf{S}=S_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ with computational rule: $\mathbf{\sigma }:\mathbf{S}=({\sigma }_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j):(S_{kl}{\mathbf{e}}_k\otimes {\mathbf{e}}_l)={\sigma }_{ij}{S}_{kl}({\mathbf{e}}_i·{\mathbf{e}}_k)\otimes ({\mathbf{e}}_j·{\mathbf{e}}_l)={\sigma }_{ij}{S}_{kl}{\delta }_{ik}{\delta }_{jl}={\sigma }_{ij}{S}_{ij}$.

For a fluid where the stress tensor is given, we substitute $\mathbf{\sigma }$ into the inequality $\mathbf{\sigma }:\mathbf{S}-\rho \dot{\psi }\ge 0$. The stress tensor is: $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{S}+4\lambda {\mathbf{S}}^2$. Substituting this into the dissipation term $\mathbf{\sigma }:\mathbf{S}=(-p\mathbf{I}+2\mu \mathbf{S}+4\lambda {\mathbf{S}}^2):\mathbf{S}$. Since ${\mathbf{S}}^2:\mathbf{S}=(\mathbf{S}·\mathbf{S}):\mathbf{S}=(S_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j·S_{kl}{\mathbf{e}}_k\otimes {\mathbf{e}}_l):S_{mn}{\mathbf{e}}_m\otimes {\mathbf{e}}_n=S_{ij}{S}_{kl}{S}_{mn}({\mathbf{e}}_i\otimes {\mathbf{e}}_j·{\mathbf{e}}_k\otimes {\mathbf{e}}_l):({\mathbf{e}}_m\otimes {\mathbf{e}}_n)=S_{ij}{S}_{kl}{S}_{mn}{\delta }_{jk}({\mathbf{e}}_i\otimes {\mathbf{e}}_l):({\mathbf{e}}_m\otimes {\mathbf{e}}_n)=S_{ik}{S}_{kl}{S}_{mn}({\mathbf{e}}_i·{\mathbf{e}}_m)({\mathbf{e}}_l·{\mathbf{e}}_n)=S_{ik}{S}_{kl}{S}_{mn}{\delta }_{im}{\delta }_{ln}=S_{ik}{S}_{kl}{S}_{il}$, and $\mathrm{tr}{\mathbf{S}}^3={\mathbf{S}}^3:\mathbf{I}=(\mathbf{S}·\mathbf{S}·\mathbf{S}):\mathbf{I}=(S_{ik}{S}_{kl}{S}_{ln}{\mathbf{e}}_i·{\mathbf{e}}_n):({\mathbf{e}}_p·{\mathbf{e}}_p)=S_{ik}{S}_{kl}{S}_{ln}({\mathbf{e}}_i·{\mathbf{e}}_p)({\mathbf{e}}_n·{\mathbf{e}}_p)=S_{ik}{S}_{kl}{S}_{li}$, noting the symmetric tensor $S_{il}=S_{li}$, leads to ${\mathbf{S}}^2:\mathbf{S}=\mathrm{tr}{\mathbf{S}}^3$. Thus, the dissipation term becomes: $-p\mathrm{tr}\left(\mathbf{S}\right)+2\mu \mathbf{S}:\mathbf{S}+4\lambda \mathrm{tr}\left({\mathbf{S}}^3\right)$. For an incompressible fluid, $\mathrm{tr}\left(\mathbf{S}\right)=0$, simplifying the dissipation term to: $\mathbf{\sigma }:\mathbf{S}=(-p\mathbf{I}+2\mu \mathbf{S}+4\lambda {\mathbf{S}}^2):\mathbf{S}=2\mu \mathbf{S}:\mathbf{S}+4\lambda \mathrm{tr}\left({\mathbf{S}}^3\right)$.

The entropy inequality requires this term to be non-negative for all possible $\mathbf{S}$. Therefore, the entropy inequality for the given stress tensor is: $2\mu \mathbf{S}:\mathbf{S}+4\lambda \mathrm{tr}\left({\mathbf{S}}^3\right)\ge 0$, which satisfy the second law of thermodynamics for all possible flow motion [33,34,35,36,33].

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 \mathbf{S}+4\lambda {\mathbf{S}}^2):\mathbf{S}{d}^3\mathbf{x}=-{\int }_{\Omega }\left[2\mu \mathbf{S}:\mathbf{S}+4\lambda \mathrm{tr}({\mathbf{S}}^3\right]d^3\mathbf{x},\hfill \end{array}$$

(3)

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

In the annals of fluid mechanics, from the works of Navier and Stokes to contemporary research, fluids have consistently been treated as Newtonian, employing the linear constitutive relation $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{S}$, while the second-order term $\lambda {\mathbf{S}}^2$ has been omitted. The exclusion of the term $\lambda {\mathbf{S}}^2$ can be attributed to two primary factors. Firstly, during the era of Newton, Navier, and Stokes, nonlinear effects were not detectable, and the necessary tools for tensor analysis, including tensor representation theory and the Cayley-Hamilton theorem[37], were not yet developed; indeed, the Cayley-Hamilton theorem was published after the formulation of the Navier-Stokes equations. Secondly, even if the expression Equation 2 were to be derived, it was generally assumed that the second-order term $\lambda {\mathbf{S}}^2$ was sufficiently insignificant to be disregarded. That is to say, the second-order constitutive relation, $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{S}+4\lambda {\mathbf{S}}^2$, is the correct and complete result, whereas its first-order constitutive relation is actually an approximation. This was not understood in the era before the development of tensor representation theory.

This study contends that the omission of the second-order term $\lambda {\mathbf{S}}^2$ in the formulation of fluid mechanics equations is a critical oversight. Given that the spatial variation of velocity, such as the gradient near a wall, can be extremely pronounced, neither the first nor the second power of the velocity gradient should be dismissed. This resonates with the scenario encountered in the boundary layer theory. Prandtl [38] emphasized, upon introducing the boundary layer concept, that the Reynolds number $Re$ in the term ${\displaystyle \frac{1}{Re}}{\nabla }^2\mathbf{v}$ must not be ignored, regardless of how high the value of $Re$ (or how low ${\displaystyle \frac{1}{Re}}$) may be, because the velocity gradient is sufficiently steep to render the entire term ${\displaystyle \frac{1}{Re}}{\nabla }^2\mathbf{v}$ significant.

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

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

(4)

wherein, the kinematic viscosity is $\nu =\mu /\rho$, and the kinematic second-order viscosity is $\kappa =\lambda /\rho$ with dimensions dim $\kappa =L^2$. 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 ·\mathbf{v}=0$, Equation 4 can be simplified to:

$${\displaystyle \frac{\partial \mathbf{v}}{\partial t}}+\mathbf{v}·(\nabla \mathbf{v})=\mathbf{b}-{\displaystyle \frac{1}{\rho }}\nabla p+\nu {\nabla }^2\mathbf{v}+\kappa \left[{(\nabla \mathbf{v}+\mathbf{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 \mathbf{v}+\mathbf{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}& =b_i-{\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}=& b_x-{\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}=& b_y-{\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}=& b_z-{\displaystyle \frac{1}{\rho }}{p}_{,z}+\nu {\nabla }^{2}w+\kappa M_z,\hfill \end{array}$$

(9)

where $M_x={[4{\left(u_{,x}\right)}^2+{(v_{,x}+u_{,y})}^2+{(w_{,x}+u_{,z})}^2]}_{,x}+{[2(u_{,x}+v_{,y})(u_{,y}+v_{,x})+(w_{,y}+v_{,z})(w_{,x}+u_{,z})]}_{,y}+{[2(u_{,x}+w_{,z})(u_{,z}+w_{,x})+(v_{,z}+w_{,y})(v_{,x}+u_{,y})]}_{,z}$, $M_y={[2(u_{,x}+v_{,y})(u_{,y}+v_{,x})+(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})+2(v_{,y}+w_{,z})(v_{,z}+w_{,y})]}_{,z}$ and $M_z={[2(u_{,x}+w_{,z})(u_{,z}+w_{,x})+(v_{,x}+u_{,y})(v_{,z}+w_{,y})]}_{,x}+{[(u_{,y}+v_{,x})(u_{,z}+w_{,x})+2(v_{,y}+w_{,z})(v_{,z}+w_{,y})]}_{,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}=& b_x-{\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}=& b_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 the two dimensional flow, we have $\mathbf{S}=v_{1,1}{\mathbf{e}}_1\otimes {\mathbf{e}}_1+{\displaystyle \frac{1}{2}}(v_{1,2}+v_{2,1}){\mathbf{e}}_1\otimes {\mathbf{e}}_2+{\displaystyle \frac{1}{2}}(v_{2,1}+v_{1,2}){\mathbf{e}}_2\otimes {\mathbf{e}}_1+v_{2,2}{\mathbf{e}}_2\otimes {\mathbf{e}}_2$ and ${\mathbf{S}}^2=\mathbf{S}·\mathbf{S}=[{\left(v_{1,1}\right)}^2+{\displaystyle \frac{1}{4}}{(v_{1,2}+v_{2,1})}^2]{\mathbf{e}}_1\otimes {\mathbf{e}}_1+{\displaystyle \frac{1}{2}}(v_{1,1}+v_{2,2})(v_{1,2}+v_{2,1}){\mathbf{e}}_1\otimes {\mathbf{e}}_2+{\displaystyle \frac{1}{2}}(v_{1,1}+v_{2,2})(v_{1,2}+v_{2,1}){\mathbf{e}}_2\otimes {\mathbf{e}}_1+[{\left(v_{2,2}\right)}^2+{\displaystyle \frac{1}{4}}{(v_{1,2}+v_{2,1})}^2]{\mathbf{e}}_2\otimes {\mathbf{e}}_2$. Noice continuity condition $v_{1,1}+v_{2,2}=0$, results to ${\mathbf{S}}^2=[{\left(v_{1,1}\right)}^2+{\displaystyle \frac{1}{4}}{(v_{1,2}+v_{2,1})}^2]{\mathbf{e}}_1\otimes {\mathbf{e}}_1+[{\left(v_{2,2}\right)}^2+{\displaystyle \frac{1}{4}}{(v_{1,2}+v_{2,1})}^2]{\mathbf{e}}_2\otimes {\mathbf{e}}_2$, thus, we have in conventional forms $\mathbf{S}:\mathbf{S}=u_x^2+v_{,y}^2+{\displaystyle \frac{1}{2}}{(u_{,y}+v_{,x})}^2$ and ${\mathbf{S}}^2:\mathbf{S}=u_{,x}^3+v_{,y}^3$. The entropy inequality can be expressed as to:

$$2\mu \mathbf{S}:\mathbf{S}+4\lambda {\mathbf{S}}^2:\mathbf{S}=2\mu \left[u_x^2+v_{,y}^2+{\displaystyle \frac{1}{2}}{(u_{,y}+v_{,x})}^2+{\displaystyle \frac{2\kappa }{\nu }}(u_{,x}^3+v_{,y}^3)\right]\ge 0.$$

(12)

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

(13)

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

3. Boundary Layer Theory and Application to Wedge Flow

In 1904, at the third International Congress of Mathematicians, Ludwig Prandtl [38] unveiled the notion of the boundary layer, proposing that the effects of frictional forces are primarily restricted to a slender layer adjacent to an object’s surface. Inside this boundary layer, both the velocity gradient and the shear stress are substantial, generating a frictional resistance that must be accounted for. The equations governing the boundary layer are notably less complex than the full Navier-Stokes equations and can be addressed numerically[39]. Notably, Prandtl identified that the flow within the boundary layer can exhibit turbulence[40,41], and in 1925, he introduced the concept of the mixing length to quantify the turbulent viscosity[42]. This concept was instrumental in formulating the logarithmic law for the wall in turbulent boundary layers and represented a pivotal effort to close the Reynolds-Averaged Navier-Stokes (RANS) equations.

In 1908, Prandtl’s student Blasius[39] 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 as shown in Figure 1, 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[39] 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[43] 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}$$

(14)

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

(15)

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 Equations 14 and 15 can be easily solved by using Maple, the solutions of velocity components $u,$ and $v$ are mainly expressed in terms of Lambert function[45]: LambertW=$W\left(z\right)$, which satisfies: $W\left(z\right)exp\left(W\right(z\left)\right)=z$. Although we can obtain the solution, it is regrettable that the integral involving the Lambert function 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 is satisfied, and the momentum conservation Equation 14 becomes

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

(16)

and corresponding boundary conditions.

Applying Sun’s similarity transformations [43]: $\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}$$

(17)

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

(18)

Substituting Equation 17 and Equation into Equation 16, 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]+{\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}$$

(19)

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 19 explicitly contains the coordinate x, which conn’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$, we have $u=Cx{f}_{,\eta },v=-\sqrt{\nu C}f$, and the Equation 19 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}$$

(20)

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 20 has get rid of coordinate x, can be solved numerically.

For its corresponding steady problem, the function f is only a function of $\eta$, the Equation 20 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}$$

(21)

and boundary conditions: $\eta =0:f=f_{,\eta }=0$ and $\eta \to \infty :f_{,\eta }=1$. Equation 21 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$.

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 2. 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(\underset{\underset{=0}{︸}}{u_{,x}+v_{,y}}\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 }_{\mathrm{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}=\rho \nu {\displaystyle \frac{Cx}{\sqrt{\nu /C}}}{f}_{,\eta \eta }\left(0\right)$:

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

(22)

This expression reveals that the resistance of the flow with $\kappa$ is much greater than that of the flow with $\kappa =0$. 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)$.

Up to this point, there is still one issue that has not been resolved, which is to verify the entropy inequality. If we apply Equation 12 to the wedge flow, noting the relationships $u=Cx{f}_{,\eta },v=-\sqrt{\nu C}f,u_{,x}=C{f}_{,\eta },u_{,y}={\displaystyle \frac{C^3}{\nu }}x{f}_{,\eta \eta },v_{,x}=0,v_{,y}=-C{f}_{,\eta }$, the entropy inequality can be expressed as follows:

$$2\mu \left[2C^2{f}_{,\eta }^2+{\displaystyle \frac{1}{2}}{\displaystyle \frac{C^3}{\nu }}{x}^2{\left(f_{,\eta \eta }\right)}^2\right]\ge 0,$$

(23)

which holds if $C>0$.

4. Steady Incompressible Flow in an Inclined Pipe and Determination of the Second Viscosity

After establishing new equations for the incompressible fluids, one remaining question is how to determine the second-order viscosity coefficient $\lambda$. Below, let’s demonstrate through an example that, when the dynamic viscosity coefficient $\mu$ of the fluid is known, the second-order viscosity coefficient $\lambda$ can be obtained using a simple experiment.

Consider steady incompressible flow in an infinitely long round pipe with a radius $R={\displaystyle \frac{D}{2}}$ inclined at an angle $\alpha$ (Figure 3). The fluid flows down the pipe due to gravity, meaning no applied pressure gradient is present in the fluid.

We use the semi-inverse method, which involves first determining the general form of the flow velocity based on the problem’s conditions, then substituting this into various expressions to see under what boundary conditions the initially assumed flow velocity holds. Due to the pipe’s axial symmetry, we can assume the flow is also axisymmetric, and thus the velocity field can be expressed as a function of radius r only, $u=0$, $v=0$, and $w=w\left(r\right)$. The no-slip condition at the wall $r=R$ is $u=v=w=0$, and there is pressure balance at the two ends of the pipe: ${\sigma }_{ij}{n}_j{|}_{x=0}=p_i^0$ and ${\sigma }_{ij}{n}_j{|}_{x=L}=p_i^L$.

Using the conversion relationships between Cartesian and polar coordinates, $r=\sqrt{x^2+y^2}$, $x=rcos\theta$, $y=rsin\theta$, and substituting $v_1=u=0$ and $v_2=v=0$, $v_3=w=w\left(r\right)$ into the strain-rate tensor $\mathbf{S}={\displaystyle \frac{1}{2}}(\nabla \mathbf{v}+\mathbf{v}\nabla )$, we obtain $\mathbf{S}={\displaystyle \frac{1}{2r}}{\displaystyle \frac{dw}{dr}}(x{\mathbf{e}}_x\otimes {\mathbf{e}}_z+y{\mathbf{e}}_y\otimes {\mathbf{e}}_z+x{\mathbf{e}}_z\otimes {\mathbf{e}}_x+y{\mathbf{e}}_z\otimes {\mathbf{e}}_y)$. Using the differential relations $r_{,1}={\displaystyle \frac{x}{r}}$, $r_{,2}={\displaystyle \frac{y}{r}}$, this constitutive relation can be expressed in matrix form as follows:

$$\left[\sigma \right]=-p\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)+{\displaystyle \frac{\mu }{r}}{\displaystyle \frac{dw}{dr}}\left(\begin{array}{ccc}0& 0& x\\ 0& 0& y\\ x& y& 0\end{array}\right)+{\displaystyle \frac{\lambda }{r^2}}{\left({\displaystyle \frac{dw}{dr}}\right)}^2\left(\begin{array}{ccc}{x}^2& xy& 0\\ xy& y^2& 0\\ 0& 0& x^2+y^2\end{array}\right).$$

(24)

For steady flow, $v_1=u=0$ and $v_2=v=0$, $v_3=w=w\left(r\right)$, the continuity condition is automatically satisfied. The momentum equations are as follows:

$$\begin{array}{cc}\hfill -\rho gcos\alpha -p_{,x}+\lambda {\left({\displaystyle \frac{x{w}_{,r}}{r}}\right)}_{,x}^2+\lambda {\left({\displaystyle \frac{xy{w}_{,r}^2}{r^2}}\right)}_{,y}& =0,\hfill \\ \hfill -p_{,y}+\lambda {\left({\displaystyle \frac{xy{w}_{,r}^2}{r^2}}\right)}_{,x}+\lambda {\left({\displaystyle \frac{y{w}_{,r}}{r}}\right)}_{,y}^2& =0,\hfill \\ \hfill \rho gsin\alpha -p_{,z}+\mu {\left({\displaystyle \frac{x{w}_{,r}}{r}}\right)}_{,x}+\mu {\left({\displaystyle \frac{y{w}_{,r}}{r}}\right)}_{,y}& =0,\hfill \end{array}$$

(25)

Using the first and second expressions of Equation 25 and noting the chain rule of calculus, we get $p_{,z}=-\rho gsin\alpha +{\displaystyle \frac{\mu }{r}}{\left(r{w}_{,r}\right)}_{,r}$, since the right-hand side does not depend on z, the left-hand side must be a constant c (which is actually the pressure gradient along the z-direction), thus $p=cz+b(r,\theta )$, where $b(r,\theta )$ is an integration function. Substituting this into Equation 25 gives the third expression, leading to the equation $c=-\rho gsin\alpha +{\displaystyle \frac{\mu }{r}}{\left(r{w}_{,r}\right)}_{,r}$, and considering the boundary conditions $w\left(0\right)=0$, $w\left(R\right)=0$, the solution can be written as

$$w=-{\displaystyle \frac{(\rho gsin\alpha +c)}{4\mu }}(R^2-r^2).$$

(26)

Bringing $p=cz+b(r,\theta )$ back into $p_{,\theta }=p_{,x}{x}_{,\theta }+p_{,y}{y}_{,\theta }=p_{,x}(-y)+p_{,y}x=-\rho gcos\alpha sin\theta$, we get $b=r\rho gcos\alpha cos\theta +h\left(r\right)$, where $h\left(r\right)$ is an integration function. Thus, the pressure is $p=cz+r\rho gcos\alpha cos\theta +h\left(r\right)$. Substituting this into $p_{,r}=p_{,x}{x}_{,r}+p_{,y}{y}_{,r}=p_{,x}{\displaystyle \frac{x}{r}}+p_{,y}{\displaystyle \frac{y}{r}}=\rho gcos\alpha cos\theta +{\displaystyle \frac{\lambda }{r}}{\left[r{\left(w_{,r}\right)}^2\right]}_{,r}$, we obtain ${\displaystyle \frac{dh}{dr}}={\displaystyle \frac{\lambda }{r}}{\left[r{\left(w_{,r}\right)}^2\right]}_{,r}$, and after integrating with Equation 26, we get $h={\displaystyle \frac{3\lambda {(sin\alpha g\rho +c)}^2{r}^2}{8{\mu }^2}}+d\left(\theta \right)$, where d is an integration function. Therefore, the pressure can be expressed as

$$p=cz+{\displaystyle \frac{3\lambda {(sin\alpha g\rho +c)}^2{r}^2}{8{\mu }^2}}+\rho grcos\alpha cos\theta +d.$$

(27)

The Cauchy stress can then be further expressed as

$$\left[\sigma \right]=-p\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)+{\displaystyle \frac{1}{2}}(c+\rho gsin\alpha )\left(\begin{array}{ccc}0& 0& x\\ 0& 0& y\\ x& y& 0\end{array}\right)+{\displaystyle \frac{{(\rho gsin\alpha +c)}^2\lambda }{4{\mu }^2}}\left(\begin{array}{ccc}{x}^2& xy& 0\\ xy& y^2& 0\\ 0& 0& x^2+y^2\end{array}\right).$$

(28)

Now, apply the boundary conditions at $z=0$ and $z=L$ to determine the undetermined constants c and d. On the cross-section at $z=0$, the shear stresses are ${\sigma }_{zr}={\sigma }_{31}=-p_1^0$ and ${\sigma }_{z\theta }={\sigma }_{32}=-p_2^0$, and the normal stress is ${\sigma }_{zz}={\sigma }_{33}=-p_z^0$. Similarly, on the cross-section at $z=L$, the shear stresses are ${\sigma }_{zr}={\sigma }_{31}=p_1^L$ and ${\sigma }_{z\theta }={\sigma }_{32}=p_2^L$, and the normal stress is ${\sigma }_{zz}=p_z^L$. Applying these conditions to p and $\sigma$ yields the boundary conditions: ${\displaystyle \frac{1}{2}}(c+\rho gsin\alpha )x=-p_1^0$, ${\displaystyle \frac{1}{2}}(c+\rho gsin\alpha )y=-p_2^0$, $-(\rho gcos\alpha cos\theta r+d)-{\displaystyle \frac{{(\rho gsin\alpha +c)}^2\lambda }{8{\mu }^2}}{r}^2=-p_z^0$, ${\displaystyle \frac{1}{2}}(c+\rho gsin\alpha )x=p_1^L$, ${\displaystyle \frac{1}{2}}(c+\rho gsin\alpha )y=p_2^L$, $cL-(\rho gcos\alpha cos\theta r+d)-{\displaystyle \frac{{(\rho gsin\alpha +c)}^2\lambda }{8{\mu }^2}}{r}^2=p_z^L$. If $p_z^0=p_z^L$, then it must be the case that $c=0$, which corresponds to the situation with no pressure gradient. Therefore, the boundary stress conditions are as follows: $p_1^0=p_x^0=-{\displaystyle \frac{1}{2}}\rho gsin\alpha x={\displaystyle \frac{1}{2}}\rho grsin\alpha cos\theta$, $p_2^0=p_y^0=-{\displaystyle \frac{1}{2}}\rho gsin\alpha y={\displaystyle \frac{1}{2}}\rho grsin\alpha sin\theta$ and $p_z^0={\displaystyle \frac{{(\rho gsin\alpha )}^2\lambda }{8{\mu }^2}}{r}^2+\rho grcos\alpha cos\theta +d$, as well as $p_1^L=p_x^L={\displaystyle \frac{1}{2}}\rho gsin\alpha x={\displaystyle \frac{1}{2}}\rho grsin\alpha cos\theta$, $p_2^L=p_y^L={\displaystyle \frac{1}{2}}\rho gsin\alpha y={\displaystyle \frac{1}{2}}\rho grsin\alpha sin\theta$ and $p_z^L=-{\displaystyle \frac{{(\rho gsin\alpha )}^2\lambda }{8{\mu }^2}}{r}^2-\rho grcos\alpha cos\theta -d$. This means that, in order to obtain the velocity form $v_1=u=0$, $v_2=v=0$ and $v_3=w=w\left(r\right)$ assumed in advance by the semi-inverse method, the above pressure boundary conditions must be provided, with the resultant force and resultant bending moment being zero, such as the force over its area of application $\int p_x^{0}rdrd\theta =0$. If we set the pressure at $r=R$ to be equal to $P_R$, i.e., ${p|}_{r=R}=P_R$, and considering there is no pressure gradient $c=0$, we can determine $d=P_R-{\displaystyle \frac{3\lambda {(\rho gsin\alpha )}^2}{8{\mu }^2}}{R}^2-\rho gRcos\alpha cos\theta$, thus obtaining the following expression for pressure:

$$P_R-p={\displaystyle \frac{3\lambda {(\rho gsin\alpha )}^2}{8{\mu }^2}}(R^2-r^2)+\rho g(R-r)cos\alpha cos\theta .$$

(29)

Since this expression includes the dynamic viscosity coefficient $\mu$ and the second-order viscosity coefficient $\lambda$, and given that $\nu$ is known, this expression can guide experiments to determine the second-order viscosity coefficient $\lambda$. For simplicity, if we place the circular pipe vertically (as in Figure 3), at this point $\alpha ={\displaystyle \frac{\pi }{2}}$, the pressure expression becomes simpler: $P_R-p={\displaystyle \frac{3\lambda {\left(\rho g\right)}^2}{8{\mu }^2}}(R^2-r^2)$. If we can measure the pressure at the axis of the circular pipe $P_0{=p|}_{r=0}$, we can obtain the following expression to predict the second-order viscosity coefficient $\lambda$:

$$\lambda ={\displaystyle \frac{8\pi }{3}}{\left({\displaystyle \frac{\mu }{\rho g}}\right)}^2{\displaystyle \frac{P_R-P_0}{\pi R^2}}.$$

(30)

With the above relation, we can determine the second-order viscosity by the fluid flowing down the pipe due to gravity alone. Especially from the structure of this expression in Equation 30, it can be seen that for a fluid with known $\mu$, its corresponding second-order viscosity coefficient $\lambda$ is directly proportional to the pressure difference between the pipe wall and its center on a unit area.

5. Discussions and Conclusions

In the light of the theoretical developments and practical applications presented herein, the quadratic term $\kappa {\mathbf{S}}^2$ emerges as a pivotal element in the depiction of turbulent flow dynamics. It is arguable that the Navier-Stokes equations, devoid of the $\kappa {\mathbf{S}}^2$ term, are fundamentally tailored to laminar flow phenomena and are insufficient for capturing the intricacies of turbulent motion. The pronounced velocity gradients inherent in turbulence necessitate the retention of the $\kappa {\mathbf{S}}^2$ term; its omission would preclude an accurate representation of turbulent behavior in the derived equations.

While the focus of this work is on the turbulence regime, the equations Equation 5 and Equation 6 that have been formulated are, in fact, universally applicable to the motion of viscous fluids. The derivation of these equations is underpinned by logical consistency at each stage, eschewing any arbitrary approximations. These equations hold potential for a wide range of incompressible viscous fluid flows, including aerodynamic gas dynamics and the acoustic emissions from turbulent processes. The same conceptual framework employed for Equation 2 can be extended to establish motion equations for compressible viscous fluids as shown in the Appendix.