11.1. Conforming Stokes approx.: Taylor-Hood element

In a first step we want to discretize the Stokes equations with the well known Taylor-Hood finite element method. To this aim we choose the spaces

\[ V_h = [H^1_{0,\Gamma\setminus\Gamma_{\text{out}}}(\Omega)]^d \cap [\mathbb{P}^{k}]^d \quad \text{and} \quad Q_h = H^1(\Omega) \cap \mathbb{P}^{k-1}. \]

This is a suitable pair of spaces for the discretization of the Stokes problem, cf. i-tutorials unit 2.6.

As a first example we want to solve a Stokes flow around a NACA2412 airfoil (naca_geometry.py).

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
from naca_geometry import *

naca_geo = OCCGeometry(occ_naca_profile(type = "2412", height=4, angle=4, h=0.05), dim=2)
mesh = Mesh(naca_geo.GenerateMesh(maxh=0.2,  grading=0.9))
mesh.Curve(3)
Draw(mesh);

The occ_naca_profile function was provided by Xaver Mooslechner. Thanks!

For the viscosity we choose \(\nu = 10^{-4}\).

nu = 1e-4

As boundary conditions we consider wall boundary conditions on the surface of the airfoil "wall" and a constant inflow velocity in front of the airfoil, i.e. the boundary "inlet". The remaining boundaries are outflow boundaries.

V = VectorH1(mesh, order = 3, dirichlet = "inlet|wall")
Q = H1(mesh, order = 2)

X = V * Q 

We set up the variational formulation of the Stokes equation on the product space X and consider a zero right hand side.

(u,p), (v,q) = X.TnT()

K = BilinearForm(X, symmetric = True)
K += (nu * InnerProduct(Grad(u), Grad(v)) - div(u) * q - div(v) * p) * dx
K.Assemble()

F = LinearForm(X).Assemble()

To impose the boundary conditions we consider a standard homogenization step via the residual. Let gfu be the GridFunction that represents the solution \(u,p\) on the product space. We impose the constant inflow value \(u_D = (1,0)\), and solve the problem.

gfu = GridFunction(X)
gfu.components[0].Set(CF((1,0)), definedon = mesh.Boundaries("inlet"))

res = gfu.vec.CreateVector()

res.data = F.vec - K.mat * gfu.vec
gfu.vec.data += K.mat.Inverse(X.FreeDofs()) * res

Draw(gfu.components[1], mesh, "p");
Draw(gfu.components[0], mesh, "u");