doc.rst

Navier-Stokes equations

Goals

Learn how to deal with mixed finite elements. Remember how fragile can numerical solutions be. Reproduce some cool physics -- Kármán vortex street.

Stokes flow around cylinder

Solve the following linear system of PDEs

$$- \nu\Delta\mathbf{u} + \nabla p &= \mathbf{0} &&\quad\text{ in }\Omega,$$$$\operatorname{div}\mathbf{u} &= 0 &&\quad\text{ in }\Omega,$$$$\mathbf{u} &= 0 &&\quad\text{ on }\Gamma_\mathrm{D},$$$$\mathbf{u} &= \mathbf{u}_\mathrm{IN} &&\quad\text{ on }\Gamma_\mathrm{IN},$$$$\nu\tfrac{\partial\mathbf{u}}{\partial\mathbf{n}} -p\mathbf{n} &= 0 &&\quad\text{ on }\Gamma_\mathrm{N}$$

using FE discretization with data

$$\Omega &= (0, 2.2)\times(0, 0.41) - B_{0.05}\left((0.2,0.2)\right),$$$$\Gamma_\mathrm{N} &= \left\{ x = 2.2 \right\} = \text{(green)},$$$$\Gamma_\mathrm{IN} &= \left\{ x = 0.0 \right\} = \text{(red)},$$$$\Gamma_\mathrm{D} &= \partial\Omega\setminus(\Gamma_\mathrm{N}\cup\Gamma_\mathrm{IN}) = \text{(black)},$$$$u_\mathrm{IN} &= \left( \frac{4Uy (0.41-y)}{0.41^2} , 0 \right),$$$$\nu &= 0.001, \qquad U = 0.3$$

where B_R(z) is a disc of radius R and center z

Task 1

Write the weak formulation of the problem and a spatial discretization by a mixed finite element method.

Task 2

Build a mesh, prepare a mesh function marking Γ_IN, Γ_N and Γ_D and plot it to check its correctness.

Hint

Use the FEniCS meshing tool mshr, see mshr documentation.

from dolfin import *
import mshr

# Discretization parameters
N_circle = 16
N_bulk = 64

# Define domain
center = Point(0.2, 0.2)
radius = 0.05
L = 2.2
W = 0.41
geometry =  mshr.Rectangle(Point(0.0, 0.0), Point(L, W)) \
           -mshr.Circle(center, radius, N_circle)

# Build mesh
mesh = mshr.generate_mesh(geometry, N_bulk)

Hint

Try yet another way to mark the boundaries by direct access to the mesh entities by vertices(mesh) <dolfin.cpp.mesh.vertices>, facets(mesh) <dolfin.cpp.mesh.facets>, cells(mesh) <dolfin.cpp.mesh.cells> mesh-entity iterators:

# Construct facet markers
bndry = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
for f in facets(mesh):
    mp = f.midpoint()
    if near(mp[0], 0.0):  # inflow
        bndry[f] = 1
    elif near(mp[0], L):  # outflow
        bndry[f] = 2
    elif near(mp[1], 0.0) or near(mp[1], W):  # walls
        bndry[f] = 3
    elif mp.distance(center) <= radius:  # cylinder
        bndry[f] = 5

# Dump facet markers to file to plot in Paraview
with XDMFFile('facets.xdmf') as f:
    f.write(bndry)

Task 3

Construct the mixed finite element space and the bilinear and linear forms together with appropriate DirichletBC <dolfin.fem.bcs.DirichletBC> object.

Hint

Use for example the stable Taylor-Hood finite elements:

# Build function spaces (Taylor-Hood)
P2 = VectorElement("P", mesh.ufl_cell(), 2)
P1 = FiniteElement("P", mesh.ufl_cell(), 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)

Hint

To define Dirichlet BC on subspace use the W.sub() <dolfin.functions.functionspace.FunctionSpace.sub> method:

bc_walls = DirichletBC(W.sub(0), (0, 0), bndry, 3)

Hint

To build the forms use:

# Define trial and test functions
u, p = TrialFunctions(W)
v, q = TestFunctions(W)

Then you can define forms on mixed space using u, p, v, q as usual.

Steady Navier-Stokes flow

Task 4

Modify the problem into the Navier-Stokes equations given by

 − νΔu + u ⋅ ∇u + ∇p = 0  in Ω

together with stokes₂--stokes₅. Compute the DFG-flow around cylinder benchmark 2D-1, laminar case, Re=20 given by navierstokes, stokes₂--stokes₅, dfg20.

Hint

As usual get rid of TrialFunctions in favour of nonlinear dependence on Function. You can split a Function on a mixed space into components:

w = Function(W)
u, p = split(w)

F = nu*inner(grad(u), grad(v))*dx + ...

Task 5

Add computation of lift and drag coefficients C_D, C_L and pressure difference p_diff as defined on the DFG 2D-1 website.

Hint

Use assemble <dolfin.fem.assembling.assemble> function to evaluate the lift and drag functionals.

Use either Function.split() <dolfin.functions.function.Function.split> or Function.sub() <dolfin.functions.function.Function.sub> to extract pressure p from solution w for evaluation. Evaluate the pressure p at point a = Point(234, 567) by calling p(a).

Task 6

Check computed pressure difference and lift/drag coefficents against the reference. Investigate if/how the lift coefficent is sensitive to changes in the discretization parameters --conduct a convergence study.

Kármán vortex street

Task 7

Consider evolutionary Navier-Stokes equations

u_t − νΔu + u ⋅ ∇u + ∇p = 0.

Prepare temporal discretization using the Crank-Nicolson scheme <theta-table> to compute a solution of unsteady, stokes₂--stokes₅, dfg20 on time interval (0, 8) but use

U = 1

instead of dfg20_6b. Plot the transient solution.

pub

Reference solution

Note

You can run FEniCS codes in parallel (using MPI) by

mpirun -n <np> python3 <yourscript>.py

where for <np> substitute number of processors to use.

To benefit from parallism you can run the unsteady Navier-Stokes part of the code below on, say, eight cores:

mpirun -n 8 python3 -c"import dfg; dfg.task_7()"

Download Code <dfg.py>

dfg.py