Navier Stokes - flow around cylinder

---

Formulation

In this benchmark example we will show how to compute the flow of a Newtonian fluid around a cylinder. Even though the inlet velocity profile is constant, the resulting velocity field may not be stationary. We will see a pattern called Kármán vortex street.

The following figure shows the computational mesh used.

The strong formulation of the time dependent Navier Stokes equations for this problem is as follows

\begin{align} \frac{\partial\mathbf{u}}{\partial t} + \mathbf{u}\cdot\nabla\mathbf{u} + \nabla p - \nu\Delta\mathbf{u} &= 0\quad\text{in }\,\Omega\times (0,T)\\ \nabla\cdot\mathbf{u}&=0\quad\text{in }\,\Omega\times (0,T)\\ \mathbf{u}(t=0)&=\mathbf{0}\quad\text{in }\,\Omega\\ \mathbf{u} &= \mathbf{u}_{\text{in}}\quad\text{on }\,\Gamma_{\text{in}}\times (0,T)\\ \mathbf{u} &= \mathbf{0}\quad\text{on }\,\Gamma_{\text{wall}}\times (0,T)\\ \mathbf{u} &= \mathbf{0}\quad\text{on }\,\Gamma_{\text{cylinder}}\times (0,T)\\ (-p \mathbb{I} + \nu \nabla \mathbf{u}) \mathbf{n} &= \mathbf{0}\quad\text{on }\,\Gamma_{\text{out}}\times (0,T), \end{align}

where \(\mathbf{u}_{\text{in}}\) is defined as

\[\begin{split}\mathbf{u}_{\text{in}} = v_{\text{avg}}\begin{pmatrix} \frac{4\,y\,(D-y)}{D^2}\\ 0 \end{pmatrix}\end{split}\]

Taking test functions \(q\) and \(\mathbf{v}\) for pressure and velocity respectively, we multiply the first equation by \(\mathbf{v}\) and the second equation by \(q\) and integrate over \(\Omega\). Then we use the integration by parts to obtain the weak formulation.

\[\int_\Omega\frac{\partial\mathbf{u}}{\partial t}\cdot\mathbf{v}\,\text{d}x+\int_\Omega(\mathbf{u}\cdot\nabla\mathbf{u})\cdot\mathbf{v}\,\text{d}x+\nu\int_\Omega\nabla\mathbf{u}\cdot\nabla\mathbf{v}\,\text{d}x-\int_\Omega p\,\text{div}\,\mathbf{v}\,\text{d}x-\int_\Omega q\,\text{div}\,\mathbf{u}\,\text{d}x = 0\]

Compared to the Stokes equations, the Navier Stokes equations include a convection term \(\int_\Omega(\mathbf{u}\cdot\nabla\mathbf{u})\cdot\mathbf{v}\,\text{d}x\). Therefore, a nonlinear solver has to be used instead of a linear solver.

Another new term is the time derivative \(\frac{\partial \mathbf{u}}{\partial t}\) which we will treat using a Crank-Nicholson time stepping scheme. First we discretize the time derivative.

\begin{equation} \frac{\partial \mathbf{u}}{\partial t}(t) \approx \frac{\mathbf{u}(t)-\mathbf{u}(t-\delta t)}{\delta t} \end{equation}

In order to use the Crank-Nicolson scheme, we make a linear combination of explicit and implicit scheme using a factor \(\theta\) in part of the equation (here we use simpler notation \(\mathbf{u}^{i} = \mathbf{u}(t)\), \(\mathbf{u}^{i-1} = \mathbf{u}(t-\delta t)\)). Notice that the last two terms are not part of the linear combination since we need to use the implicit scheme on them to have stability.

\begin{equation} \begin{split} \int_\Omega\frac{\mathbf{u}^i-\mathbf{u}^{i-1}}{\delta t}\cdot\mathbf{v}\,\text{d}x\\ +\theta\left(\int_\Omega(\mathbf{u}^i\cdot\nabla\mathbf{u}^i)\cdot\mathbf{v}\,\text{d}x+\nu\int_\Omega\nabla\mathbf{u}^i\cdot\nabla\mathbf{v}\,\text{d}x\right)\\+(1-\theta)\left(\int_\Omega(\mathbf{u}^{i-1}\cdot\nabla\mathbf{u}^{i-1})\cdot\mathbf{v}\,\text{d}x\\+\nu\int_\Omega\nabla\mathbf{u}^{i-1}\cdot\nabla\mathbf{v}\,\text{d}x\right)\\ -\int_\Omega p^i\,\text{div}\,\mathbf{v}\,\text{d}x-\int_\Omega q\,\text{div}\,\mathbf{u}^i\,\text{d}x = 0 \end{split} \end{equation}

Implementation

First, we need to import all neccessary packages.

import dolfin as df
from time import time

Then we read mesh and boundary markers from the file.

mesh = df.Mesh()
hdf = df.HDF5File(mesh.mpi_comm(), "bench_csg.h5", "r")
hdf.read(mesh, "/mesh", False)

bndry = df.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
hdf.read(bndry, "/boundaries")

with df.XDMFFile("results/boundaries.xdmf") as xdmf:
    xdmf.write(bndry)

Next, we define the function spaces V, P and create the mixed function space W. The MINI element is used in this example.

dim = mesh.geometry().dim()
U = df.FiniteElement("CG", mesh.ufl_cell(), 1)
B = df.FiniteElement("Bubble", mesh.ufl_cell(), dim + 1)
P = df.FiniteElement("CG", mesh.ufl_cell(), 1)
V = df.VectorElement(df.NodalEnrichedElement(U, B))
W = df.FunctionSpace(mesh, df.MixedElement([V, P]))

Let us define some parameters and constants such are the average velocity at the inlet, kinematic viscosity, timestep size and time T.

velocity = 1.5
nu = df.Constant(0.001)
dt = 0.1
t_end = 15

Let us define the boundary conditions. The boundary markings are inlet: 1, outlet: 2, wall: 3, cylinder: 5. We prescribe zero vector on the wall and the cylinder using Dirichlet boundary conditions. Inlet velocity can be defined using an expression object.

inlet_expression = df.Expression(("v * 4.0 * x[1] * (D - x[1]) / pow(D,2)", "0.0"), D=0.41, v=velocity, degree=2)
bc_in = df.DirichletBC(W.sub(0), inlet_expression, bndry, 1)

zero_vector = df.Constant((0.0, 0.0))
bc_walls = df.DirichletBC(W.sub(0), zero_vector, bndry, 3)
bc_sphere = df.DirichletBC(W.sub(0), zero_vector, bndry, 5)

bcs = [bc_in, bc_walls, bc_sphere]

Now we define the variational form for functions \(\mathbf{u}, p\) with help of test functions \(\mathbf{v}, q\) and functions \(\mathbf{u}_0, p_0\) representing the velocity and pressure fields from previous timestep.

v, q = df.TestFunctions(W)
w = df.Function(W)
w0 = df.Function(W)
u, p = df.split(w)
u0, p0 = df.split(w0)

a = lambda u, v: df.inner(df.grad(u)*u, v)*df.dx + nu*df.inner(df.grad(u), df.grad(v))*df.dx
b = lambda q, v: q*df.div(v)*df.dx

F1 = a(u, v) - b(p, v) - b(q, u)
F0 = a(u0, v) - b(p, v) - b(q, u)

theta = df.Constant(0.5)
F = df.Constant(1.0/dt)*df.inner((u-u0), v)*df.dx + theta*F1 + (1-theta)*F0

We create a nonlinear variational problem object using our variational form and a nonlinear variational solver. Optionally, we can set up some of the solver parameters such as absolute and relative tolerance for convergence.

problem = df.NonlinearVariationalProblem(F, w, bcs, J=df.derivative(F, w))
solver = df.NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['linear_solver'] = 'mumps'
solver.parameters['newton_solver']['absolute_tolerance'] = 1e-12
solver.parameters['newton_solver']['relative_tolerance'] = 1e-12

Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.

Setup XDMF files to save results for viewing in Paraview.

ufile = df.XDMFFile("results/u.xdmf")
pfile = df.XDMFFile("results/p.xdmf")
ufile.parameters["flush_output"] = True
pfile.parameters["flush_output"] = True

Since we have time dependent problem, we need to solve the problem for each timestep using a loop. We start at time \(t=0\), then for each timestep, we assign the previous solution \(w\) to \(w_0\) and solve for the new \(w\). It is also a good idea to save the solution at each timestep into the XDMF file.

tick = time()
t = 0.0
(u, p) = w.split(deepcopy=True)
u.rename("v", "velocity")
p.rename("p", "pressure")
df.assign(u, w.sub(0))
df.assign(p, w.sub(1))
ufile.write(u, t)
pfile.write(p, t)
while t < t_end:
    w0.assign(w)
    t += dt
    print("t={:.1f}s".format(t))
    solver.solve()
    df.assign(u, w.sub(0))
    df.assign(p, w.sub(1))
    ufile.write(u, t)
    pfile.write(p, t)
print("ellapsed = ", time() - tick, "s")

Complete Code

Run the code using the command python3 NavierStokes2D.py.

import dolfin as df
from time import time

mesh = df.Mesh()
hdf = df.HDF5File(mesh.mpi_comm(), "bench_csg.h5", "r")
hdf.read(mesh, "/mesh", False)

bndry = df.MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
hdf.read(bndry, "/boundaries")

with df.XDMFFile("results/boundaries.xdmf") as xdmf:
    xdmf.write(bndry)

dim = mesh.geometry().dim()
U = df.FiniteElement("CG", mesh.ufl_cell(), 1)
B = df.FiniteElement("Bubble", mesh.ufl_cell(), dim + 1)
P = df.FiniteElement("CG", mesh.ufl_cell(), 1)
V = df.VectorElement(df.NodalEnrichedElement(U, B))
W = df.FunctionSpace(mesh, df.MixedElement([V, P]))

velocity = 1.5
nu = df.Constant(0.001)
dt = 0.1
t_end = 15

inlet_expression = df.Expression(("v * 4.0 * x[1] * (D - x[1]) / pow(D,2)", "0.0"), D=0.41, v=velocity, degree=2)
bc_in = df.DirichletBC(W.sub(0), inlet_expression, bndry, 1)

zero_vector = df.Constant((0.0, 0.0))
bc_walls = df.DirichletBC(W.sub(0), zero_vector, bndry, 3)
bc_sphere = df.DirichletBC(W.sub(0), zero_vector, bndry, 5)

bcs = [bc_in, bc_walls, bc_sphere]

v, q = df.TestFunctions(W)
w = df.Function(W)
w0 = df.Function(W)
u, p = df.split(w)
u0, p0 = df.split(w0)

a = lambda u, v: df.inner(df.grad(u)*u, v)*df.dx + nu*df.inner(df.grad(u), df.grad(v))*df.dx
b = lambda q, v: q*df.div(v)*df.dx

F1 = a(u, v) - b(p, v) - b(q, u)
F0 = a(u0, v) - b(p, v) - b(q, u)

theta = df.Constant(0.5)
F = df.Constant(1.0/dt)*df.inner((u-u0), v)*df.dx + theta*F1 + (1-theta)*F0
J = df.derivative(F,w)

problem = df.NonlinearVariationalProblem(F, w, bcs, J)
solver = df.NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['linear_solver'] = 'mumps'
solver.parameters['newton_solver']['absolute_tolerance'] = 1e-12
solver.parameters['newton_solver']['relative_tolerance'] = 1e-12

ufile = df.XDMFFile("results/u.xdmf")
pfile = df.XDMFFile("results/p.xdmf")
ufile.parameters["flush_output"] = True
pfile.parameters["flush_output"] = True

tick = time()
t = 0.0
(u, p) = w.split(True)
u.rename("v", "velocity")
p.rename("p", "pressure")
df.assign(u, w.sub(0))
df.assign(p, w.sub(1))
ufile.write(u, t)
pfile.write(p, t)
while t < t_end:
    w0.assign(w)
    t += dt
    print("t={:.1f}s".format(t))
    solver.solve()
    df.assign(u, w.sub(0))
    df.assign(p, w.sub(1))
    ufile.write(u, t)
    pfile.write(p, t)
print("ellapsed = ", time() - tick, "s")