We are solving the unsteady Navier Stokes equations
Find $(u,p):[0,T] \to (H_{0,D}^1)^d \times L^2$, s.t.
\begin{align} \int_{\Omega} \partial_t u \cdot v + \int_{\Omega} \nu \nabla u \nabla v + u \cdot \nabla u v - \int_{\Omega} \operatorname{div}(v) p &= \int f v && \forall v \in (H_{0,D}^1)^d, \\ - \int_{\Omega} \operatorname{div}(u) q &= 0 && \forall q \in L^2, \\ \quad u(t=0) & = u_0 \end{align}import netgen.gui
%gui tk
import tkinter
from math import pi
from ngsolve import *
from netgen.geom2d import SplineGeometry
We again consider the benchmark setup of from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry is a 2D channel with a circular obstacle which is positioned (only slightly) off the center of the channel. The geometry:
The left part (red) is the inflow boundary,
The right part (green) is the outflow boundary.
The viscosity is set to $\nu = 10^{-3}$.
from netgen.geom2d import SplineGeometry
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")
mesh = Mesh( geo.GenerateMesh(maxh=0.08))
mesh.Curve(3)
# viscosity
nu = 0.001
We use a Taylor-Hood discretization of degree $k$. The finite element space for the velocity components $u_x$ and $u_y$ and the pressure are:
\begin{align} u_x,u_y \in V &= \{ v \in H^1(\Omega) | v|_T \in \mathcal{P}^k(T) \} \\ p \in Q &= \{ q \in H^1(\Omega) | q|_T \in \mathcal{P}^{k-1}(T) \} \end{align}k = 3
V = H1(mesh,order=k, dirichlet="wall|cyl|inlet")
Q = H1(mesh,order=k-1)
X = FESpace([V,V,Q])
As initial data we use the solution to the Stokes problem:
Find $u=(u_x,u_y) \in V \times V$, $p \in Q$ so that
\begin{align} \int_{\Omega} \nu \nabla u \nabla v - \int_{\Omega} \operatorname{div}(v) p &= \int f v && \forall v \in V \times V, \\ - \int_{\Omega} \operatorname{div}(u) q &= 0 && \forall q \in Q, \end{align}with boundary conditions:
u = GridFunction(X)
velocity = CoefficientFunction (u.components[0:2])
Draw(velocity,mesh,"u")
Draw(u.components[2],mesh,"p")
# parabolic inflow at bc=1:
uin = 1.5*4*y*(0.41-y)/(0.41*0.41)
u.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))
Redraw()
# a sketch of the inflow profile:
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
s = np.arange(0.0, 0.42, 0.01)
bvs = 6*s*(0.41-s)/(0.41)**2
plt.plot(s, bvs)
plt.xlabel('y')
plt.ylabel('$u_x(0,y)$')
plt.title('inflow profile')
plt.grid(True)
plt.show()
To solve the Stokes (and later the Navier-Stokes) problem, we introduce the bilinear form:
$$ a((u,p),(v,q)) := \int_{\Omega} \nu \nabla u \nabla v - \operatorname{div}(v) p - \operatorname{div}(u) q $$ux,uy,p = X.TrialFunction()
vx,vy,q = X.TestFunction()
div_u = grad(ux)[0]+grad(uy)[1]
div_v = grad(vx)[0]+grad(vy)[1]
stokes = nu*grad(ux)*grad(vx)+nu*grad(uy)*grad(vy) \
-div_u*q-div_v*p
a = BilinearForm(X)
a += SymbolicBFI(stokes)
a.Assemble()
f = LinearForm(X)
f.Assemble()
inv_stokes = a.mat.Inverse(X.FreeDofs())
res = f.vec.CreateVector()
res.data = f.vec - a.mat*u.vec
u.vec.data += inv_stokes * res
Redraw()
For the time integration we consider an semi-implicit Euler method (IMEX) where the convection is treated only explicitly and the remaining part implicitly:
Find $(u^{n+1},p^{n+1}) \in X = V \times V \times Q$, s.t. for all $(v,q) \in X = V \times V \times Q$ \begin{equation} \underbrace{m(u^{n+1},v) + \Delta t ~\cdot~a((u^{n+1},p^{n+1}),(v,q))}_{ \to M^\ast} ~=~ m(u^{n},v) - \Delta t c(u^{n}; u^{n},v) \end{equation}
with \begin{equation} m(u,v) = \int u \cdot v \end{equation} and \begin{equation} c(w,u,v) = \int w \cdot \nabla u \cdot v \end{equation}
We prefer the incremental form (as it homogenizes the boundary conditions): \begin{equation} m(\Delta u^{n+1},v) + \Delta t ~\cdot~a((\Delta u^{n+1},p^{n+1}),(v,q)) ~=~ \Delta t (-a((u^{n},p^n),(v,q)) -c(u^{n}; u^{n},v)) \end{equation}
dt = 0.001
# matrix for implicit Euler
mstar = BilinearForm(X)
mstar += SymbolicBFI(ux*vx+uy*vy + dt*stokes)
mstar.Assemble()
inv = mstar.mat.Inverse(X.FreeDofs())
velocity = CoefficientFunction (u.components[0:2])
conv = LinearForm(X)
conv += SymbolicLFI(velocity*grad(u.components[0])*vx)
conv += SymbolicLFI(velocity*grad(u.components[1])*vy)
t = 0
tend = 0
# implicit Euler/explicit Euler splitting method:
tend += 1
while t < tend-0.5*dt:
print ("\rt=", t, end="")
conv.Assemble()
res.data = a.mat * u.vec + conv.vec
u.vec.data -= dt * inv * res
t = t + dt
Redraw()
In the previously used benchmark the stresses on the obstacle are evaluated, the so-called drag and lift coefficients.
The force acting on the obstacle is $$ F_{\circ} = (F_D,F_L) = \int_{\Gamma_{cyl}} \sigma_n = \int_{\Gamma_{cyl}} (-\nu \frac{\partial u}{\partial n} + p I) \cdot n $$
The drag/lift coefficients are $$ c_D = \frac{2 }{\bar{u} L} F_D, \quad c_L = \frac{2 }{\bar{u} L} F_L. $$ where $\bar{u} = 1$ is the mean inflow velocity and $L = 0.41$ is the channel width.
We use the residual of our discretization to compute the forces. On the boundary degrees of freedoms of the disk we overwrote the equation by prescribing the (homogeneous) boundary conditions. The equations related to these dofs describe the force (im)balance at this boundary.
Testing the residual (functional) with the characteristic function on that boundary (in the x- or y-direction) we obtain the integrated stresses (in the x- or y-direction):
We define the functions which are characteristic functions (in terms of boundary evaluations)
drag_x_test = GridFunction(X)
drag_x_test.components[0].Set(CoefficientFunction(-20.0), definedon=mesh.Boundaries("cyl"))
drag_y_test = GridFunction(X)
drag_y_test.components[1].Set(CoefficientFunction(-20.0), definedon=mesh.Boundaries("cyl"))
We will collect drag and lift forces over time and therefore create empty arrays
and reset initial data (and time)
time_vals = []
drag_x_vals = []
drag_y_vals = []
# restoring initial data
res.data = f.vec - a.mat*u.vec
u.vec.data += inv_stokes * res
Redraw()
t = 0
tend = 0
With the same discretization as before.
Remarks:
# implicit Euler/explicit Euler splitting method:
tend += 1
while t < tend-0.5*dt:
print ("\rt=", t, end="")
conv.Assemble()
res.data = a.mat * u.vec + conv.vec
u.vec.data -= dt * inv * res
t = t + dt
Redraw()
time_vals.append( t )
drag_x_vals.append(InnerProduct(res, drag_x_test.vec) )
drag_y_vals.append(InnerProduct(res, drag_y_test.vec) )
#print(drag)
# Plot drag/lift forces over time
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(time_vals, drag_x_vals)
plt.xlabel('time')
plt.ylabel('drag')
plt.title('drag')
plt.grid(True)
plt.show()
plt.plot(time_vals, drag_y_vals)
plt.xlabel('time')
plt.ylabel('lift')
plt.title('lift')
plt.grid(True)
plt.show()