Benchmark: Flow around a cylinder

23.3. Benchmark: Flow around a cylinder#

Popular benchmark examples are proposed by Schäfer and Turek: https://wwwold.mathematik.tu-dortmund.de/lsiii/cms/papers/SchaeferTurek1996.pdf

Find velocity $u : Ω \times [0, T] \to R^{d}$ and pressure $p : Ω \times [0, T] \to R$ such that

\begin{array}{r} \begin{array}{ccccl} \frac{\partial u}{\partial t} - ν Δ u + u \nabla u & + & \nabla p & = & f \\ div u & = & 0 \end{array} \end{array}

from ngsolve import *
from ngsolve.webgui import Draw

geometry for 2D testcase:

from netgen.occ import *

shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()
shape.edges.name="wall"
shape.edges.Min(X).name="inlet"
shape.edges.Max(X).name="outlet"
Draw (shape);

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.07)).Curve(3)
Draw (mesh);

Higher order Taylor-Hood element pairing:

V = VectorH1(mesh,order=3, dirichlet="wall|cyl|inlet")
Q = H1(mesh,order=2)
X = V*Q

u,p = X.TrialFunction()
v,q = X.TestFunction()

nu = 0.001  # viscosity
stokes = (nu*InnerProduct(grad(u), grad(v))+ \
    div(u)*q+div(v)*p - 1e-10*p*q)*dx

a = BilinearForm(stokes).Assemble()

# nothing here ...
f = LinearForm(X).Assemble()

# gridfunction for the solution
gfu = GridFunction(X)

solve Stokes equation for initial condition, parabolic inflow at inlet:

umean = 0.3
uin = CoefficientFunction( (umean*4*y*(0.41-y)/(0.41*0.41), 0) )
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

inv_stokes = a.mat.Inverse(X.FreeDofs())

res = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

Draw (gfu.components[0], mesh);

implicit/explicit time-stepping:

\frac{u_{n + 1} - u_{n}}{τ} - ν Δ u_{n + 1} + \nabla p_{n + 1} = f - u_{n} \nabla u_{n}

and

div u_{n + 1} = 0

See high order IMEX timestepping by Ascher, Ruuth and Spiteri: https://www.sciencedirect.com/science/article/pii/S0168927497000561

tau = 0.001 # timestep

mstar = BilinearForm(u*v*dx+tau*stokes).Assemble()
inv = mstar.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")

the non-linear convective term $\int u \nabla u v$

conv = BilinearForm(X, nonassemble=True)
conv += (Grad(u) * u) * v * dx

implicit Euler/explicit Euler splitting method:

t = 0; i = 0
tend = 5
vel = gfu.components[0]
scene = Draw (gfu.components[0], mesh, min=0, max=0.4, autoscale=False)

with TaskManager():
    while t < tend:
        res = conv.Apply(gfu.vec) + a.mat*gfu.vec
        gfu.vec.data -= tau * inv * res    

        t = t + tau; i = i + 1
        if i%10 == 0: 
            scene.Redraw()
            # print(f"t = {t}", end='\r')

23.3.1. Exercises#

Try other examples of the benchmark (2D and 3D)
Evaluate drag and lift coefficients (see https://docu.ngsolve.org/latest/i-tutorials/unit-3.2-navierstokes/navierstokes.html)

Benchmark: Flow around a cylinder

Contents

23.3. Benchmark: Flow around a cylinder#

23.3.1. Exercises#