Navier-Stokes Equation

In the previous section ({doc}./ContinuumMechanics), we looked a very general way use conservation principles to construct equations describing how sources, fluxes and material transport of varous quantities balance one another through time. This approach can also be used to generate conservation equations for additional physical variables (for example, angular momentum, electric current).

One thing is still missing, though, and that is how we express the flux term for a quantity in terms of gradients in the quantity itself. We did this without any comment in the conservation of heat energy because we used Fourier’s law \(\mathbf{F} = -k \nabla T\) to express the heat flux as a linear function of the temperature gradient. This is an empirical relationship that contains a physical constant that we are not able to derive from first principles. This missing relationship is known as a constitutive law.

Viscosity

For a fluid, the constitutive relationship relating momentum flux to momentum gradients introduces a new physical parameter, viscosity which relates stress to rates of change in velocity gradients and a second viscosity relating to rates of volume change.

If we make the assumption that the fluid is incompressible (\(\nabla \cdot \mathbf{v}\)) which we already mentioned in deriving the temperature equation, we only have to worry about the viscosity that relates to velocity gradients

\[ \sigma_{ij} = \eta \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right) - p\delta_{ij} \]

(An expanded derivation for this is in Landau & Lifschitz[REF] and they also talk extensively about rate of volume change effects).

If the viscosity is constant, then we can substitute the constitutive law into the stress-divergence term of the momentum conservation equation. In index notation once again,

\[ \begin{split} \nabla \cdot \boldsymbol{\sigma} & = \frac{\partial}{\partial x_j} \eta \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) - \frac{\partial p}{\partial x_j} \delta_{ij} \\ & = \eta \frac{\partial^2 v_i}{\partial x_j \partial x_j} + \eta \frac{\partial^2 v_j}{\partial x_i \partial x_j} - \frac{\partial p}{\partial x_i}\\ & = \eta \nabla^2 \mathbf{v} + \eta \textcolor[rgb]{0.0,0.7,0.0}{ \nabla (\nabla \cdot \mathbf{v})} - \nabla p \end{split} \]

the second (green) term in this final form must vanish because of the incompressibility assumption, so the momentum conservation equation becomes

This is the Navier-Stokes equation. Once again some new notation has shown up uninvited. The Laplacian operator \(\nabla^2\) is defined (in a scalar context) as

and in a vector context as

\[ \begin{split} \nabla^2 \mathbf{u} & = \nabla \nabla \cdot \mathbf{u} - \nabla \times (\nabla \times \mathbf{u}) \\ & = \mathbf{i} \nabla \cdot \nabla u_x + \mathbf{j} \nabla \cdot \nabla u_y + \mathbf{k} \nabla \cdot \nabla u_z \mathrm{\hspace{1cm} (Cartesian)} \end{split} \]

In Cartesian coordinates, the Laplacian operator has a simple form, and the vector Laplacian is simply the scalar operator applied in each direction. In other coordinate systems this operator becomes substantially more complicated but this notation is handy because the equations always look the same no matter the coordinate system.

Boussinesq Approximation, Equation of State, Density Variations

The equation which relates pressure, temperature and density is known as the equation of state. For the equations derived so far, we have specified an incompressible fluid and that means the density can’t change at all.

Can you see the problem ? We have also included a source term for momentum which relies on gravity acting on density variations. This conflict is typical in fluid mechanics: simplifying assumptions if taken to their logical limit imply only trivial solutions (no motions are possible, that sort of thing)

We have to do something to get around this problem and the simplest thing to do is assume density changes are typically small relative to the magnitude of the density. Any terms in the equations which are multiplied by density can therefore assume that density is a large, constant value. Terms which contain gradients of density or density variations (where the constant terms cancel out) need to consider the equation of state.

We call this the Boussinesq approximation and it is only good for nearly- incompressible fluids. In the Navier-Stokes equation, the hydrostatic pressure does not influence the velocity field at all. Only differences in density drive fluid flow, and so the only term where we need to worry about density not being constant is in the gravitational body forces.

In the case of density variations due to temperature, the equation of state is this:

\[ \rho = \rho_0 \left(1 - \alpha ( T-T_0 )\right) \tag{10.1}\] where \(\rho_0\) is the density at a reference temperature \(T_0\). \(\alpha\) is the coefficient of thermal expansion. \(\alpha\) is generally much smaller than one which justifies making the approximation in the first place.

The energy and momentum conservation equations are now coupled through the buoyancy forces,

\[ g\rho\hat{\mathbf{z}} = g \rho_0 \left(1 - \alpha(T-T_0)\right) \]

Density variations due to pressure produce a perfectly vertical, isotropic forcing term on the momentum conservation equation. In the steady state case, this is balanced by the hydrostatic pressure gradient. (The isotropic term does not contribute at all the the deviatoric part of the stress equation and thus cannot induce steady flow). This is the reason we can ignore the vertical density gradient due to the weight of the fluid.

Another density variation is the one that results from variation in chemical composition from one fluid element to another. An example would be where two immiscible fluids live in the same region. Now density variations might be large – the fluid domains must be considered separately if we still want to assume small density changes due to temperature.

Fluid Motion: Eulerian v. Lagrangian

We have seen the \(\mathbf{v} \cdot \nabla\) operator a number of times now. The presence of this term causes major complications when we want to find solutions to the equations of continuum mechanics since it introduces a strong non-linearity. It is this term which is behind turbulence in high speed flows.

This term describes the advection of a quantity which accounts for the passive transport of something (temperature, momentum), by the motion of the fluid. Advection also presents some serious headaches in numerical methods and has spawned a vast literature devoted to efficient and accurate solution methods.

One way to side-step
the complexity of advection is to consider an elemental volume of space which moves with the fluid. The surface flux term from equation Equation 9.1 vanishes at once if there is no motion of material in or out of the boundary. In reality, we are just sweeping the problem under the rug – it doesn’t actually go away altogether – but it can be useful to do this because every once in a while that term will cancel out before we need to actually do anything with it ! And, sometimes, when it comes to numerical methods, different views of the equations can lead to better approximations.

Mathematically, we introduce a new notation:

\[ \frac{D \phi}{D t} = \frac{d}{dt} \phi\left[ x_1(t),x_2(t),x_3(t),t \right] \]

where the change in the reference position \((x_1(t),x_2(t),x_3(t))\) is governed by the local flow velocity:

\[ \frac{d x_1}{d t} = v_1 \;\;\;\; \frac{d x_2}{d t} = v_2 \;\;\;\; \frac{d x_3}{d t} = v_3 \]

and this keeps the reference point moving with the fluid. Differentiating gives

\[ \frac{D \phi}{D t} = \frac{\partial \phi}{\partial t} + \frac{\partial \phi}{\partial x_1}\frac{d x_1}{d t} + \frac{\partial \phi}{\partial x_2}\frac{d x_2}{d t} + \frac{\partial \phi}{\partial x_3}\frac{d x_3}{d t} \]

which leads to

\[ \frac{D \phi}{D t} = \frac{\partial \phi}{\partial t} + v_1 \frac{\partial \phi}{\partial x_1} + v_2 \frac{\partial \phi}{\partial x_2} + v_3 \frac{\partial \phi}{\partial x_3} \]

which is equivalent to

\[ \frac{D \phi}{D t} = \frac{\partial \phi}{\partial t} + (\mathbf{v} \cdot \nabla) \phi \]

If we think of \(\phi\) as the concentration of a dye in the fluid, then the above is a conservation equation assuming the dye does not diffuse and has no sources or sinks.

Viewed from a reference frame locked to a particular fluid element, the energy conservation equation becomes

\[ \frac{D T}{Dt} = \kappa \nabla^2 T + \frac{H}{C_p} \]

and the momentum conservation equation now becomes

\[ \frac{D \mathbf{v} }{D t} = \eta \nabla^2 \mathbf{v} - \nabla P - g\rho\hat{\mathbf{z}} \]

This is a considerably more compact way of writing the equations, but it is now necessary to use a coordinate system which is locked into the fluid and rapidly deforms as the fluid flows. Before long, the coordinate system is unimaginably complex – the advection problem returns in another guise. This formulation is known as the Lagrangian formulation and contrasts with the Eulerian viewpoint which is fixed in space.

From the numerical point of view, however, this approach can have some distinct advantages. The computer can often track the distorted coordinate system far better than it can handle successive applications of the \(\mathbf{v} \cdot \nabla\) operator at a fixed point in space. We will return to this point later.

Reading Material

G. F. Davies. Dynamic Earth. Cambridge University Press, New York, 1999.

J. Grotzinger, T. H. Jordan, F. Press, and R. Siever. Understanding Earth. W. H. Freeman & Co, 5 edition, 2006. URL http://bcs.whfreeman.com/understandingearth5e/.

L. D. Landau and E. M. Lifshitz. Fluid Mechanics, volume 6 of Course of Theoretical Physics. Pergamon Press, 1959.

G. Schubert, D. L. Turcotte, and P. Olson. Mantle Convection in Earth and Planets. Cambridge University Press, UK, 2001.

D.L. Turcotte and G. Schubert. Geodynamics. John Wiley and Sons, New York, 1982.

Exercises

This exercises section was newly added and adapted from MTH3360 Fluid Mechanics course materials. The exercises have been modified to include geodynamic applications and use the book’s notation conventions, but may need further refinement for pedagogical flow and difficulty progression.

Changes made: - Adapted 8 exercises from ICFM Chapter 2 (Equations of Motion) - Added 2 new exercises on Boussinesq approximation and Reynolds number for geodynamic scenarios - Converted LaTeX notation to Quarto format with global math macros - Added detailed solutions to 5 exercises in collapsible callout blocks - Modified examples to use geodynamic contexts where appropriate

Note

Exercises marked with ⭐ have worked solutions provided.

Exercise 1: Material Derivative ⭐

Under what condition is the advective rate-of-change equal to the total rate-of-change? In other words, when does: \[ \frac{D\mathbf{f}}{Dt} = (\mathbf{u}\cdot \nabla) \mathbf{f} \]

Explain the physical meaning of this condition.

Exercise 2: Isothermal Atmosphere ⭐

Show that, in hydrostatic equilibrium, the pressure and density in an isothermal (constant temperature) atmosphere vary with height according to: \[ p(z) = p(0)\exp(-z/H_s) \quad \text{and} \quad \rho(z) = \rho(0)\exp(-z/H_s) \]

where \(H_s = RT/g\), \(R = 287\) J K\(^{-1}\) kg\(^{-1}\) is the specific gas constant, and \(z\) points vertically upwards.

Show that for realistic values of \(T\) in the troposphere, the e-folding height scale \(H_s\) is on the order of 8 km.

Exercise 3: Vector Differential Operators ⭐

The vector differential operator del (nabla) is defined as: \[ \nabla= \left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right) \]

Express in full Cartesian form the quantities: - \(\nabla\cdot\mathbf{u}\) (divergence) - \(\nabla\times\mathbf{u}\) (curl) - \(\mathbf{u}\cdot \nabla\) (directional derivative) - \(\nabla\cdot \nabla\equiv \nabla^2\) (Laplacian)

Identify and explain the physical meaning of each operator.

Exercise 4: Material Derivative and Pollutant Transport ⭐

Consider a homogeneous, incompressible fluid with velocity \(\mathbf{u}= (\alpha x, -\alpha y, 0)\) where \(\alpha\) is constant. The concentration of a pollutant in the fluid is: \[ c(x,y,z,t) = c_0 + \gamma t + \beta x^2 y e^{-\alpha t} \]

where \(c_0\), \(\gamma\), and \(\beta\) are constants.

  1. What is the (time) rate of change of concentration at a particular point (Eulerian perspective)?

  2. How does the concentration change following fluid parcels (Lagrangian perspective)?

  3. Consider a parcel at point \((1, 1, 0)\) at \(t = 0\). What is its concentration at \(t = 1\)?

Exercise 5: Commutation of Derivatives

Show that: \[ \frac{\partial}{\partial x}\left(\frac{Df}{Dt}\right) = \frac{D}{Dt}\left(\frac{\partial f}{\partial x}\right) + \frac{\partial \mathbf{u}}{\partial x} \cdot \nabla f \]

Clarify the meaning of the final term using index notation.

Hint: Expand \(D/Dt\) and use the product rule carefully.

Exercise 6: Incompressibility Condition

Is the velocity field \(\mathbf{u}= (\cos x, xz^2, z\sin x)\) kinematically admissible for an incompressible fluid? Justify your answer.

Exercise 7: Pressure Field in Inviscid Flow ⭐

Find the pressure field in the inviscid, incompressible flow with velocity field \(\mathbf{u}= (kx, -ky, 0)\), where \(k\) is a constant. Assume the flow is steady and the only body force is gravity.

Hint: Use the Euler equation (Navier-Stokes without viscosity).

Exercise 8: Cylindrical Polar Coordinates ⭐

Show that the continuity equation in cylindrical polar coordinates \((r, \theta, z)\) is: \[ \frac{1}{r}\frac{\partial(ru_r)}{\partial r} + \frac{1}{r}\frac{\partial u_\theta}{\partial \theta} + \frac{\partial u_z}{\partial z} = 0 \]

where \(\mathbf{U} = (u_r, u_\theta, u_z)\) is the velocity in cylindrical coordinates.

Hint: Transform from Cartesian coordinates using the chain rule.

Exercise 9: Boussinesq Approximation Application

Consider mantle convection where the temperature difference between the surface and CMB is \(\Delta T = 3000\) K, thermal expansivity \(\alpha = 3 \times 10^{-5}\) K\(^{-1}\), and reference density \(\rho_0 = 4000\) kg m\(^{-3}\).

  1. Calculate the density variation \(\Delta \rho = \rho_0 \alpha \Delta T\)

  2. What percentage of the reference density is this?

  3. Is the Boussinesq approximation valid for mantle convection?

  4. Under what conditions might the Boussinesq approximation break down in planetary interiors?

Exercise 10: Reynolds Number Estimation

Calculate the Reynolds number for the following geodynamic scenarios:

  1. Mantle convection: \(U = 50\) mm/yr, \(L = 2900\) km, \(\nu = 10^{17}\) m\(^2\)/s

  2. Magma flow in a dike: \(U = 1\) m/s, \(L = 1\) m, \(\nu = 10^3\) m\(^2\)/s

  3. Core convection: \(U = 0.1\) mm/s, \(L = 3000\) km, \(\nu = 10^{-6}\) m\(^2\)/s

Comment on whether inertial effects are important in each case.

Solutions to Selected Exercises

The material derivative is: \[ \frac{D\mathbf{f}}{Dt} = \frac{\partial \mathbf{f}}{\partial t} + (\mathbf{u}\cdot \nabla)\mathbf{f} \]

For the advective term to equal the total derivative, we need: \[ \frac{\partial \mathbf{f}}{\partial t} = 0 \]

Physical meaning: This occurs when the flow is steady (time-independent at fixed points in space). In steady flow, all changes experienced by a fluid parcel are due to its motion through spatially varying fields, not due to temporal changes at fixed locations.

In hydrostatic equilibrium: \(\frac{dp}{dz} = -\rho g\)

For an ideal gas: \(p = \rho RT\) where \(T\) is constant.

Therefore: \(\rho = \frac{p}{RT}\)

Substituting into the hydrostatic equation: \[ \frac{dp}{dz} = -\frac{pg}{RT} = -\frac{p}{H_s} \]

where \(H_s = RT/g\) is the scale height.

Integrating: \[ p(z) = p(0)\exp\left(-\frac{z}{H_s}\right) \]

And from the equation of state: \[ \rho(z) = \frac{p(z)}{RT} = \rho(0)\exp\left(-\frac{z}{H_s}\right) \]

For the troposphere, take \(T \approx 250\) K: \[ H_s = \frac{287 \times 250}{9.81} \approx 7.3 \text{ km} \]

This is indeed on the order of 8 km.

Divergence (scalar): \[ \nabla\cdot\mathbf{u}= \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} \] Index notation: \(\nabla_k u_k\) or \(\partial u_k/\partial x_k\)

Physical meaning: Rate of volume expansion per unit volume.

Curl (vector): \[ \nabla\times\mathbf{u}= \left(\frac{\partial w}{\partial y} - \frac{\partial v}{\partial z}, \frac{\partial u}{\partial z} - \frac{\partial w}{\partial x}, \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}\right) \] Index notation: \(\varepsilon_{ijk}\partial u_k/\partial x_j\)

Physical meaning: Local rotation or vorticity of the fluid.

Directional derivative (scalar operator): \[ \mathbf{u}\cdot \nabla= u\frac{\partial}{\partial x} + v\frac{\partial}{\partial y} + w\frac{\partial}{\partial z} \] Index notation: \(u_k \partial/\partial x_k\)

Physical meaning: Rate of change along the direction of velocity.

Laplacian (scalar): \[ \nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \] Index notation: \(\partial^2/\partial x_k \partial x_k\)

Physical meaning: Measure of how much a field differs from its average over neighboring points.

Given: \(\mathbf{u}= (\alpha x, -\alpha y, 0)\) and \(c = c_0 + \gamma t + \beta x^2 y e^{-\alpha t}\)

(a) Eulerian rate of change (at fixed point): \[ \frac{\partial c}{\partial t} = \gamma - \alpha\beta x^2 y e^{-\alpha t} \]

(b) Lagrangian rate of change (following parcels): \[ \frac{Dc}{Dt} = \frac{\partial c}{\partial t} + \mathbf{u}\cdot \nabla c \]

First calculate \(\nabla c\): \[ \nabla c = \left(2\beta xy e^{-\alpha t}, \beta x^2 e^{-\alpha t}, 0\right) \]

Then: \[ \begin{aligned} \mathbf{u}\cdot \nabla c &= \alpha x(2\beta xy e^{-\alpha t}) + (-\alpha y)(\beta x^2 e^{-\alpha t}) \\ &= 2\alpha\beta x^2y e^{-\alpha t} - \alpha\beta x^2 y e^{-\alpha t} \\ &= \alpha\beta x^2 y e^{-\alpha t} \end{aligned} \]

Therefore: \[ \frac{Dc}{Dt} = \gamma - \alpha\beta x^2 y e^{-\alpha t} + \alpha\beta x^2 y e^{-\alpha t} = \gamma \]

The concentration increases at constant rate \(\gamma\) following fluid parcels!

(c) Concentration at \(t = 1\):

First find the particle path. The trajectory satisfies: \[ \frac{dx}{dt} = \alpha x, \quad \frac{dy}{dt} = -\alpha y \]

Solutions: \(x(t) = x_0 e^{\alpha t}\), \(y(t) = y_0 e^{-\alpha t}\)

Starting at \((x_0, y_0, 0) = (1, 1, 0)\) at \(t = 0\): - At \(t = 1\): \(x(1) = e^\alpha\), \(y(1) = e^{-\alpha}\)

Since \(Dc/Dt = \gamma\): \[ c(t) = c(0) + \int_0^t \gamma \, dt' = c(0) + \gamma t \]

Initial concentration: \(c(0) = c_0 + \beta(1)^2(1)e^0 = c_0 + \beta\)

Final concentration: \(c(1) = c_0 + \beta + \gamma\)

Given: \(\mathbf{u}= (kx, -ky, 0)\) (steady, inviscid, incompressible)

Check incompressibility: \[ \nabla\cdot\mathbf{u}= k - k + 0 = 0 \quad \checkmark \]

Calculate the advective acceleration: \[ (\mathbf{u}\cdot \nabla)\mathbf{u}= \left(kx\frac{\partial}{\partial x} - ky\frac{\partial}{\partial y}\right)(kx, -ky, 0) \]

\[ = (k^2 x, k^2 y, 0) = \nabla\left(\frac{k^2x^2}{2} + \frac{k^2y^2}{2}\right) \]

Euler equation (inviscid Navier-Stokes): \[ (\mathbf{u}\cdot \nabla)\mathbf{u}= -\frac{1}{\rho}\nabla p - g\hat{\mathbf{k}} \]

Therefore: \[ \nabla\left(\frac{k^2x^2}{2} + \frac{k^2y^2}{2}\right) = -\frac{1}{\rho}\nabla p - \nabla(gz) \]

Integrating: \[ p = p_0 - \rho\left(\frac{k^2x^2}{2} + \frac{k^2y^2}{2} + gz\right) \]

where \(p_0\) is the pressure at the origin.

Transform from Cartesian \((x, y, z)\) to cylindrical \((r, \theta, z)\): \[ x = r\cos\theta, \quad y = r\sin\theta, \quad z = z \]

Velocity transformation: \[ u_r = u\cos\theta + v\sin\theta, \quad u_\theta = -u\sin\theta + v\cos\theta, \quad u_z = w \]

Using the chain rule: \[ \frac{\partial}{\partial x} = \cos\theta\frac{\partial}{\partial r} - \frac{\sin\theta}{r}\frac{\partial}{\partial \theta} \]

\[ \frac{\partial}{\partial y} = \sin\theta\frac{\partial}{\partial r} + \frac{\cos\theta}{r}\frac{\partial}{\partial \theta} \]

Substitute into the Cartesian continuity equation \(\partial u/\partial x + \partial v/\partial y + \partial w/\partial z = 0\) and use the velocity transformations. After simplification using trigonometric identities: \[ \frac{1}{r}\frac{\partial(ru_r)}{\partial r} + \frac{1}{r}\frac{\partial u_\theta}{\partial \theta} + \frac{\partial u_z}{\partial z} = 0 \]