This page was generated from wta/navierstokes.ipynb.

Navier Stokes Equations

Find velocity \(u : \Omega \times [0,T] \rightarrow R^d\) and pressure \(p : \Omega \times [0,T] \rightarrow R\) such that

\[\begin{split}\begin{array}{ccccl} \frac{\partial u}{\partial t} - \nu \Delta u + u \nabla u & + & \nabla p & = & f \\ \operatorname{div} u & & & = & 0 \end{array}\end{split}\]

from netgen.geom2d import SplineGeometry
from ngsolve import *
from ngsolve.webgui import Draw

Schäfer-Turek benchmark geometry:

geo = SplineGeometry()
geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc="cyl", maxh=0.02)
mesh = Mesh(geo.GenerateMesh(maxh=0.07))

mesh.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 = FESpace([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

a = BilinearForm(X)
a += stokes*dx
a.Assemble()

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

# gridfunction for the solution
gfu = GridFunction(X)

parabolic inflow at inlet:

uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))
# Draw (Norm(gfu.components[0]), mesh, "velocity", sd=3)
Draw (gfu.components[0], mesh, "vel")

solve Stokes problem for initial conditions:

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

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

scene = Draw (gfu.components[0], mesh)

implicit/explicit time-stepping:

\[\frac{u_{n+1}-u_n}{\tau} - \nu \Delta u_{n+1} + \nabla p_{n+1} = f - u_n \nabla u_n\]

and

\[\operatorname{div} u_{n+1} = 0\]

tau = 0.001 # timestep parameter

mstar = BilinearForm(X)
mstar += (u*v+tau*stokes)*dx
mstar.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
tend = 1
i = 0
gfut = GridFunction(V, multidim=0)
vel = gfu.components[0]
with TaskManager():
    while t < tend:
        # print ("t=", t, end="\r")

        conv.Apply (gfu.vec, res)
        res.data += a.mat*gfu.vec
        gfu.vec.data -= tau * inv * res

        t = t + tau
        i = i + 1
        if i%10 == 0:
            scene.Redraw()
            gfut.AddMultiDimComponent(vel.vec)

Draw (gfut, mesh, interpolate_multidim=True, animate=True)