Notebook 5: Stokes Equation

Binder

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

#|  echo: false  # Hide in html version

# This is required to fix pyvista 
# (visualisation) crashes in interactive notebooks (including on binder)

import nest_asyncio
nest_asyncio.apply()

#| output: false # Suppress warnings in html version

import underworld3 as uw
import numpy as np
import sympy

uw.__file__

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 = meshball.Gamma / sympy.sqrt(meshball.Gamma.dot(meshball.Gamma))

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.sym * 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"

v_theta_fn_xy = r * meshball.CoordinateSystem.rRotN.T * sympy.Matrix((0,1))
v_theta_fn_xy

# 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)

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")

penalty = 1e5
Gamma_N = meshball.Gamma / sympy.sqrt(meshball.Gamma.dot(meshball.Gamma))
stokes.add_natural_bc(penalty * Gamma_N.dot(v_soln.sym) *  Gamma_N, "Upper")

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

stokes.solve(verbose=False)

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)

dv = uw.function.evaluate(ns * v_theta_fn_xy, v_soln.coords).squeeze() / vnorm
v_soln.array[...] -= dv.reshape(v_soln.array.shape)

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

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")

#| fig-cap: "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"


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

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

stokes.view()

meshball.CoordinateSystem.spherical.coords