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)

 Generate Mesh from spline geometry
 Boundary mesh done, np = 86
 CalcLocalH: 86 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 load internal triangle rules
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters
 Curve elements, order = 3
 Update clusters
 Curving edges
 Curving faces

BaseWebGuiScene

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)

assemble VOL element 482/482

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")

setvalues element 86/86

BaseWebGuiScene

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)

call pardiso ... done

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")

assemble VOL element 482/482

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)

BaseWebGuiScene