H(div)-conforming Stokes

45. H(div)-conforming Stokes#

[Master thesis Lehrenfeld 2010, Lehrenfeld+Schöberl, 2016]

  • Discretize velocity \(u\) in \(H(\operatorname{div})\), with Raviart-Thomas or BDM elements. This allows to use pressure space \(Q_h = \operatorname{div} V_h\). This the discrete velocity is exactly divergence free.

  • Since RT/BDM elements have only continuous normal components, they cannot be used directly for the viscosity term. So we use HDG for gluing together also the tangential components.

from ngsolve import *
from ngsolve.webgui import Draw

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"

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.07)).Curve(3)
Draw (mesh);

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)
V = VT*VF*Q

u, uhat, p = V.TrialFunction()
v, vhat, q = V.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

# Thesis Christoph Lehrenfeld, page 71
a = BilinearForm (V)
a += nu*InnerProduct(Grad(u), Grad(v)) * dx
a += nu*InnerProduct(Grad(u)*n, tang(vhat-v)) * dS
a += nu*InnerProduct(Grad(v)*n, tang(uhat-u)) * dS
a += nu*4*order**2/h * InnerProduct(tang(vhat-v), tang(uhat-u)) * dS
a += div(u)*q*dx + div(v)*p*dx
a.Assemble();

f = LinearForm(V).Assemble()

invstokes = a.mat.Inverse(V.FreeDofs())

gfu = GridFunction(V)
uin = CoefficientFunction((1.5*4*y*(0.41-y)/(0.41*0.41),0))
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

res = f.vec - a.mat*gfu.vec
gfu.vec.data += invstokes * res

Draw (gfu.components[0], mesh, "vel") # velocity
Draw (gfu.components[2], mesh, "pressure")  # pressure
Draw (div(gfu.components[0]), mesh, "divvel"); # div velocity