Stokes and continuity equations in 2D

Stokes

The equation for the conservation of mass, for an incompressible fluid with negligible inertial forces (in comparison to gravitational acceleration), i.e. the Stokes equation, is

\[\frac{\partial \sigma'_{ij}}{\partial x_j} - \frac{\partial P}{\partial x_i} = -\rho g_i\]

In 2D, this gives us two equations,

(1)\[\begin{split}\begin{align} \frac{\partial \sigma'_{xx}}{\partial x} + \frac{\partial \sigma'_{xz}}{\partial z} - \frac{\partial P}{\partial x} & = -\rho g_x \\ \frac{\partial \sigma'_{zz}}{\partial z} + \frac{\partial \sigma'_{zx}}{\partial x} - \frac{\partial P}{\partial z} & = -\rho g_z \end{align}\end{split}\]

In the case of incompressible fluid a rheological law relates the deviatoric stress to strain rate:

(2)\[\sigma'_{ij} = 2\eta\dot\epsilon_{ij}\]

This is called the Newtonian viscous rheology, where the deviatoric stress is linearly dependent on the strain rate. Non-Newtonian rheologies, on the other hand, relate these with an exponent \((\sigma'_{ij})^n = 2\eta\dot\epsilon_{ij}\).

The strain rate is defined as

(3)\[\dot\epsilon_{ij} = \frac{1}{2}\left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right)\]

Applying (3) in (2) gives

\[\sigma'_{ij} = \eta\left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right)\]

which can be plugged in the Stokes equations (1) to get

(4)\[\begin{split}\begin{align} \frac{\partial}{\partial x} \left( \eta\left( \frac{\partial v_x}{\partial x} + \frac{\partial v_x}{\partial x} \right) \right) + \frac{\partial}{\partial z} \left( \eta\left( \frac{\partial v_x}{\partial z} + \frac{\partial v_z}{\partial x} \right) \right) - \frac{\partial P}{\partial x} & = -\rho g_x \\ \frac{\partial}{\partial z} \left( \eta\left( \frac{\partial v_z}{\partial z} + \frac{\partial v_z}{\partial z} \right) \right) + \frac{\partial}{\partial x} \left( \eta\left( \frac{\partial v_z}{\partial x} + \frac{\partial v_x}{\partial z} \right) \right) - \frac{\partial P}{\partial z} & = -\rho g_z \end{align}\end{split}\]

If viscosity is constant in respect to \(x\) and \(z\), this simplifies to

(5)\[\begin{split}\begin{align} \eta\left( \frac{\partial^2 v_x}{\partial x^2} + \frac{\partial^2 v_x}{\partial z^2} \right) - \frac{\partial P}{\partial x} & = -\rho g_x \\ \eta\left( \frac{\partial^2 v_z}{\partial z^2} + \frac{\partial^2 v_z}{\partial x^2} \right) - \frac{\partial P}{\partial z} & = -\rho g_z \end{align}\end{split}\]

As can be seen here, the Stokes equation relates fluid velocity to pressure gradients and body forces (gravitational forces); the flow is driven by either one or both.

Exercise

In what geological process would the fluid (rock) flow be driven by

  • gravitational forces
  • pressure gradients?

Exercise

Discretize equation (5).

  • What are the unknowns in the equation?
  • Can you rearrange the equations to solve for those unknowns explicitly?

Continuity

The continuity equation states that for an infinitesimally small volume of material, the incoming flux of new material is balanced by same amount of material going out. This applies for incompressible fluids.

(6)\[\frac{\partial v_x}{\partial x} + \frac{\partial v_z}{\partial z} = 0\]

Exercise

Discretize equation (6).

  • What are the unknowns in the equation?
  • Can you rearrange the equation to solve for those unknowns explicitly?

Stokes and continuity

In two-dimensional problems we have two Stokes equations (so called x-Stokes and z-Stokes here), and one continuity equation. So for each grid point we can formulate three equations. With these three equations we can solve for the three unknowns using an implicit formulation.

Note that the Stokes/continuity equations do not include time: they are used to find out the velocity field caused by density differences and external forcing (pressure gradients). If one knows the velocity and has chosen a \(\Delta t\), one can calculate where the material would advect during the next time step.

Exercise

Looking at the two Stokes equations (5) and the continuity equation (6), and their discretized versions,

  • how many boundary conditions do you need to solve them?
  • what would the coefficients for the implicit solution look like (matrix \(\mathbf{M}\))? How large would \(\mathbf{M}\) be?

Boundary conditions

Similarly to a constant temperature boundary condition for the heat equation, one can define constant velocity boundary conditions for the Stokes equation

\[\begin{split}\begin{align} v_x = \textrm{constant} & \quad \textrm{and} \\ v_z = \textrm{constant} & \end{align}\end{split}\]

This is so called “no-slip” boundary condition. Another often used boundary condition is so called free-slip. E.g. at boundary \(x=x_0\) this would be stated as

\[\begin{split}\begin{align} v_x = 0 & \quad \textrm{and} \\ \frac{\partial v_z}{\partial x} = 0 \end{align}\end{split}\]

If the velocity perpendicular to the boundary does not change across the boundary, the boundary condition is so called “stress free” boundary condition:

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

This would permit the outflow and inflow of material across the boundary.