Published on

Navier-Stokes Equations

Authors
  • avatar
    Name
    Jyotir Sai
    Twitter
    Engineering Student

The Navier-Stokes equations are fundamental equations describing fluid flow. I'll be going over the assumptions and math used to derive these equations. You should be familiar with calculus, differential equations, and some basic physics as this makes up the bulk of the derivation.

Eulerian Description of Fluid Flow

The Navier-Stokes equations form the foundation of modern fluid mechanics. They are essentially conservation of momentum equations derived from Newton's second law on a fluid. However, there are several concepts that need to be established on the journey towards Navier-Stokes. Let's first talk about the Eulerian description of flow. Rather than following a particular fluid particle throughout its motion (Lagrange description), we will define an imaginary control volume and evaluate the particles that pass in and out of this control volume. Therefore, instead of defining the velocity of a specific fluid particle, we define the velocity field of the fluid that passes through the control volume.

V=u(x,y,z,t)i+v(x,y,z,t)j+w(x,y,z,t)k\vec{V}=u(x,y,z,t)\vec{i}+v(x,y,z,t)\vec{j}+w(x,y,z,t)\vec{k}

Fluid elements passing through the control volume at a particular location and time will have the velocity defined above. The derivative of the velocity yields acceleration, but this is only for Lagrangian particles. When using the Eulerian description, we must use the material derivative to calculate acceleration from the velocity field. The material derivative is defined as:

DFDt=Ft+(V)F\frac{DF}{Dt}=\frac{\partial F}{\partial t}+\left (\vec{V}\cdot \vec{\nabla } \right )F

The second term is the divergence of the velocity field which is defined as:

divV=V=(x,y,z)(Vx,Vy,Vz)=Vxx+Vyy+Vzzdiv \bold V = \nabla\cdot \bold V=\left ( \frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z} \right )\cdot \left ( V_{x}, V_{y}, V_{z} \right ) = \frac{\partial V_{x}}{\partial x} + \frac{\partial V_{y}}{\partial y}+\frac{\partial V_{z}}{\partial z}

The material acceleration is therefore:

a(x,y,z,t)=DVDt=Vt+(V)V\vec{a}(x,y,z,t)=\frac{D\vec{V}}{Dt}=\frac{\partial \vec{V}}{\partial t}+\left ( \vec{V}\cdot \vec{\nabla} \right )\vec{V}
a(x,y,z,t)=Vt+uVx+vVy+wVz\vec{a}(x,y,z,t)=\frac{\partial \vec{V}}{\partial t}+u\frac{\partial \vec{V}}{\partial x}+v\frac{\partial \vec{V}}{\partial y}+w\frac{\partial \vec{V}}{\partial z}

The material derivative offers a link between a Langragian and Eulerian descriptions of fluid flow.

Reynold's Transport Theorem (RTT)

In fluid mechanics, it is more common to work with control volumes rather than a specific system of particles. In a system approach, we assume a closed system i.e. mass cannot cross our system boundary. The control volume approach uses an open system where mass can go in and out of the boundaries. Most conservation laws are derived using a system rather than a control volume.

The above figure serves to illustrate the difference between a system and a control volume. The system encompasses all of the mass and will expand with it. The control volume does not expand with the mass but it does record the mass entering and leaving its boundaries (along with the mass inside). Reynold's Transport Theorem relates a system to a control volume. This allows us to derive conservation laws in terms of a control volume. We can see from the above figure that the change in the system is equal to the change in the "stuff" inside the control volume plus the net change going in and out of the control volume.

dBsysdt=dBCVdt+B˙outB˙in\frac{dB_{sys}}{dt}=\frac{dB_{CV}}{dt}+\dot{B}_{out}-\dot{B}_{in}

BB simply represents an extensive property. An extensive property is a property that is proportional to the size or quantity of matter in a system such as mass, volume, or entropy. An intensive property does not depend on the size or quantity of matter and is related to an extensive property via b=Bmb=\frac{B}{m}. The B˙\dot{B} term can be re-written as bm˙b\dot{m} where m˙\dot{m} is the mass flow rate. The mass flow rate can be further broken down into m˙=ρVA\dot{m} = \rho \vec{V} A where ρ\rho is density, V\vec{V} is velocity, and AA is area.

The net flux (B˙outB˙in\dot{B}_{out}-\dot{B}_{in}) can be written as:

B˙outB˙in=CSbρVndA\dot{B}_{out}-\dot{B}_{in}=\int_{CS}^{}b\rho\vec{V}\cdot \vec{n} dA

n\vec{n} is the normal vector and its dot product with V\vec{V} ensures that we get the normal velocity. We integrate over the entire differential surface area dAdA (inlet or outlet area). We write it this way so we can account for multiple inlets and outlets of differing sizes. For the control volume, BCVB_{CV} can be written as cvbρdV\int_{cv}^{}b\rho d\boldsymbol{\mathbf{V}} where dVd\boldsymbol{\mathbf{V}} is the differential volume. The original expression relating the system to the control volume is now:

dBsysdt=ddtCVbρdV+CSbρVndA\frac{dB_{sys}}{dt}=\frac{d}{dt}\int_{CV}^{}b\rho d\boldsymbol{\mathbf{V}}+\int_{CS}^{}b\rho\vec{V}\cdot \vec{n} dA

This is Reynold's Transport Theorem for a fixed control volume. For a moving control volume we simply replace the velocity in the above expression with the relative velocity Vrel\vec{V_{rel}} equal to:

Vrel=VVCS\vec{V_{rel}}=\vec{V}-\vec{V_{CS}}

VCS\vec{V_{CS}} is the velocity of the control surface.

Conservation of Mass

In Reynold's Transport Theorem, we looked at a generic extensive property called BB and its intensive property bb. In the case of mass conservation, the extensive property is equal to mass (B=mB=m) and therefore the intensive property is equal to b=mm=1b=\frac{m}{m}=1. Let's write RTT in terms of mass:

dmsysdt=ddtCVρdV+CSρVndA\frac{dm_{sys}}{dt}=\frac{d}{dt}\int_{CV}^{}\rho d\boldsymbol{\mathbf{V}}+\int_{CS}^{}\rho\vec{V}\cdot \vec{n} dA

The left hand side of the above equation is the rate of change of mass in the system. However, we know that from a system perspective the mass does not change. Unlike in a control volume, the system will continually transform along with the mass. Therefore, the change in mass of the system is equal to 0.

0=ddtCVρdV+CSρVndA0=\frac{d}{dt}\int_{CV}^{}\rho d\boldsymbol{\mathbf{V}}+\int_{CS}^{}\rho\vec{V}\cdot \vec{n} dA

We can simplify the above expression by writing both terms as a volume integral.

0=CV[ρt+(ρV)]dVcv0=\int_{CV}^{}\left [ \frac{\partial \rho}{\partial t}+\nabla\cdot \left ( \rho \vec{V} \right ) \right ]dV_{cv}

In order for the above expression to equal 0, the expression in the integrand must also equal 0:

ρt+(ρV)=0\frac{\partial \rho}{\partial t}+\nabla\cdot \left ( \rho \vec{V} \right ) = 0

This is the general equation for conservation of mass, also known as the continuity equation. It can be further expanded into the following using cartesian coordinates:

ρt+Vxx+Vyy+Vzz=0\frac{\partial \rho}{\partial t}+\frac{\partial V_{x}}{\partial x}+\frac{\partial V_{y}}{\partial y}+\frac{\partial V_{z}}{\partial z}= 0

Conservation of Momentum

For this case, our extensive property is momentum (mVm\vec{V}), so from RTT we have:

d(mV)sysdt=ddtCVVρdV+CSρVVndA\frac{d(m\vec{V})_{sys}}{dt}=\frac{d}{dt}\int_{CV}^{}\vec{V}\rho d\boldsymbol{\mathbf{V}}+\int_{CS}^{}\rho\vec{V}\vec{V}\cdot \vec{n} dA

First, let's take a look at the left hand side. Our mass in the system is not changing so we can take out that term from the derivative. That leaves us with the derivative of the velocity which is just acceleration. Therefore, we have mass times acceleration which is simply equal to the net force on the system.

F=ddtCVVρdV+CSρVVndA\sum F=\frac{d}{dt}\int_{CV}^{}\vec{V}\rho d\boldsymbol{\mathbf{V}}+\int_{CS}^{}\rho\vec{V}\vec{V}\cdot \vec{n} dA

The forces acting on the fluid includes surface (ex. pressure) and body (ex. gravity) forces. We can simplify the left hand side by assuming a Newtonian fluid with the following stress tensor:

T=(p+23μVkxk)δij+μ(Vixj+Vjxi)T=-\left ( p+\frac{2}{3}\mu \frac{\partial V_{k}}{\partial x_{k}} \right )\delta_{ij}+\mu\left ( \frac{\partial V_{i}}{\partial x_{j}}+\frac{\partial V_{j}}{\partial x_{i}} \right )

where pp is pressure, μ\mu is the dynamic viscosity of the fluid, and δij\delta_{ij} is the Kronecker delta. The Kronecker delta is equal to 1 when i=ji=j and equal to 0 when iji \neq j. If we assume that the fluid is incompressible (which is a reasonable assumption), then the stress tensor can be evaluated to the following.

T=[P+2μuxμ(uy+vx)μ(uz+wx)μ(vx+uy)P+2μvyμ(vz+wy)μ(wx+uz)μ(wy+vz)P+2μwz]T=\begin{bmatrix} -P+2\mu \frac{\partial u}{\partial x} & \mu\left ( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right ) & \mu\left ( \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x} \right ) \\ \mu\left ( \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y} \right ) & -P+2\mu \frac{\partial v}{\partial y} & \mu\left ( \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y} \right ) \\ \mu\left ( \frac{\partial w}{\partial x}+\frac{\partial u}{\partial z} \right ) & \mu\left ( \frac{\partial w}{\partial y}+\frac{\partial v}{\partial z} \right ) & -P+2\mu \frac{\partial w}{\partial z} \end{bmatrix}

The left hand side is now equal to:

F=Fsurfaces+Fbody=CSTndA+CVρbdV\sum F = \sum F_{surfaces} + \sum F_{body} = \int_{CS}^{}T \cdot \vec{n}dA + \int_{CV}^{}\rho \vec{b}dV

The surface forces are inside a surface integral. b\vec{b} is the body forces per unit mass. The momentum expression is now:

CSTndA+CVρbdV=ddtCVVρdV+CSρVVndA\int_{CS}^{}T \cdot \vec{n}dA + \int_{CV}^{}\rho \vec{b}dV = \frac{d}{dt}\int_{CV}^{}\vec{V}\rho d\boldsymbol{\mathbf{V}}+\int_{CS}^{}\rho\vec{V}\vec{V}\cdot \vec{n} dA

We have terms in surface and volume integrals so we will write everything to be inside a volume integral using the divergence theorem.

CSρVVndA=CV(ρVV)dV\int_{CS}^{}\rho\vec{V}\vec{V}\cdot \vec{n} dA = \int_{CV}^{} \vec{\nabla} \cdot \left ( \rho \vec{V}\vec{V} \right )dV
CSTndA=CVTdV\int_{CS}^{}T \cdot \vec{n}dA = \int_{CV}^{} \vec{ \nabla} \cdot T dV

Now we can rearrange everything inside a volume integral.

CV[t(ρV)+(ρVV)ρbT]dV=0\int_{CV}^{}\left [ \frac{\partial}{\partial t} \left ( \rho \vec{V} \right) + \vec{\nabla} \cdot \left ( \rho \vec{V} \vec{V} \right) - \rho \vec{b} - \vec{ \nabla} \cdot T \right ]dV = 0

The expression inside the integrand must also equal zero for the above to hold true. Therefore, we can simplify to the following:

t(ρV)+(ρVV)=ρb+T\frac{\partial}{\partial t} \left ( \rho \vec{V} \right) + \vec{\nabla} \cdot \left ( \rho \vec{V} \vec{V} \right) = \rho \vec{b} + \vec{ \nabla} \cdot T

The expressions on the left can be expanded into the following.

t(ρV)=ρvt+Vρt\frac{\partial}{\partial t} \left ( \rho \vec{V} \right) = \rho \frac{\partial \vec{v}}{\partial t}+\vec{V}\frac{\partial \rho}{\partial t}
(ρVV)=V(ρV)+ρ(V)V\vec{\nabla} \cdot \left ( \rho \vec{V} \vec{V} \right) = \vec{V}\vec{\nabla}\cdot\left ( \rho \vec{V} \right)+\rho\left (\vec{V}\cdot\vec{\nabla} \right)\vec{V}

Plugging the above back into our expression we have:

ρVt+Vρt+V(ρV)+ρ(V)V=ρb+T\rho \frac{\partial \vec{V}}{\partial t}+\vec{V}\frac{\partial \rho}{\partial t} + \vec{V}\vec{\nabla}\cdot\left ( \rho \vec{V} \right)+\rho\left (\vec{V}\cdot\vec{\nabla} \right)\vec{V} = \rho \vec{b} + \vec{ \nabla} \cdot T

On the left hand side, we have the expression Vρt+V(ρV)\vec{V}\frac{\partial \rho}{\partial t} + \vec{V}\vec{\nabla}\cdot\left ( \rho \vec{V} \right) which can be written as:

V(ρt+(ρV))\vec{V}\left (\frac{\partial \rho}{\partial t}+ \vec{\nabla}\cdot\left ( \rho \vec{V} \right)\right)

The expression in brackets is simply our conservation of mass expression from before which is equal to 0! So now we have:

ρVt+ρ(V)V=ρb+T\rho \frac{\partial \vec{V}}{\partial t}+\rho\left (\vec{V}\cdot\vec{\nabla} \right)\vec{V}=\rho \vec{b} + \vec{ \nabla} \cdot T

The left hand side can be reduced even further. We can write ρVt+ρ(V)V\rho \frac{\partial \vec{V}}{\partial t}+\rho\left (\vec{V}\cdot\vec{\nabla} \right)\vec{V} as ρ[Vt+(V)V]\rho \left [\frac{\partial \vec{V}}{\partial t}+\left (\vec{V}\cdot\vec{\nabla} \right)\vec{V}\right] where the expression in brackets is the definition of the material derivative. Therefore, we have:

ρDvDt=T+ρb\rho \frac{D\vec{v}}{Dt}=\vec{\nabla}\cdot T+\rho \vec{b}

If the only body force that acts on the CV is gravity then we can write b=g\vec{b}=\vec{g}. Let's expand the above expression into cartesian coordinates.

xcomponent:ρDuDt=T11x+T21y+T31z+ρgxx-component: \rho \frac{Du}{Dt}=\frac{\partial T_{11}}{\partial x}+\frac{\partial T_{21}}{\partial y}+\frac{\partial T_{31}}{\partial z}+\rho g_{x}
ycomponent:ρDvDt=T12x+T22y+T32z+ρgyy-component: \rho \frac{Dv}{Dt}=\frac{\partial T_{12}}{\partial x}+\frac{\partial T_{22}}{\partial y}+\frac{\partial T_{32}}{\partial z}+\rho g_{y}
zcomponent:ρDwDt=T13x+T23y+T33z+ρgzz-component: \rho \frac{Dw}{Dt}=\frac{\partial T_{13}}{\partial x}+\frac{\partial T_{23}}{\partial y}+\frac{\partial T_{33}}{\partial z}+\rho g_{z}

Let's go over evaluating the derivatives of the stress tensor terms. We'll start with the x component.

T11x=Px+2μ2ux2\frac{\partial T_{11}}{\partial x} = \frac{-\partial P}{\partial x}+2\mu\frac{\partial^{2}u}{\partial x^{2}}
T21y=μy(vx+uy)\frac{\partial T_{21}}{\partial y} = \mu \frac{\partial}{\partial y} \left ( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right)
T31z=μz(wx+uz)\frac{\partial T_{31}}{\partial z} = \mu \frac{\partial}{\partial z} \left ( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right)

The x-component can now be written as:

ρDuDt=Px+2μ2ux2+μy(vx+uy)+μz(wx+uz)+ρgx\rho \frac{Du}{Dt}=\frac{-\partial P}{\partial x}+2\mu\frac{\partial^{2}u}{\partial x^{2}}+\mu \frac{\partial}{\partial y} \left ( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right)+\mu \frac{\partial}{\partial z} \left ( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z} \right)+\rho g_{x}

Which can be rearranged into:

ρDuDt=Px+μ[x(ux+vy+wz)+2ux2+2uy2+2uz2]+ρgx\rho \frac{Du}{Dt}=\frac{-\partial P}{\partial x}+\mu \left [\frac{\partial}{\partial x} \left (\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \right) +\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}} \right]+\rho g_{x}

The term in parentheses is simply 0 due to mass conservation and the squared partial derivatives are the Laplacian (2\nabla^{2}) of the uu.

ρDuDt=Px+μ2u+ρgx\rho \frac{Du}{Dt}=\frac{-\partial P}{\partial x}+\mu \nabla^{2} u+\rho g_{x}

The y and z terms can similarly be reduced to:

ρDvDt=Py+μ2v+ρgy\rho \frac{Dv}{Dt}=\frac{-\partial P}{\partial y}+\mu \nabla^{2} v+\rho g_{y}
ρDwDt=Pz+μ2w+ρgz\rho \frac{Dw}{Dt}=\frac{-\partial P}{\partial z}+\mu \nabla^{2} w+\rho g_{z}

The three expressions can be written as one vector equation:

ρDVDt=P+ρg+μ2V\rho \frac{D\vec{V}}{Dt}=-\vec{\nabla}P+\rho\vec{g}+\mu\nabla^{2}\vec{V}

This final expression is known as the Navier-Stokes equation. The above equation needs to be used in conjunction with the continuity equation when solving problems. This equation is a non-linear, 2nd-order partial differential equation. The solutions to this equation are simplified cases with specific boundary conditions. There is no general solution to the above equation and it is in fact one of the millenium prize problems.