32. 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 ngsolve import *
from ngsolve.webgui import Draw

Schäfer-Turek benchmark geometry:

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:

uin = CoefficientFunction( (1.5*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}{\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

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
gfut = GridFunction(V, multidim=0)
vel = gfu.components[0]
scene = Draw (gfu.components[0], mesh, min=0, max=2, 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()
        if i%50 == 0: 
            gfut.AddMultiDimComponent(vel.vec)
            # print(f"t = {t}", end='\r')

In the multidim - GridFunction gfut we have collected several time-steps, which can be animated now by the visualization:

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