Navier-Stokes equations¶

Stokes flow around cylinder¶

Solve the following linear system of PDEs

(1)¶\[ \begin{align}\begin{aligned}- \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}\end{aligned}\end{align} \]

using FE discretization with data

(2)¶\[ \begin{align}\begin{aligned}\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\end{aligned}\end{align} \]

where \(B_R(\mathbf{z})\) is a disc of radius \(R\) and center \(\mathbf{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 \(\Gamma_\mathrm{IN}\), \(\Gamma_\mathrm{N}\) and \(\Gamma_\mathrm{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), facets(mesh), cells(mesh) 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 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() 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

(3)¶\[- \nu\Delta\mathbf{u} + \mathbf{u}\cdot\nabla\mathbf{u} + \nabla p = 0 \quad\text{ in }\Omega\]

together with (1)\(_2\)–(1)\(_5\). Compute the DFG-flow around cylinder benchmark 2D-1, laminar case, Re=20 given by (3), (1)\(_2\)–(1)\(_5\), (2).

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_\mathrm{D}\), \(C_\mathrm{L}\) and pressure difference \(p_\mathrm{diff}\) as defined on the DFG 2D-1 website.

Hint

Use assemble function to evaluate the lift and drag functionals.

Use either Function.split() or 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

(4)¶\[u_t - \nu\Delta\mathbf{u} + \mathbf{u}\cdot\nabla\mathbf{u} + \nabla p = 0.\]

Prepare temporal discretization using the Crank-Nicolson scheme to compute a solution of (4), (1)\(_2\)–(1)\(_5\), (2) on time interval \((0,8)\) but use

\[U = 1\]

instead of (2)\(_{6b}\). Plot the transient solution.

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

Show/Hide Code

Download Code

from dolfin import *
import mshr
import matplotlib.pyplot as plt


def build_space(N_circle, N_bulk, u_in):
    """Prepare data for DGF benchmark. Return function
    space, list of boundary conditions and surface measure
    on the cylinder."""

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

    # 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

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

    # Prepare Dirichlet boundary conditions
    bc_walls = DirichletBC(W.sub(0), (0, 0), bndry, 3)
    bc_cylinder = DirichletBC(W.sub(0), (0, 0), bndry, 5)
    bc_in = DirichletBC(W.sub(0), u_in, bndry, 1)
    bcs = [bc_cylinder, bc_walls, bc_in]

    # Prepare surface measure on cylinder
    ds_circle = Measure("ds", subdomain_data=bndry, subdomain_id=5)

    return W, bcs, ds_circle


def solve_stokes(W, nu, bcs):
    """Solve steady Stokes and return the solution"""

    # Define variational forms
    u, p = TrialFunctions(W)
    v, q = TestFunctions(W)
    a = nu*inner(grad(u), grad(v))*dx - p*div(v)*dx - q*div(u)*dx
    L = inner(Constant((0, 0)), v)*dx

    # Solve the problem
    w = Function(W)
    solve(a == L, w, bcs)

    return w


def solve_navier_stokes(W, nu, bcs):
    """Solve steady Navier-Stokes and return the solution"""

    # Define variational forms
    v, q = TestFunctions(W)
    w = Function(W)
    u, p = split(w)
    F = nu*inner(grad(u), grad(v))*dx + dot(dot(grad(u), u), v)*dx \
        - p*div(v)*dx - q*div(u)*dx

    # Solve the problem
    solve(F == 0, w, bcs)

    return w


def solve_unsteady_navier_stokes(W, nu, bcs, T, dt, theta):
    """Solver unsteady Navier-Stokes and write results
    to file"""

    # Current and old solution
    w = Function(W)
    u, p = split(w)

    w_old = Function(W)
    u_old, p_old = split(w_old)

    # Define variational forms
    v, q = TestFunctions(W)
    F = ( Constant(1/dt)*dot(u - u_old, v)
          + Constant(theta)*nu*inner(grad(u), grad(v))
          + Constant(theta)*dot(dot(grad(u), u), v)
          + Constant(1-theta)*nu*inner(grad(u), grad(v))
          + Constant(1-theta)*dot(dot(grad(u_old), u_old), v)
          - p*div(v)
          - q*div(u)
        )*dx
    J = derivative(F, w)

    # Create solver
    problem = NonlinearVariationalProblem(F, w, bcs, J)
    solver = NonlinearVariationalSolver(problem)
    solver.parameters['newton_solver']['linear_solver'] = 'mumps'

    f = XDMFFile('velocity_unteady_navier_stokes.xdmf')
    u, p = w.split()

    # Perform time-stepping
    t = 0
    while t < T:
        w_old.vector()[:] = w.vector()
        solver.solve()
        t += dt
        f.write(u, t)


def save_and_plot(w, name):
    """Saves and plots provided solution using the given
    name"""

    u, p = w.split()

    # Store to file
    with XDMFFile("results_{}/u.xdmf".format(name)) as f:
        f.write(u)
    with XDMFFile("results_{}/p.xdmf".format(name)) as f:
        f.write(p)

    # Plot
    plt.figure()
    pl = plot(u, title='velocity {}'.format(name))
    plt.colorbar(pl)
    plt.figure()
    pl = plot(p, mode='warp', title='pressure {}'.format(name))
    plt.colorbar(pl)


def postprocess(w, nu, ds_circle):
    """Return lift, drag and the pressure difference"""

    u, p = w.split()

    # Report drag and lift
    n = FacetNormal(w.function_space().mesh())
    force = -p*n + nu*dot(grad(u), n)
    F_D = assemble(-force[0]*ds_circle)
    F_L = assemble(-force[1]*ds_circle)

    U_mean = 0.2
    L = 0.1
    C_D = 2/(U_mean**2*L)*F_D
    C_L = 2/(U_mean**2*L)*F_L

    # Report pressure difference
    a_1 = Point(0.15, 0.2)
    a_2 = Point(0.25, 0.2)
    try:
        p_diff = p(a_1) - p(a_2)
    except RuntimeError:
        p_diff = 0

    return C_D, C_L, p_diff


def tasks_1_2_3_4():
    """Solve and plot alongside Stokes and Navier-Stokes"""

    # Problem data
    u_in = Expression(("4.0*U*x[1]*(0.41 - x[1])/(0.41*0.41)", "0.0"),
                      degree=2, U=0.3)
    nu = Constant(0.001)

    # Discretization parameters
    N_circle = 16
    N_bulk = 64

    # Prepare function space, BCs and measure on circle
    W, bcs, ds_circle = build_space(N_circle, N_bulk, u_in)

    # Solve Stokes
    w = solve_stokes(W, nu, bcs)
    save_and_plot(w, 'stokes')

    # Solve Navier-Stokes
    w = solve_navier_stokes(W, nu, bcs)
    save_and_plot(w, 'navier-stokes')

    # Open and hold plot windows
    plt.show()


def tasks_5_6():
    """Run convergence analysis of drag and lift"""

    # Problem data
    u_in = Expression(("4.0*U*x[1]*(0.41 - x[1])/(0.41*0.41)", "0.0"),
                      degree=2, U=0.3)
    nu = Constant(0.001)

    # Push log levelo to silence DOLFIN
    old_level = get_log_level()
    warning = LogLevel.WARNING if cpp.__version__ > '2017.2.0' else WARNING
    set_log_level(warning)

    fmt_header = "{:10s} | {:10s} | {:10s} | {:10s} | {:10s} | {:10s}"
    fmt_row = "{:10d} | {:10d} | {:10d} | {:10.4f} | {:10.4f} | {:10.6f}"

    # Print table header
    print(fmt_header.format("N_bulk", "N_circle", "#dofs", "C_D", "C_L", "p_diff"))

    # Solve on series of meshes
    for N_bulk in [32, 64, 128]:
        for N_circle in [N_bulk, 2*N_bulk, 4*N_bulk]:

            # Prepare function space, BCs and measure on circle
            W, bcs, ds_circle = build_space(N_circle, N_bulk, u_in)

            # Solve Navier-Stokes
            w = solve_navier_stokes(W, nu, bcs)

            # Compute drag, lift
            C_D, C_L, p_diff = postprocess(w, nu, ds_circle)
            print(fmt_row.format(N_bulk, N_circle, W.dim(), C_D, C_L, p_diff))

    # Pop log level
    set_log_level(old_level)


def task_7():
    """Solve unsteady Navier-Stokes to resolve
    Karman vortex street and save to file"""

    # Problem data
    u_in = Expression(("4.0*U*x[1]*(0.41 - x[1])/(0.41*0.41)", "0.0"),
                      degree=2, U=1)
    nu = Constant(0.001)
    T = 8

    # Discretization parameters
    N_circle = 16
    N_bulk = 64
    theta = 1/2
    dt = 0.2

    # Prepare function space, BCs and measure on circle
    W, bcs, ds_circle = build_space(N_circle, N_bulk, u_in)

    # Solve unsteady Navier-Stokes
    solve_unsteady_navier_stokes(W, nu, bcs, T, dt, theta)


if __name__ == "__main__":

    tasks_1_2_3_4()
    tasks_5_6()
    task_7()