Navier-Stokes Derivation

The traditional derivation of the Navier-Stokes equations starts by looking at a fluid parcel and the different fluxes over the surface in the integral form. The integral form is preferred as it is more general than the differential form: For the latter one has to assume differentiability and thus it is not valid for flow discontinuities such as shocks in compressible fluids. Many of these derivations are though somewhat incomplete and often even skip the Stoke's hypothesis. Some are obscuring stuff even further by introducing concepts like the dissipation function directly without any real intermediate steps. I will try to include those steps in the following derivation.

I will make use a certain nomenclature:

• $mathcal{D}$ denotes the spatial dimensions, in the most general case of a three-dimensional flow given by ${ x, y, z }$

• every definition of a universal quantity is introduced by $:=$

The general conservation equations of a continuum in a moving reference frame

We start with the model idea of a continuum, continuous blobs of mass that completely and at least locally uniformly fill space. We therefore assume that limit value for density, stresses and forces must exist (unlike in dilute gases where this clearly breaks down). After balancing integral fluxes and applying Gauss' divergence theorem or directly balancing differential changes across a unit volume moving with the flow (Lagrangian specification) we are left with a system of following equations. The equations describes the conservation of mass eqref{1}, momentum eqref{2} and energy eqref{3} on a continuum level respectively.

$$frac{partial rho}{partial t} + sumlimits_{j in mathcal{D}} frac{partial (rho u_j )}{partial x_j }=0 tag{1}label{1}$$

$$frac{partial (rho u_i )}{partial t} + sumlimits_{j in mathcal{D}}frac{partial (rho u_i u_j )}{partial x_j} = sumlimits_{j in mathcal{D}} frac{partial sigma_{ij}}{partial x_j } + rho g_i tag{2}label{2}$$

$$frac{partial (rho e)}{partial t} + sumlimits_{j in mathcal{D}} frac{partial (rho u_j e)}{partial x_j} = - sumlimits_{j in mathcal{D}} frac{partial q_j}{partial x_j} + sumlimits_{i, j in mathcal{D}} frac{partial (sigma _{ij} u_i)}{partial x_j} + sumlimits_{j in mathcal{D}} rho u_j g_j tag{3}label{3}$$

where the total energy is given by the combination of internal $e_{in}$ and macroscopic energy $e := e_{in} + sumlimits_{j in mathcal{D}} frac{u_j u_j}{2}$. The local heat flux $q_i$ is generally assumed to be proportional to the gradient of the transported quantity, in this case the temperature, according to Fourier's law

$$q_i = - k frac{partial T}{partial x_i}. tag{4}label{4}$$

While equations eqref{1} and eqref{3} are scalar, equation eqref{2} is a vector equation with $dim(mathcal{D})$ dimensions. Furthermore the energy equation eqref{3} might be further simplified by subtracting the mechanical energy using equations eqref{2}.

Note the similar structure of the conservation equations: We are dealing with a property $phi$ which changes according to $frac{D Phi}{D t} = frac{partial Phi}{partial t} + frac{partial (phi u_i)}{partial x_i} = s$ where $s$ is a source term. For example what we perceive as a change of momentum in a particular direction $i$ due to temporal change or transport of the momentum in every potential direction $j$ with a velocity $u_j$ is equal to our sources of momentum - the change of stress pointing in direction $i$ (not only normal but also due to shearing - thus the sum over $j$) and the specific body force in the direction $i$.

Material law for a single dimension

So far we have made no assumptions on the precise material law yet. In the end a continuum is nothing but a dense more or less uniform chunk of mass that has properties of a mass oscillator (Kelvin-Voigt model). We pretend that we are moving with the center of mass and the volume of interest expands and contracts exerting forces similar to a combination of a linear spring and a linear dampener. Neglecting external forces such as gravity we end up with the equation of motion for a linear one-dimensional continuum element

$$F = c , x + d , dot{x} - p A.$$

Dividing by area $A$ we yield the stress $sigma$

$$sigma := frac{F}{A} = frac{c , x + d , dot{x}}{A} - p,$$

which can be further rewritten introducing the Young's modulus $E := c , L / A$, the viscosity $mu := d , L / A$, the displacement $epsilon := x / L$ and its time derivative $dot{epsilon} := dot{x} / L$ to

$$sigma = E , epsilon + mu , dot{epsilon} - p.$$

This model still includes two parts of deformation: A reversible elastic - characteristic for a lot of solid materials - and a time-dependent irreversible deformation - usually found in fluids and is thus referred to as a viscoelastic material.

Material law for multiple dimensions

Extending this model to multiple dimensions requires characteristic invariant deformation measures that do not change with rigid body motions (translation and rotation). The symmetric part (the anti-symmetric part results in rigid body rotations) of the gradients (this excludes rigid body translations) of the corresponding quantities $vec x$ and $vec u$ is thus a plausible choice

$$underline{epsilon} := frac{1}{2} left( vec nabla otimes vec x + (vec nabla otimes vec x)^T right) phantom{abcdefg} underline{S} := frac{1}{2} left( vec nabla otimes vec u + (vec nabla otimes vec u)^T right).$$

(I am not a big fan of vector annotation with the outer product so we will drop them in a second and replace them by index annotation!)

Now we need a multi-linear relationship between two different tensors of rank two, $sigma_{ij}$ and $epsilon_{kl}$ and $S_{kl}$. This is given by a tensor of rank four as follows

$$sigma_{ij} = underbrace{-p delta_{ij}}_{sigma_{ij}^{(0)}} + underbrace{sumlimits_{k,l} C_{ijkl} epsilon_{kl} + sumlimits_{k,l} D_{ijkl} S_{kl}}_{tau_{ij}}$$

where $sigma_{ij}^{(0)}$ is the stress distribution in the resting state of the continuum of interest, $tau_{ij}$ the stresses arising from deformation, while $C_{ijkl}$ and $D_{ijkl}$ are fourth rank tensor, with $3^4=81$ coefficients each in their most general form, that establish the linear relationship between deformations and stresses.

In the classical theory of linear elasticity it is assumed that the viscous part may be neglected, resulting in the total stress tensor

$$sigma_{ij} = -p delta_{ij} + underbrace{sumlimits_{k,l} C_{ijkl} epsilon_{kl}}_{tau_{ij}}.$$

This corresponds to a linear spring that stores all the supplied energy and releases it at a later point and is the common approach in strength of materials, another continuum science.

For an isotropic material, where there is no preferred direction, $C_{ijkl}$ and $D_{ijkl}$ clearly have to be a isotropic tensors of rank four. In order to derive the form of an isotropic tensor of rank four $underline{T}$ we introduce a scalar $s$ that is obtained from vectors $vec a$, $vec b$, $vec c$, $vec d$ and a tensor $underline{T}$ as begin{equation} s = sumlimits_{i,j,k,l} T_{ijkl} a_i b_j c_k d_l. end{equation} Without making any assumptions on the fourth-order tensor $underline{T}$ the scalar $s$ would depend linearly on the magnitude of every single vector $in { vec a, vec b, vec c, vec d }$ and their relative orientation in space. If we now assume an isotropic tensor the precise direction of the four vectors should not affect the scalar $s$ but instead only the orientation of the vectors to another given by the dot product begin{equation} begin{aligned} s = sumlimits_{i,j,k,l} T_{ijkl} a_i b_j c_k d_l stackrel{!}{=} alpha left(vec a cdot vec b right)left(vec c cdot vec d right) + beta left(vec a cdot vec c right)left(vec b cdot vec d right) + gamma left(vec a cdot vec d right)left(vec c cdot vec b right) = = sumlimits_{i,j} left( alpha , a_i b_i c_j d_j + beta , a_i c_i b_j d_j + gamma , a_i d_i c_j b_j right) = sumlimits_{i,j,k,l} underbrace{left( alpha , delta_{ij} delta_{kl} + beta , delta_{ik} delta_{jl} + gamma , delta_{il} delta_{jk} right)}_{T_{ijkl}} a_i b_j c_k d_l. end{aligned} end{equation}

The resulting isotropic tensor of rank four

$$ T_{ijkl} = alpha , delta_{ij} delta_{kl} + beta , delta_{ik} delta_{jl} + gamma , delta_{il} delta_{jk} $$

can be further simplified by enforcing Cauchy's second law of motion, the conservation of angular momentum, which results in the symmetry of the Cauchy stress tensor

$$tau_{ij} = tau_{ji}$$

and due to

$$tau_{ij} = sumlimits_{k,l} C_{ijkl} epsilon_{kl} = sumlimits_{k,l} left( alpha delta_{ij} delta_{kl} + beta delta_{ik} delta_{jl} + gamma delta_{il} delta_{jk} right) epsilon_{kl}$$

$$tau_{ji} = sumlimits_{k,l} C_{jikl} epsilon_{kl} = sumlimits_{k,l} left( alpha delta_{ji} delta_{kl} + beta delta_{jk} delta_{il} + gamma delta_{jl} delta_{ik} right) epsilon_{kl}$$

(exchanging the two indices $i$ and $j$ due to the symmetry) it follows by comparison that only two coefficients $alpha$ and $beta = gamma$ are actually independent.

This way the isotropic elasticity tensor degenerates to

$$C_{ijkl} = lambda delta_{ijkl} + mu (delta_{ik} delta_{jl} + delta_{il} delta_{jk})$$

where the two independent coefficients are termed the first $lambda$ and second $mu$ Lame's parameter.

Analogously one can rewrite the constitutive equation for a viscous Newtonian fluid, neglecting the elastic contribution, with the viscosity tensor $D_{ijkl}$ and the strain-rate tensor $S_{ij}$, the symmetric part of velocity gradient, given by

$$S_{ij} := frac{1}{2} left( vec nabla otimes vec u + (vec nabla otimes vec u)^T right) = frac{1}{2} left( frac{partial u_i}{partial x_j} + frac{partial u_j}{partial x_i} right) $$

to

$$sigma_{ij} = sigma_{ij}^{(0)} + sumlimits_{k,l} D_{ijkl} S_{kl} = - p delta_{ij} + tau_{ij}.$$

This basically corresponds to a viscous linear dampener that dissipates energy! (For non-Newtonian fluids, where the flow properties are different from the aforementioned Newtonian fluid, e.g. the viscosity is a function of the shear-rate, other approaches have to be used. Even though almost all fluids show more or less non-Newtonian behaviour, most, like water and air, can be approximated as Newtonian.)

The stress distribution in the resting state corresponds to the hydrostatic pressure $p_0$ in a resting and to the thermodynamic pressure $p$ in a moving fluid. Generally it is assumed that rotational viscosity, the rate at which angular momentum differences are equilibrated, may be neglected resulting in symmetric viscous stress tensor $tau_{ij} = tau_{ji}$. The two resulting material-dependent coefficients are referred to as shear viscosity $mu$ and dilatational viscosity $lambda$

$$tau_{ij}=mu left(frac{partial u_i}{partial x_j}+frac{partial u_j}{partial x_i} right) + lambda sumlimits_{k in mathcal{D}} frac{partial u_k}{partial x_k} delta_{ij}.$$

Introducing the mean mechanical pressure $overline{p}$, analogously to the hydrostatic stress in linear elasticity theory $pi$,

$$overline{p} := - frac{1}{dim(mathcal{D})} sumlimits_{j in mathcal{D}} sigma_{jj} = - frac{1}{3} left( sigma_{11} + sigma_{22} + sigma_{33} right) = p - underbrace{left( lambda + frac{2}{3} mu right)}_kappa sumlimits_{k in mathcal{D}} frac{partial u_k}{partial x_k}$$

leads to a surprising result: Unless either the divergence of the velocity (incompressible flow) or the term $kappa := lambda + frac{2}{3} mu$, often referred to as bulk viscosity, are zero, the mechanical pressure is not equivalent to the thermodynamic pressure: The bulk viscosity introduces additional dissipation during a change of volume, where shear forces are not present.

Stokes simply assumed a vanishing bulk viscosity (Stokes' hypothesis) which left him with a stress tensor according to

$$sigma_{ij} = - p delta_{ij} + 2 mu S_{ij} - frac{2}{3} mu sumlimits_{k in mathcal{D}} S_{kk} delta_{ij}. tag{5}label{5}$$

Buresti suggests this should rather be seen as $overline{p} approx p$. And even though this assumption is widely used throughout fluid dynamics it seems as if it would only hold for mono-atomic gases, while in the case of poly-atomic gases molecular interactions are probably responsible for a thermodynamic pressure that deviates from the mechanical pressure.

Three possible basic deformations of a continuum element can be identified: angular and linear deformation as well as volumetric dilatation. While the latter is caused by the first two terms in the stress tensor the last term results in angular deformation for $i neq j$ and linear deformation for $i = j$.

The Navier-Stokes-Fourier-equations

Equations eqref{1} - eqref{5} leaves us with the the full Navier-Stokes-Fourier equations. As they can be pretty unhandy to write - by inserting all the involved laws - most literature introduces certain concepts like the aforementioned dissipation function.

Starting from there one may again simplify the equation system by assuming certain material parameters (density $rho$, viscosity $nu$ or thermal conductivity $k$) as constant or neglect terms depending on their relative order of magnitude.