Navier Stokes Equations

5. Navier Stokes Equations#

Find velocity $u : \Omega \times [0,T] \rightarrow {\mathbb R}^d$ and pressure $p : \Omega \times [0,T] \rightarrow {\mathbb 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

from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
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"
DrawGeo (shape)
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.07)).Curve(3)
Draw (mesh);

H(div)-conforming HDG for Stokes

order=4
VT = HDiv(mesh, order=order, dirichlet="wall|cyl|inlet")
VF = TangentialFacetFESpace(mesh, order=order, dirichlet="wall|cyl|inlet")
Q = L2(mesh, order=order-1)
X = VT*VF*Q

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

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size
dS = dx(element_boundary=True)

nu = 1e-3

def tang(vec):
    return vec - (vec*n)*n

stokes = nu*InnerProduct(Grad(u), Grad(v)) * dx \
    + nu*InnerProduct(Grad(u)*n, tang(vhat-v)) * dS \
    + nu*InnerProduct(Grad(v)*n, tang(uhat-u)) * dS \
    + nu*4*order*order/h * InnerProduct(tang(vhat-v), tang(uhat-u)) * dS \
    + div(u)*q*dx + div(v)*p*dx -1e-10/nu*p*q*dx

a = BilinearForm (stokes).Assemble()
f = LinearForm(X).Assemble()

solve Stokes problem for initial condition:

gfu = GridFunction(X)

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(), inverse="sparsecholesky")

res = f.vec.CreateVector()
res.data = 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(X)
mstar += u*v*dx+tau*stokes
mstar.Assemble()
inv = mstar.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")

the non-linear convective term $\int u \nabla u v$ with up-wind:

conv = BilinearForm(X, nonassemble = True)

conv += -InnerProduct(Grad(v)*u, u)*dx

un = u*n
upwind = IfPos(un, un*u, un*u.Other(bnd=uin))

conv += InnerProduct (upwind, v) * dx(element_boundary=True)

time-loop with IMEX time-stepping

t = 0; i = 0
tend = 10
scene = Draw (gfu.components[0], mesh, min=0, max=2, autoscale=False, order=3)

with TaskManager():
    while t < tend:
        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()