Notebook 5: Stokes Equation

In this example, we solve the Stokes equation for buoyancy-driven, incompressible flow in a creeping, viscous fluid.

\[ -\nabla \cdot \left[ \boldsymbol{\tau} - p \mathbf{I} \right] = \rho \mathbf{g} \]

\[ \nabla \cdot \mathbf{u} = 0 \]

\(p\) is pressure. The constitutive law expresses the deviatoric stress ( \(\boldsymbol{\tau}\) ) in terms of the velocity (\(\mathbf{u}\)) gradients ( \(\nabla \mathbf{u}\) ) through a viscosity tensor, \(\eta\). \(p\) is pressure.

\[ \boldsymbol{\tau} = \frac{\eta}{2}\left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) \]

A classical problem in geodynamics, is to solve this equation for harmonic forcing at one specific depth (the Greens function or kernels) for which analytic solutions are available in various geometries (e.g. Parsons & Daly, 1981).

We have pre-built meshes in underworld3 that have embedded internal surfaces where we can apply the harmonic forcing. This allows us to compute the flow kernels. We’ll use a cylindrical annulus this time and we’ll use free-slip boundary conditions throughout

import underworld3 as uw
import numpy as np
import sympy

res = 0.075
r_o = 1.0
r_int = 0.825
r_i = 0.55

meshball = uw.meshing.AnnulusInternalBoundary(radiusOuter=r_o, 
                                              radiusInternal=r_int, 
                                              radiusInner=r_i, 
                                              cellSize_Inner=res,
                                              cellSize_Internal=res*0.5,
                                              cellSize_Outer=res,
                                              centre=False,)
# meshball.view()

x, y = meshball.CoordinateSystem.X
r, th = meshball.CoordinateSystem.xR
unit_rvec = meshball.CoordinateSystem.unit_e_0


# Orientation of surface normals
Gamma_N = unit_rvec

Solver setup

We can obtain unit vectors in the natural coordinate system (here \(r\), \(\theta\)) as mesh.CoordinateSystem.unit_e_0, mesh.CoordinateSystem.unit_e_1. There is a null space if we apply the boundary conditions exactly, and so we define a function to represent the null space.

We can set solver options via the petsc_options interface on the solver. PETSc provides an enormous variety of solver configurations and options that can be set. Their Guide to Stokes Equation is an illustration of the complexity of the problem and the difficulty in trying to set up a good one-size-fits-all default configuration, particularly when moving to parallel examples.

# Create a density structure / buoyancy force
# gravity will vary linearly from zero at the centre
# of the sphere to (say) 1 at the surface


# Null space in velocity (constant v_theta) expressed in x,y coordinates
v_theta_fn_xy = r * meshball.CoordinateSystem.rRotN.T * sympy.Matrix((0,1))

# Mesh variables for the unknowns

v_soln = uw.discretisation.MeshVariable("V0", meshball, 2, degree=2, varsymbol=r"{v_0}")
p_soln = uw.discretisation.MeshVariable("p", meshball, 1, degree=1, continuous=True)

stokes = uw.systems.Stokes(
    meshball, 
    velocityField=v_soln, 
    pressureField=p_soln,
)

stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.shear_viscosity_0 = 1
stokes.tolerance = 1.0e-6

stokes.petsc_options.setValue("ksp_monitor", None)
stokes.petsc_options.setValue("snes_monitor", None)
stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"

uw.systems.Stokes.view()

Stokes Equation Solver

This class provides functionality for a discrete representation of the Stokes flow equations assuming an incompressibility (or near-incompressibility) constraint.

\[ -\nabla \cdot \color{Blue}{\underbrace{\Bigl[ \boldsymbol{\tau} - p \mathbf{I} \Bigr]}_{\mathbf{F}}} = \color{Maroon}{\underbrace{\Bigl[ \mathbf{f} \Bigl] }_{\mathbf{h}}} \]

\[ \underbrace{\Bigl[ \nabla \cdot \mathbf{u} \Bigr]}_{\mathbf{h}_p} = 0 \]

The flux term is a deviatoric stress ( \(\boldsymbol{\tau}\) ) related to velocity gradients ( \(\nabla \mathbf{u}\) ) through a viscosity tensor, \(\eta\), and a volumetric (pressure) part \(p\)

\[ \mathbf{F}: \quad \boldsymbol{\tau} = \frac{\eta}{2}\left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) \]

The constraint equation, \(\mathbf{h}_p = 0\) gives incompressible flow by default but can be set to any function of the unknown \(\mathbf{u}\) and \(\nabla\cdot\mathbf{u}\)

Properties

• The unknowns are velocities \(\mathbf{u}\) and a pressure-like constraint parameter \(\mathbf{p}\)

• The viscosity tensor, \(\boldsymbol{\eta}\) is provided by setting the constitutive_model property to one of the scalar uw.constitutive_models classes and populating the parameters. It is usually a constant or a function of position / time and may also be non-linear or anisotropic.

• \(\mathbf f\) is a volumetric source term (i.e. body forces) and is set by providing the bodyforce property.

• An Augmented Lagrangian approach to application of the incompressibility constraint is to penalise incompressibility in the Stokes equation by adding $ $ when the weak form of the equations is constructed. (this is in addition to the constraint equation, unlike in the classical penalty method). This is activated by setting the penalty property to a non-zero floating point value which adds the term in the sympy expression.

• A preconditioner is usually required for the saddle point system and this is provided through the saddle_preconditioner property. The default choice is \(1/\eta\) for a scalar viscosity function.

Notes

• For problems with viscoelastic behaviour, the flux term contains the stress history as well as the stress and this term is a Lagrangian quantity that has to be tracked on a particle swarm.

• The interpolation order of the pressureField variable is used to determine the integration order of the mixed finite element method and is usually lower than the order of the velocityField variable.

• It is possible to set discontinuous pressure variables by setting the p_continous option to False

# The class instance has details of how the general class has been
# constructed for the particular case (boundary conditions, constitutive models etc)
# If class_documentation=True is used, then the description above is printed first.

stokes.view(class_documentation=False)

Class: <class ‘underworld3.systems.solvers.SNES_Stokes’>

Underworld / PETSc General Saddle Point Equation Solver

Primary problem:

<IPython.core.display.Latex object>

$ = 0 $

Constraint:

$ = 0 $

Where:

\(\quad\)\(\displaystyle \uplambda\)\(=\)\(\displaystyle 0\)

\(\quad\)\(\displaystyle \eta\)\(=\)\(\displaystyle 1\)

Boundary Conditions

Type	Boundary	Expression

This solver is formulated as a 2 dimensional problem with a 2 dimensional mesh

Natural boundary conditions

In FEM, natural boundary conditions are fluxes specified at a boundary (through surface integrals). We can also apply integrals to internal surfaces. We need to compute the vector components of the boundary condition (normal / tangential) and supply them in the Cartesian frame.

The bouyancy force on the internal surface is already radial, so this translates to

    stokes.add_natural_bc(-t_init * unit_rvec, "Internal")

To set a no-normal-flow boundary condition, we need to penalise the radial velocity at the boundaries. Symbolically this is

\[ \mathbf{f} = \lambda \left( \mathbf{v} \cdot \Gamma_N \right) \,\, \Gamma_N \]

where \(\lambda\) is a large penalty value. This translates into sympy code as

    stokes.add_natural_bc(
                penalty * Gamma_N.dot(v_soln.sym) *  Gamma_N, "Lower"
                )

stokes.bodyforce = sympy.Matrix([0,0])

# Note, the thermal bouyancy field is localised in the radius using a 
# gaussian solely for the purposes of plotting. It is automatically
# a delta function when applied through a surface integral

t_init = sympy.sin(5*th) * sympy.exp(-1000.0 * ((r - r_int) ** 2)) 

stokes.add_natural_bc(-t_init * unit_rvec, "Internal")
stokes.add_natural_bc(10000 * Gamma_N.dot(v_soln.sym) *  Gamma_N, "Upper")

if r_i != 0.0:
    stokes.add_natural_bc(10000 * Gamma_N.dot(v_soln.sym) *  Gamma_N, "Lower")

stokes.solve()

  0 SNES Function norm 2.307507502179e-01
    Residual norms for Solver_7_ solve.
    0 KSP Residual norm 8.742204758196e+00
    1 KSP Residual norm 7.602202411042e-07
  1 SNES Function norm 4.872437719294e-08

Removal of the null space

We can use the uw.maths.Integral to compute the inner product of the null space and the velocity solution. It is not zero, so we remove it !

# Null space evaluation

I0 = uw.maths.Integral(meshball, v_theta_fn_xy.dot(v_theta_fn_xy))
vnorm = I0.evaluate()
                                  
I0.fn = v_theta_fn_xy.dot(v_soln.sym)
ns = I0.evaluate()

print(ns/vnorm, vnorm)

with meshball.access(v_soln):
    dv = uw.function.evaluate(ns * v_theta_fn_xy, v_soln.coords) / vnorm
    v_soln.data[...] -= dv 

ns = I0.evaluate()
print(ns/vnorm, vnorm)

0.0005624988973658309 1.4249504223200373
3.7567961262338314e-19 1.4249504223200373

if uw.mpi.size == 1:
    import pyvista as pv
    import underworld3.visualisation as vis

    pvmesh = vis.mesh_to_pv_mesh(meshball)
    pvmesh.point_data["P"] = vis.scalar_fn_to_pv_points(pvmesh, p_soln.sym)
    pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v_soln.sym)
    pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, t_init)
    
    skip = 1
    points = np.zeros((meshball._centroids[::skip].shape[0], 3))
    points[:, 0] = meshball._centroids[::skip, 0]
    points[:, 1] = meshball._centroids[::skip, 1]
    point_cloud = pv.PolyData(points)

    pvstream = pvmesh.streamlines_from_source(
        point_cloud, vectors="V", 
        integration_direction="both", 
        integrator_type=45,
        surface_streamlines=True,
        initial_step_length=0.01,
        max_time=0.5,
        max_steps=500, 
    )
   

    pl = pv.Plotter(window_size=(750, 750))

    pl.add_mesh(
        pvmesh,
        cmap="coolwarm",
        edge_color="Grey",
        edge_opacity=0.33,
        scalars="T",
        show_edges=True,
        use_transparency=False,
        opacity=1.0,
        show_scalar_bar=False
    )


    pl.add_mesh(pvstream, opacity=0.3, show_scalar_bar=False, cmap="Greens", render_lines_as_tubes=False)

    pl.export_html("html5/stokes_annulus_plot.html")

from IPython.display import IFrame
IFrame(src="html5/stokes_annulus_plot.html", width=500, height=400)

Interactive Image: Annulus mesh of triangular elements on which we evaluated Stokes flow driven by an internal delta function buoyancy. Boundary conditions are free slip, imposed using a penalty on the radial velocity at the boundary

References

Parsons, B., & Daly, S. (1983). The relationship between surface topography, gravity anomalies, and temperature structure of convection. Journal of Geophysical Research, 88(B2), 1129. https://doi.org/10/frq6gv