12. Instationary Navier-Stokes

We consider the approximation of the full non-linear Navier-Stokes equations. As a first step we use an \(H(\operatorname{div})\) conforming approach and consider an IMEX operator splitting, see NGSolve docu - Navier Stokes.

As example we want to solve a flow through a Tesla valve (tesla.py), see US patent 1329559 “Valvular Conduit”.

The original drawing by N. Tesla in his patent:

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from tesla import GetValve

ddim = 2

# N controls the number of valves (one valve consisting of an upper and lower loop)
valve = GetValve(N = 3, dim = ddim, R = 0.5, alpha = 25, Winlet = 1, beta = 180, L1 = 6.4, L2 = 7, Linlet = 5, closevalve = True)
mesh = Mesh(OCCGeometry(valve, dim=ddim).GenerateMesh(maxh = 0.5))
mesh.Curve(3)

Draw(mesh);

We consider wall boundary conditions everywhere, and as force and viscosity we choose

\[\begin{split} (f,v) = r \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \qquad \text{and} \qquad \nu = 0.005, \end{split}\]

where the factor \(r \in \{-1,1\}\) will determine the flow direction (left to right, or right to left). Our goal is to evaluate the mass flux through the inlet boundary

\[ flux(u_h) = \int_{\Gamma_{in}} u_h \cdot n, \]

over a period of \(T_{end} = 30\).

# viscosity
nu = 0.005
# time step
tau = 0.001
# end time
tend = 30
# approximation order of the velocity
order = 2

For this example we want to make use of the hodivfree flag for the HDiv space. In NGSolve the hierarchical basis functions are split into four groups

• The first group corresponds to the lowest order Raviart-Thomas \(\mathcal{RT}_0\) space whose basis functions are associated to facets. Further note that the divergence of these functions is element wise constant, \(\operatorname{div}{\mathcal{RT}_0} = \mathbb{P}^0 = Q_h^0\), which corresponds to the lowest order pressure space.

• The second group are the remaining facet basis functions which are all divergence free.

• The third group represents n-bubbles (functions with a vanishing normal trace) with a non-zero divergence. More precisely, they correspond to the remaining high order part of the pressures, \(Q_h^0/Q_h^0 = \mathbb{P}^{k-1}/\mathbb{P}^0\).

• The last group are also n-bubbles, but are all divergence free.

Since we want our velocity solution \(u_h\) to be divergence-free, we already know beforehand that we can exclude the third group of basis functions, hodivfree=True. Although this motivates to remove also the lowest order Raviart-Thomas functions (producing a divergence) this can not be done since they are non-local. Accordingly, we still need the lowest-order pressure as a Lagrange multiplier in our system.

Note

More information on the basis function in NGSolve can be found in the PhD thesis of Sabine Zaglmayr “High order finite element methods for electromagnetic field computation”. For the \(H(div)\)-space we also refer to “P.L Lederer, C. Lehrenfeld, J. Schöberl, Hybrid discontinuous Galerkin methods with relaxed H(div)-conformity for incompressible flows, https://doi.org/10.1137/17M1138078”

VT =   HDiv(mesh,order=order, dirichlet="wall", hodivfree=True)
Vhat = TangentialFacetFESpace(mesh,order=order, dirichlet="wall")

Q = L2(mesh,order=0)

X = VT * Q * Vhat

uh, vh = X.TnT()
u, p, uhat = uh
v, q, vhat = vh

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

def tang(u): # tangential part of a vector field
    return u - (u*n)*n
def tjump(u): # hybrid DG jump of the tangential part
    return tang(u[0]-u[2])

We want to use an IMEX splitting method (see NGSolve docu - Navier Stokes). Hence, let \((u,\hat u, p)^n\) be the coefficient vector of the discrete solution \(u_h^n = u_h(t_n), \hat u_h^n = \hat u_h(t_n)\) and \(p_h^n = p_h(t_n)\) for an aequidistant mesh \(0 = t_0 < \ldots < t_N = T_{end}\). We want to solve the system

\[\begin{split} \underbrace{\left( \begin{pmatrix} M & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} + \tau K\right)}_{= M^*} \begin{pmatrix} u \\ \hat u \\ p \end{pmatrix}^{n+1} = \begin{pmatrix} u \\ \hat u \\ p \end{pmatrix}^{n} + \tau \left(f - K \begin{pmatrix} u \\ \hat u \\ p \end{pmatrix}^{n+1} - C(u^n)u^n \right), \end{split}\]

where \(M\) is the mass matrix on VT and \(C(u^n)u^n\) is an explicit evaluation of the non-linear convection term (where we integrated by parts)

\[ [C(u^n)u^n]_i = \sum_T -\int_T u_h^n \otimes u_h^n : \nabla \phi_i + \int_{\partial T} u_h^n \cdot n (u_h^n)^{up} \cdot \phi_i, \]

where \(\phi_i\) is a basis function of the space VT, and \((u_h^n)^{up}\) is the upwinded velocity.

#small regularization so that we use the sparse cholesky solver
stokes = (nu*InnerProduct(Grad(u), Grad(v))- div(u)*q - div(v)*p - nu * 1e-10*p*q)*dx
deb = dx(element_boundary = True)

stokes += nu * Grad(u) * n * tjump(vh) * deb
stokes += nu * Grad(v) * n * tjump(uh) * deb
stokes += nu * 10 * order*order/h * InnerProduct(tjump(uh),tjump(vh)) * deb

K = BilinearForm(stokes, symmetric=True).Assemble()
mstar = BilinearForm(u*v*dx + tau*stokes, symmetric=True, condense = True).Assemble()

Sinv = mstar.mat.Inverse(X.FreeDofs(True), inverse = "sparsecholesky")
inv = ((IdentityMatrix(mstar.mat.height) + mstar.harmonic_extension) @ (Sinv + mstar.inner_solve) @ (IdentityMatrix(mstar.mat.height) + mstar.harmonic_extension_trans))        

The right hand side is only assembled once and will me multiplied with the factor \(r\) in the solving routine.

f = LinearForm(CF((1,0)) * v * dx(definedon = mesh.Materials("valve"))).Assemble()   

For the evaluation of the convection we make use of a geometry-free implementation (see also NGSolve docu) that, on top, is evaluated on the graphics card. This is achieved via the additional flag nonlinear_matrix_free_bdb=True when defining the BilinearForm. In the DG-upwind definition of the convection we need to have access to the neighboring valuve u.Other(), which is not feasible on the graphic card. To this end we consider a reformulation of the boundary term (we omit the superscript \(n\) for the time dependence here)

\[ \int_{\partial T} u_h \cdot n \, u_h^{up} \cdot \phi_i = \int_{\partial T_{out}} u_h \cdot n \, u_h \cdot \phi_i \int_{\partial T_{in}} u_h \cdot n \, u_h^{other} \cdot \phi_i, \]

in the case of an inflow boundary \(T_{in}\). Since \(u_h\) is normal continuous we have

\[ u_h^{other} = (u_h \cdot n) n + (u_h^{other})_t. \]

Now assume that you have stored the mean value of the tangential projection of \(u_h\) in an auxiliary function \(\tilde u\) defined on Vhat, thus \( (\tilde u_h)_t = \frac{1}{2}(u_h + u_h^{other})_t\), then we see

\[ (u_h^{other})_t = (u_h + u_h^{other})_t - (u_h)_t = 2 (\tilde u_h) - (u_h)_t. \]

Note, that in the later equation we do not need to have access to the values of \(u_h\) of the neighboring element. Now with \((u_h)_t = u_n - (u_h \cdot n) n\), thus \((u_h^{other})_t = 2 (\tilde u_h) - u_h + (u_h \cdot n) n\), we get

\[\begin{align*} \int_{\partial T_{in}} (u_h \cdot n) \, u_h^{other} \cdot \phi_i &= \int_{\partial T_{in}} u_h \cdot n \Big( (u_h \cdot n) n + (u_h^{other})_t\Big) \cdot \phi_i \\ &= \int_{\partial T_{in}} (u_h \cdot n ) (u_h \cdot n) (\phi_i \cdot n) + (u_h \cdot n) (u_h^{other})_t \cdot \phi_i\\ &= \int_{\partial T_{in}} (u_h \cdot n) (u_h \cdot n) (\phi_i \cdot n) + (u_h \cdot n) \Big( 2 (\tilde u_h) - u_h + (u_h \cdot n) n\Big)\cdot \phi_i\\ &= \int_{\partial T_{in}} 2 (u_h \cdot n) (u_h \cdot n) (\phi_i \cdot n) + (u_h \cdot n) \Big( 2 (\tilde u_h) - u_h\Big)\cdot \phi_i, \end{align*}\]

The conv operator will be applied to a GridFunction where the proper values of \(u_h\) and \(\tilde u_h\) have to be stored. Note the line with conv.Assemble() where the matrices for the geometry free implementations are built.

conv = BilinearForm(X, nonlinear_matrix_free_bdb=True)
conv += -InnerProduct(grad(v)*u, u) * dx(bonus_intorder=order)

un = u.Operator("normalcomponent")
uhat_t = uhat.Operator("tangentialcomponent")
vn = v.Operator("normalcomponent")
conv += IfPos(un, un * u * v, 2 * un * un * vn + un * ((2*uhat_t-(u)) * v) ) * dx(element_boundary = True, bonus_intorder=order)
conv.Assemble()
conv_operator = conv.mat

In order to calculate the mean value \(\tilde u_h\) we set up the \(L^2\)-projection problem (for a given \(u_h\))

\[ \sum_T \int_{\partial T} (\tilde u_h)_t (\hat v_h)_t = \sum_T \int_{\partial T} (u_h)_t (\hat v_h)_t \quad \forall \hat v_h \in \hat V_h. \]

Since the tangential component of functions in Vhat are single valued on facets (in contrast to \(u_h\)) we have on a common facet

\[ 2 (\tilde u_h)_t (\hat v_h)_t = \Big((u_h)_t + (u_h^{other})_t \Big)(\hat v_h)_t. \]

mfacet = BilinearForm(X)
mfacet += tang(uhat) * tang(vhat) * dx(element_boundary=True)
mfacet.Assemble()

fdofs = BitArray(X.ndof)
fdofs.Clear()
fdofs[VT.ndof + Q.ndof:X.ndof] = 1

mfinv = mfacet.mat.Inverse(fdofs, inverse = "sparsecholesky")
average_rhs = BilinearForm(X, nonassemble = True)
average_rhs += tang(u) * tang(vhat) * dx(element_boundary = True)

t = 0
cnt = 0

# list to store the mass flux
massflux = []
massflux_reverse = []
time = []

# solution functions
gfu = GridFunction(X)
gfu_reverse = GridFunction(X)
gfu.vec[:]= 0
gfu_reverse.vec[:]= 0

average_vec = gfu.vec.CreateVector()
gfu_average = GridFunction(X)

# helping vectors
res = gfu.vec.CreateVector()
convvec = gfu.vec.CreateVector()
convvec[:] = 0.0
convvec_reverse = gfu.vec.CreateVector()
convvec_reverse[:] = 0.0


gfut = GridFunction(gfu.components[0].space, multidim=0)
gfut_reverse = GridFunction(gfu.components[0].space, multidim=0)

gfut.AddMultiDimComponent(gfu.vec)
gfut_reverse.AddMultiDimComponent(gfu_reverse.vec)

with TaskManager():
    while t < tend:   

        # solve for the forward operator
        average_rhs.Apply(gfu.vec, average_vec)
        gfu_average.vec.data = mfinv * average_vec
        gfu_average.components[0].vec.data = gfu.components[0].vec
        convvec.data = conv_operator * gfu_average.vec

        res.data = f.vec - convvec - K.mat * gfu.vec
        gfu.vec.data += tau * inv * res    

        # solve for the backward operator
        average_rhs.Apply(gfu_reverse.vec, average_vec)
        gfu_average.vec.data = mfinv * average_vec
        gfu_average.components[0].vec.data = gfu_reverse.components[0].vec
        convvec_reverse.data = conv_operator * gfu_average.vec

        res.data = -f.vec - convvec_reverse - K.mat * gfu_reverse.vec
        gfu_reverse.vec.data += tau * inv * res    
        
        t = t + tau
        time.append(t)
    
        flux = Integrate(gfu.components[0]*n, mesh, definedon = mesh.Boundaries("inlet"))
        flux_reverse = Integrate(gfu_reverse.components[0]*n, mesh, definedon = mesh.Boundaries("inlet"))
        massflux.append(abs(flux))
        massflux_reverse.append(abs(flux_reverse))
        cnt += 1

        if cnt % int(1.5/(tau)) == 0:
            gfut.AddMultiDimComponent(gfu.components[0].vec)
            gfut_reverse.AddMultiDimComponent(gfu_reverse.components[0].vec)

In the following plots we can see the solutions at \(t = T_{end}\). We observe that for the flow from left to right the solution is much more turbulent and the valve is “activated” since the fluid stops it self. On the other hand, for the flow from right to left, the solution is nearly laminar and much faster. This can also be observed in the the maxs flux ploted over time below.

Draw (gfut, mesh, interpolate_multidim=True, animate=True,
      min=0, max=4, autoscale=False, vectors = {"grid_size": 100});
Draw (gfut_reverse, mesh, interpolate_multidim=True, animate=True,
      min=0, max=4, autoscale=False, vectors = {"grid_size": 100});

import matplotlib.pyplot as plt

plt.plot(time, massflux, "-", color = "orange", label = "left_to_right")
plt.plot(time, massflux_reverse, "-", color = "blue", label = "right_to_left")
plt.legend(loc = "lower right")
plt.show()