Derivation of Navier-Stokes Equation in Cartesian Coordinates

Problems raised in the differenential momentum equation can be solved by introducing relationships between stress components and equations between stress components and velocity gradients. Since the velocity gradient in flow field has the meaning of strain rate, the relationship between velocity gradient and stress is the same as the relationship between strain and stress in solid mechanics. In solid mechanics, this relationship is called Hooke's law, and according to it, stress is proportional to the strain of an elastic material. Without these laws, problems in solid mechanics cannot be solved. In electrical engineering, there is a similar law, i.e. Ohm's law, and in heat transfer, Fourier's law plays a similar role. These relationships between solid and fluid mechanics contain knowledge about how an object will behave after being acted on by a force.

Let us derive equations between the velocity gradient and stress components discovered by Stokes in 1845. First, the constitutive model of viscous fluid flow defines the relationship between the shear strain rate D' and the deviatoric stress σ'. γ=2D' hense,

\[\boldsymbol\sigma'=2\mu{\bf D}'\]

where μ is the viscosity coefficient of the fluid. Since the strain rate tensor D=(1/3)tr(D)+D'=(1/3)∇·V+D' and the stress tensor σ=-pI+σ', substitute them into the above equation and organize it:

\[\boldsymbol\sigma=-p{\bf I}+2\mu\left\{{\bf D}-{1\over3}\text{tr}(\bf D){\bf I}\right\}=-p{\bf I}+2\mu{\bf D}-{2\over3}\mu\nabla\cdot{\bf VI}\cdots(1)\]

where p is the hydrostatic pressure, tr(D) is the diagonal sum of the strain rate tensor D, V is the velocity vector, ∇ is the del operator, and I is the unit matrix. In order to express the above equation as components, we express the stress, strain rate and unit matrix in matrix notation as shown below.

\[\begin{align}\boldsymbol\sigma&=\begin{bmatrix}\sigma_{xx}&\tau_{xy}&\tau_{xz}\\\tau_{yx}&\sigma_{yy}&\tau_{yz}\\\tau_{zx}&\tau_{zy}&\sigma_{zz}\end{bmatrix}\\{\bf D}&=\begin{bmatrix}\frac{\partial u}{\partial x}&{1\over2}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)&{1\over2}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\\{1\over2}\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)&\frac{\partial v}{\partial y}&{1\over2}\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\\{1\over2}\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)&{1\over2}\left(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}\right)&\frac{\partial w}{\partial z}\end{bmatrix}\\{\bf I}&=\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}\end{align}\]

The u, v and w in strain rate tensor are the x, y and z component in velocity vector V respectively.

Substituting these matrices into equation (1) and comparing eash component, the relationship between stresses and velocity gradients is as follows. However, these expressions can only be applied to Newtonian fluids.

\[\begin{align}\sigma_{xx}&=-p+2\mu\frac{\partial u}{\partial x}-{2\over3}\mu\nabla\cdot{\bf V}\\\sigma_{yy}&=-p+2\mu\frac{\partial v}{\partial y}-{2\over3}\mu\nabla\cdot{\bf V}\\\sigma_{zz}&=-p+2\mu\frac{\partial w}{\partial z}-{2\over3}\mu\nabla\cdot{\bf V}\\\tau_{xy}&=\tau_{yx}=\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\\\tau_{yz}&=\tau_{zy}=\mu\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\\\tau_{zx}&=\tau_{xz}=\mu\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)\end{align}\]

The above equations for shear stresses can be seen as a generalization of τ=μdu/dy. Combining the above equation with the differential momentum equations (4), we obtain the following equations of motion.

\[\small\begin{split}\rho\frac{Du}{Dt}&=-\frac{\partial p}{\partial x}+\frac{\partial}{\partial x}\left\{\mu\left(2\frac{\partial u}{\partial x}-{2\over3}\nabla\cdot{\bf V}\right)\right\}+\frac{\partial}{\partial y}\left\{\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right\}+\frac{\partial}{\partial z}\left\{\mu\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)\right\}+\rho f_x\\\rho\frac{Dv}{Dt}&=-\frac{\partial p}{\partial y}+\frac{\partial}{\partial y}\left\{\mu\left(2\frac{\partial v}{\partial y}-{2\over3}\nabla\cdot{\bf V}\right)\right\}+\frac{\partial}{\partial z}\left\{\mu\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\right\}+\frac{\partial}{\partial x}\left\{\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right\}+\rho f_y\\\rho\frac{Dw}{Dt}&=-\frac{\partial p}{\partial z}+\frac{\partial}{\partial z}\left\{\mu\left(2\frac{\partial w}{\partial z}-{2\over3}\nabla\cdot{\bf V}\right)\right\}+\frac{\partial}{\partial x}\left\{\mu\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)\right\}+\frac{\partial}{\partial y}\left\{\mu\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\right\}+\rho f_z\end{split}\]

The equations of motion in these three directions are called the Navier-Stokes equations. The Navier-Stokes equations is an equation that opened a new chapter in fluid mechanics by being published by Navier and Stokes in 1827 and 1845, repectively, after a period of stagnation in fluid mechanics. This equation is very important in modern fluid mechanics, and almost all mathematical solutions to flow are based on this equation.

The above equations are very general and can also be applied to compressible flow accompanied by a change in viscosity coefficient. If the viscosity coefficient is assumed to be constant, the Navier-Stokes equations is simplified as follows by substituting \(\nabla\cdot{\bf V}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\) and rewrite the terms.

\[\begin{align}\frac{Du}{Dt}&=-\frac{\partial p}{\rho\partial x}+\nu\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}\right)+\frac{\nu\partial}{3\partial x}(\nabla\cdot{\bf V})+f_x\\\frac{Dv}{Dt}&=-\frac{\partial p}{\rho\partial y}+\nu\left(\frac{\partial^2v}{\partial x^2}+\frac{\partial v}{\partial y^2}+\frac{\partial^2v}{\partial z^2}\right)+\frac{\nu\partial}{3\partial y}(\nabla\cdot{\bf V})+f_y\\\frac{Dw}{Dt}&=-\frac{\partial p}{\rho\partial z}+\nu\left(\frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}+\frac{\partial^2w}{\partial z^2}\right)+\frac{\nu\partial}{3\partial z}(\nabla\cdot V)+f_z\end{align}\]

where ν=μ/ρ is the kinematic viscosity. Using vector notation, the above equations can be expressed as a single equation as follows.

\[\frac{D{\bf V}}{Dt}=-\frac{\nabla p}{\rho}+\nu\nabla^2{\bf V}+{\nu\over3}\nabla(\nabla\cdot{\bf V})+{\bf f}\cdots(2)\]

where DV/Dt has the Eulerian exression of acceleration, a vector expression given by \({\bf a}=({\bf V}\cdot\nabla){\bf V}+\frac{\partial{\bf V}}{\partial t}\). The terms on the right side represent the effects of pressure, shear force, expansion force due to compressibility, and gravity, respectively.

In the case of incompressible flow, the continuity equation, ∇·V=0 is established, so the equation (2) can be reduced as follows.

\[\frac{D{\bf V}}{Dt}=-\frac{\nabla p}{\rho}+\nu\nabla^2{\bf V}+{\bf f}\cdots(3)\]

In the equations of motion in the three directions, the unknowns are the velocity components in each direcition(u, v, w or \(v_r,\,v_\theta,\,v_z\) in cylindrical coordinates), pressure p, and density ρ, but in the case of incompressible flow, a constant value ρ can be known in advance, hence the unknowns are a total of four. However, since three equations of motion and one equation of continuity can be used to analyze flow, incompressible flow can be analyzed using only the equations of motion and continuity equations. To date, only laminar flow can be solved directly from the Navier-Stokes equations.

Boundary conditions are needed to find solutions of given differential equations. Therefore, let us learn about two important boundary conditions that are generally applied to fluid flow.

In cases where part or all of the of the boundary surface of the flow field consists of a solid, i.e. the flows such as those in Figure 1, it can be assumed that the velocity distribution is as shown in the Figure. In each case, the fluid particles located on a solid surface have no velocity relative to the solid surface. This is the boundary condition most often used to obtain flow velocity distribution (no slip condition). In the case of Figure 1(b), it can be written as follows that this condition, the flow velocity on a stationary solid surface, is also zero.

\[v_z(R,\theta,z)=0\]

Another boundary condition is that when two fluids flow in contact, the distribution of shear stress and velocity is continuous at the interface.

Figure 1

[Example] An incompressible fluid with a constant viscosity coefficient flows between two parallel plates(Figure 2). One plate is stationary and the other is moving with velocity Vp. It is assumed that there is a pressure gradient in the x-direction within the flow field, and that there are no velocity components in the y and z directions.

Figure 2

[Solution] The given flow is a steady flow, and there are no velocity components in the y and z directions, hence

\[\frac{\partial u}{\partial t}=0\ \text{and}\ v=w=0\]

Additionally, the following equation is obtained from the continuity equation of incompressible flow.

\[\frac{\partial u}{\partial x}=0\]

Since the given flow is a two-dimensional flow with no change in the z direction, all rates in z direction become zero.

Applying the above conditions to the equation of motion (3) of incompressible flow, they are simplified as below.

\[\begin{align}0&=-\frac{\partial p}{\partial x}+\mu\frac{\partial^2u}{\partial y^2}&\cdots(a)\\0&=-\frac{\partial p}{\partial y}&\cdots(b)\\0&=-\frac{\partial p}{\partial z}&\cdots(c)\end{align}\]

Equations (b) and (c) simply tell us that pressure is a function of x only. Therefore, if we rewire only equation (a), it is:

\[\frac{d^2u}{dy^2}=\frac{dp}{\mu dx}\qquad\cdots(d)\]

If you repeat integration,

\[\frac{du}{dy}=\frac{dp}{\mu dx}y+c_1\]

and,

\[u=\frac{dp}{2\mu dx}y^2+c_1+c_2\qquad\cdots(e)\]

where the pressure gradient dp/dx is assumed to be constant, and c1 and c2 are constants. To determine  these constants, the boundary conditions along the two plates are as below.

\[\begin{split}u&=V_p\ \text{at}\ y=\ \ \ a\\u&=0\ \ \ \text{at}\ y=-a\end{split}\]

Substituting these conditions into the equation (e), we obtain:

\[c_1=\frac{V_p}{2a},\qquad c_2=\frac{V_p}{2}-\frac{dp}{2\mu dx}a^2\]

Therefore, the velocity distribution equation can be obtained from the equation (e).

\[u={V_p\over2}\left(1+{y\over a}\right)-\frac{a^2dp}{2\mu dx}\left(1-{y^2\over a^2}\right)\qquad\cdots(f)\]

Figure 3 shows the velocity distribution equation (f) depending on the extent of pressure gradient.

Figure 3