Navier Stoke Equation

Geodynamics

(Lecture 8 Part 1)

Navier-Stoke Equation

Navier-Stokes equation, in fluid mechanics, is a set of two partial differential equations that describes the flow of incompressible fluids.

(248)\[\rho\left[\frac{\partial \vec{v}}{\partial t} + \vec{v}\cdot \nabla \vec{v}\right] = \rho g - \nabla p +\mu \nabla^{2} \vec{v}\]

and

(249)\[\rho (\nabla \cdot\vec{v}) = 0\]

where \(\rho\) is density and \(p\) is pressure (scalars); \(g\) is gravity and \(v\) is velocity (vectors).

The first equation ((248)) states that the motion of fluid is predominantly governed by gravity, pressure and viscosity. The second equation ((249)) represents the mass continuity.

Equation (248) and (249) is governed by and derived from the Newton’s second law of motion (250) and Law of Conservation of Mass ((251)) respectively.

(250)\[\sum \vec{F}=m\vec{a}=\frac {m\vec{v}}{t}\]

(251)\[\dot{m}_{in} = \dot{m}_{out}\]

The mass continuity equation is easy to obtain by equating all the mass flow rates into and out of the differential control volume. Hence, the following derivation is focused solely on the governing equation ((248)). We can start from the Newton’s second law of motion:

(252)\[\sum \vec{F}= \vec{F_{p}}+\vec{F_{v}}+\vec{F_{g}} =m\vec{a}\]

where \(\vec{F_{p}}\) is force due to pressure and \(\vec{F_{v}}\) is force due to viscosity; they are known as the internal forces. There are many external forces due to gravity, electromagnetism and etc. However, we consider the gravity \(\vec{F_{g}}\) as the only external force for simplicity.

The equation is developed by summing up the forces of an infinitesimal control volume \(V\) (\(V=d_{x}d_{y}d_{z}\)) in Cartesian coordinate system, and setting the result equal to \(m\vec{a}\). This is demonstrated in the figure below for viscosity forces.

Viscosity Forces

Figure 1: Infinitesimal control volume \(V\) (left) and Cartesian coordinate system (right). Examples of stresses acting on different planes and directions are also shown.

We define x-planes are the planes that are perpendicular to x-axis; As shown above, we can see that longitudinal stress \(\sigma_{xx}\) represents the stress that acts on x-plane in the direction of x-axis; shear stress \(\tau_{xy} \) represents the stress that acts on x-plane in the direction of y-axis. \(\langle\sigma_{xx} \rangle(x+dx)\) is the longitudinal stress acting on x-plane at \(x+dx\).

Due to the redundancy in 3D derivation, we shall analysis the stresses just in x-direction and conclude the resulting equation hold true for the other two directions by symmetry. Figure 2 shows the stresses in x-direction only.

Figure 2: Longitudinal stresses (\(\sigma_{xx}\)) in (a) and shear stresses (\(\tau_{xy} \)) in (b) & (c) in x-direction from six planes. Shear stresses are due to viscosity and longitudinal stresses are created by pressure.

Note that the longitudinal stresses \(\sigma_{xx}\) acting on two x-planes are not the same. We can apply Taylor Series to interpolate the \(\sigma_{xx}\) evaluated at \(x+dx\) based on the value of \(\sigma_{xx}\) at \(x\):

(253)\[\langle\sigma_{xx} \rangle(x+dx) = \sigma_{xx}(x) + \frac{\partial \langle\sigma_{xx}\rangle(x)}{\partial x}\cdot dx \]

where \(\tau_{yx}(y+dy)\) and \(\tau_{zx}(z+dz)\) can be expressed by the same logic.

Once we obtain the stress and pressure, we can get the force by multiplying the surface area (\(A=dx\,dy=dx\,dz= dy\,dz\)):

(254)\[F_i^x = A(-\sigma_{xx}(x)+\langle\sigma_{xx} \rangle(x+dx)-\tau_{yx}(y)+\langle\tau_{yx} \rangle(y+dy)-\tau_{zx}(z)+\langle\tau_{zx} \rangle(z+dz))\]

Plugging in (253) into (254) to replace \(\langle\sigma_{xx} \rangle(x+dx)\), we have

(255)\[F_i^x = A(\frac{\partial \sigma_{xx}(x) }{\partial x}\cdot dx+\frac{\partial \tau_{yx}(y) }{\partial y}\cdot dy+\frac{\partial \tau_{zx}(z) }{\partial z}\cdot dz) = V(\frac{\partial \sigma_{xx}(x) }{\partial x} +\frac{\partial \tau_{yx}(y) }{\partial y}+\frac{\partial \tau_{zx}(z) }{\partial z})=ma\]

since \(V=dxdydz\).

Dividing \(V\) at the both sides, we have the viscosity forces per unit volume in the x-direction as follow

(256)\[\frac{\partial \sigma_{xx}(x) }{\partial x} +\frac{\partial \tau_{yx}(y) }{\partial y}+\frac{\partial \tau_{zx}(z) }{\partial z} = \frac{m}{V}a = \rho \cdot \frac{Du}{Dt}\]

Now, the question becomes how to express normal stress \(\sigma_{xx}\) and shear stress \(\tau_{yx}\) and \(\tau_{zx}\)?

2D Stress

Let’s consider the shear stress \(\tau\) in a 2d problem

Figure 3: If rotational equilibrium, then net torque equals zero. In another word, \(\tau_{xy}\) is a reaction stress of \(\tau_{yx}\), and they are the same.

If we assume the change of angle \(\delta \alpha\) (strain) is tiny, then we can say that

(257)\[\delta \alpha \approx \tan(\delta \alpha) = \frac{\delta x}{dy} =\frac{(u_0 - (u_0+\delta u))\cdot dt}{dy}=\frac{\delta u \,dt}{dy}\]

where \(u_0\) is the original horizontal velocity of the block.

We can define the rate of the change of the angle \(\delta \alpha\) as strain rate (\(\frac{\delta \alpha}{dt}\))

(258)\[\frac{\delta \alpha}{dt}=\frac{\delta u\,dt}{dy\,dt}=\frac{du}{dy}\]

Same logic applies to the rate of change of angle \(\beta\)

(259)\[\frac{\delta \beta}{dt}=\frac{\delta v\,dt}{dx\,dt}=\frac{dv}{dx}\]

So the total rate of change of angle (\(\alpha+\beta\)) or the total strain rate due to shear stress is

(260)\[\frac{(\delta \alpha + \delta \beta)}{dt} = \frac{du}{dy}+\frac{dv}{dx}\]

Since we know that shear stress \(\tau\) equals viscosity \(\mu\) multiplies strain rate, we can have the following expression

(261)\[\tau_{xy} = \tau_{yx} = \mu\left(\frac{du}{dy}+\frac{dv}{dx}\right)\]

since we know that \(\tau_{xy}\) and \(\tau_{yx}\) is one thing, they are the reactional existence of one another.

As for normal stress \(\sigma_{xx}\) or \(\sigma_{yy}\), we have the strain from both sides, so we have doubled the strain rate from a single side.

Figure 4

(262)\[\sigma_{xx}=2\mu\dot{e}=2\mu\frac{\delta x}{dx}\div t = 2\mu\frac{(u_0 + \delta u - u_0)dt}{dx}\div dt=2\mu\frac{du}{dx}\]

(263)\[\sigma_{yy}=2\mu\frac{dv}{dy}\]

3D Stress

If we apply the same logic, we can have the normal and shear stresses in 3D problem due to viscosity alone

(264)\[\sigma_{xx}=2\mu\frac{du}{dx}\]

(265)\[\sigma_{yy}=2\mu\frac{dv}{dy}\]

(266)\[\sigma_{zz}=2\mu\frac{dw}{dz}\]

(267)\[\tau_{xy}=\tau_{yx}=\mu\left(\frac{du}{dy}+\frac{dv}{dx}\right)\]

(268)\[\tau_{yz}=\tau_{zy}=\mu\left(\frac{dw}{dy}+\frac{dv}{dz}\right)\]

(269)\[\tau_{xz}=\tau_{zx}=\mu\left(\frac{dw}{dx}+\frac{du}{dz}\right)\]

where \(w\) is the z-direction velocity.

Finally, if we plug in equations (264), (267) and (269) into the equation (256), we find that the viscous forces in the x-directions per unit volume of the element is, for constant viscosity,

(270)\[2\mu\frac{\delta^2u}{\delta x^2}+\mu\left(\frac{\delta^2 u}{\delta y^2}+\frac{\delta^2v}{\delta x \, \delta y}+\frac{\delta^2 w}{\delta x \, \delta z}+\frac{\delta^2 u}{\delta z^2}\right)\]

All of these expressions can be further simplified by substituting the continuity equation (267) differentiated wrt. \(x\).

(271)\[\frac{du}{dx}+\frac{dv}{dy}+\frac{dw}{dz}=0\]

(272)\[\frac{d^2u}{dx^2}=-\frac{d^2v}{dy\,dx}-\frac{d^2 w}{dz\,dx}\]

Plugging in equation (272) into (270), we arrive at

(273)\[\mu\left(\frac{\delta^2u}{\delta x^2}+\frac{\delta^2u}{\delta y^2}+\frac{\delta^2u}{\delta z^2}\right)\]

and same logic applies to the y and z direction

(274)\[\mu\left(\frac{\delta^2v}{\delta x^2}+\frac{\delta^2v}{\delta y^2}+\frac{\delta^2v}{\delta z^2}\right)\]

and

(275)\[\mu\left(\frac{\delta^2w}{\delta x^2}+\frac{\delta^2w}{\delta y^2}+\frac{\delta^2w}{\delta z^2}\right)\]

as the expressions for the net viscous forces per unit volume in the y and z directions, respectively.

Pressure Forces

It is not hard to find that the net pressure force on the element in the x-direction per unit area/volume of the fluid element is

(276)\[\frac{p(x)\delta y - p(x+\delta x)\delta y}{\delta x \, \delta y}=-\frac{[p(x+\delta x)-p(x)]}{\delta x}\]

which by virtue of a simple Taylor series expansion is

(277)\[-\frac{\delta p}{\delta x}\]

Same logic applies to the y and z-directions

(278)\[-\frac{\delta p}{\delta y}\]

(279)\[-\frac{\delta p}{\delta z}\]

Navier-Stoke Equation Final Form

With the easy-to-proof gravity force, we can put the three forces together to obtain the final Navier-Stoke equation we stated in the equation (223). However, let me further explain the acceleration \(a\) term on the left-hand-side of the equation. The acceleration measures the change of velocity. In a vector field for velocity, even if a particle stays at the same position, the velocity may change from time to time. For example, you can feel the water velocity changes through time when you stay in the same spot of a river. On the other hand, the velocity may also change from place to place. For example, you might find that the velocity of the water is faster at the entrence of a river than, say, at the centre of the river.

Applying the chain rule we can derive the temporal and convective acceleration in x-direction

(280)\[a_x =\frac{\text{d}u}{\text{d}t}=\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}\frac{\partial x}{\partial t}+\frac{\partial u}{\partial y}\frac{\partial y}{\partial t}+\frac{\partial u}{\partial z}\frac{\partial z}{\partial t}\]

So put all together, we should have understood the Navier-Stoke equation in 3D vector form in (248).