This page was generated from wta/navierstokes.ipynb.
Navier Stokes Equations¶
Find velocity u:Ω×[0,T]→Rd and pressure p:Ω×[0,T]→R such that
∂u∂t−νΔu+u∇u+∇p=fdivu=0
[1]:
from netgen.geom2d import SplineGeometry
from ngsolve import *
from ngsolve.webgui import Draw
Schäfer-Turek benchmark geometry:
[2]:
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)
[2]:
Higher order Taylor-Hood element pairing:
[3]:
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:
[4]:
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")
[4]:
solve Stokes problem for initial conditions:
[5]:
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:
un+1−unτ−νΔun+1+∇pn+1=f−un∇un
and
divun+1=0
[6]:
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 ∫u∇uv
[7]:
conv = BilinearForm(X, nonassemble = True)
conv += (Grad(u) * u) * v * dx
implicit Euler/explicit Euler splitting method:
[8]:
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)
[9]:
Draw (gfut, mesh, interpolate_multidim=True, animate=True)
[9]: