Nitsche method for Stokes flow.

Question

Hello everybody. I am trying to simulate fluid flow around an object using the incompressible Stokes flow equations. I want to define a free-slip boundary condition around an "obstacle" with a curved surface in the middle of the domain, and to do that I'm trying to use the Nitsche method. Here is a paper describing the Nitsche method (the relevant material is on pages 5-6). And here somebody actually implemented the Nitsche method in FEniCS, but I am unable to run it with my current version (2016.2). I tried to structure my own program using this person's implementation.

So here's the problem. I am trying to benchmark my program using an analytical solution, but my simulation results are nowhere close. In fact, the simulation results seem to better match the case where a no-slip boundary condition is applied around the "obstacle" in the middle of the domain. I am unable to figure out where the discrepancy is coming from. Does anybody have experience with fluid mechanics who can help me out?

I'll post my code below, but here are some relevant details:

• The domain under consideration is pretty simple. It's just a 10 x 10 square with a circular hole in the center, and the radius of the hole is 1. I want to apply a free-slip boundary condition around this hole so that the velocity normal to the surface is zero.

• I'm using Taylor-Hood elements, so the problem is defined in a mixed function space with quadratic interpolation for velocity and linear interpolation for pressure.

• For Dirichlet boundary conditions, I have imposed free slip on the top and bottom edges of the domain. On the left and right edges, I have set the velocity as Constant((1,0)).

And here is the full code:

# import Python modules
from __future__ import division
from fenics import *

# suppress FEniCS output to screen
set_log_active(False)

# load mesh
mesh = Mesh()
hdf5 = HDF5File(mpi_comm_world(),"mesh/cylinder.h5","r")
hdf5.read(mesh,"mesh",False)
hdf5.close()

# define function spaces
P1 = FiniteElement("CG", mesh.ufl_cell(), 1) # first order scalar finite element
P2 = VectorElement("CG", mesh.ufl_cell(), 2) # second order vector finite element
V = FunctionSpace(mesh, MixedElement([P2,P1])) # mixed finite element

# classes for domain edge boundaries
class left_edge(SubDomain):
    def inside(self, x, on_boundary): return near(x[0],0) and on_boundary
class right_edge(SubDomain):
    def inside(self, x, on_boundary): return near(x[0],10) and on_boundary
class bottom_edge(SubDomain):
    def inside(self, x, on_boundary): return near(x[1],0) and on_boundary
class top_edge(SubDomain):
    def inside(self, x, on_boundary): return near(x[1],10) and on_boundary
class hole_edge(SubDomain):
    def inside(self, x, on_boundary): return ((x[0]-5)**2+(x[1]-5)**2 <= 1.1) and on_boundary

left = left_edge()
right = right_edge()
bottom = bottom_edge()
top = top_edge()
hole = hole_edge()

# mark boundaries as "ds"
boundaries = FacetFunction("size_t",mesh,0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
bottom.mark(boundaries, 3)
top.mark(boundaries, 4)
hole.mark(boundaries, 5)
ds = Measure("ds",subdomain_data=boundaries)

# Dirichlet BC's
free_slip_top = DirichletBC(V.sub(0).sub(1), Constant(0), top)
free_slip_bottom = DirichletBC(V.sub(0).sub(1), Constant(0), bottom)
flow_left = DirichletBC(V.sub(0), Constant((1,0)), left)
flow_right = DirichletBC(V.sub(0), Constant((1,0)), right)
BC = [free_slip_top, free_slip_bottom, flow_left, flow_right]

# symmetric gradient operator
def D(x):
    return sym(nabla_grad(x))

# linear viscosity
eta = Constant(1)

# stress
def sig(x,y):
    return 2*eta*D(x) - y*Identity(2)

# trial and test functions for velocity and pressure
u, p = TrialFunctions(V) # trial functions in (P2, P1) space
v, q = TestFunctions(V) # test functions in (P2, P1) space

# loading terms
f = Constant((0, 0))
g = Constant(0)

# penalty terms
h = CellSize(mesh)
beta = Constant(10)

# variational form
n = FacetNormal(mesh) # normal vector
t = dot(sig(u,p),n) # traction
s = dot(sig(v,q),n) # traction

# stokes flow main terms
a = (2*eta*inner(D(v),D(u)) - p*div(v) + q*div(u))*dx
L = inner(f,v)*dx

# nitsche terms
a += (
    - inner(dot(t,n),dot(v,n))
    - inner(dot(s,n),dot(u,n))
    + beta/h*inner(dot(u,n),dot(v,n))
    )*ds(5)
L += -inner(g,inner(s,n))*ds(5)

# solve the variational form
w = Function(V)
solve(a==L, w, BC)
u, p = w.split(deepcopy=True)

# store data
hdf5 = HDF5File(mpi_comm_world(),"output/data.h5","w")
hdf5.write(mesh,"mesh")
hdf5.write(u,"velocity")
hdf5.write(p,"pressure")
hdf5.close()

Thank you in advance.

Comments

Did you check whether the facets are all marked right? If the 1.1 tolerance is too high for your particular mesh you may find a whole layer of element boundaries and mess up your weak constraints. EDIT: sorry, I overlooked that you check for on_boundary. My comment is probably not helpful then.

Answer 1

First you should note that your system is singular.

Since you are prescribing the normal component of the velocity field on all boundaries of the mesh, then the pressure is defined up to a constant.

Krylov methods may be able to handle this singular system, provided the right hand side is orthogonal to the null space of the operator).
On the other hand, direct solver (like the one I'm expecting dolfin is using as default) may return the wrong solution.

Possible walkaround are:

1. Prescribe a natural (i.e. do nothing) boundary condition on the outflow boundary.

2. Add a lagrangian multiplier lmbda in FunctionSpace(mesh, 'R', 0) and modify the weak form to impose a zero-mean pressure.

3. (Mathematically hugly, but ok in practice) Define a pointwise boundary condition for the pressure, so that you fix the value of the pressure at an arbitrary degree of freedom.

Let me know if this help.

Umbe

Comments

Thanks for the response Umbe. Unfortunately, workarounds 1 and 3 did not work for me. The velocity solution did not change when I fixed the pressure at a point. I am also not certain how to impose the zero-mean constraint in the variational form. Here is the method I have tried:

Function spaces:

P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
P2 = VectorElement("CG", mesh.ufl_cell(), 2)
R0 = FiniteElement("R", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, MixedElement([P2,P1,R0]))

Trial and test functions:

u, p, rho = TrialFunctions(V) # trial functions in (P2, P1) space
v, q, lam = TestFunctions(V) # test functions in (P2, P1) space

Variational form:

# main terms
a = (2*eta*inner(D(v),D(u)) - p*div(v) + q*div(u))*dx
L = inner(f,v)*dx

# lagrange terms
a += (p*lam + rho*q)*dx

# nitsche terms
a += (
    - inner(dot(t,n),dot(v,n))
    - inner(dot(s,n),dot(u,n))
    + beta/h*inner(dot(u,n),dot(v,n))
        )*ds(5)
L += -inner(g,inner(s,n))*ds(5)

Is this correct?
Thanks again for your help.

---

Jimmy, "workaround" 1 should definitely work. Just set a zero Neumann BC for
the traction (dot(n, sigma) = 0) at, for instance, the right boundary. To
do this, simply delete the Dirichlet BC defined on that boundary. This BC fixes the
pressure (integral of the pressure over that boundary is zero) and gives you a physically meaningful solution.
I just ran your code and the velocity field looks rather good? (It won't change that much when you
use the Neumann BC). There is definitely "slip" along the cylinder.

---

I did get a velocity solution when using workaround 1, however it does not seem to match the analytical solution that I posted for free-slip around the cylinder. I do not think that the free-slip boundary condition is being correctly applied to the cylinder in the middle of the domain, and that is the reason I said the workaround did not work.

Answer 2

The solution given in your link is for potential flow (incompressible, inviscid and irrotational flow)which can approximate a very high Reynolds number (Re) flow but can't take into account skin friction and flow separation. The Stokes flow approximation is for very low Re number flows.

What I would suggest to solve in fenics is the potential flow equation (section 3.6.2):

http://web.mit.edu/2.20/www/lectures/2016/Lecture_9_2016.pdf

This equation has the same analytic solutions that are shown in your link

What you have to do is solve a Poisson equation:

https://fenicsproject.org/olddocs/dolfin/1.5.0/python/demo/documented/poisson/python/documentation.html

And take the gradient of the solution to get the velocity field (pg 5):

http://people.cs.uchicago.edu/~ridg/FEniCScourses/chiphoecrse.pdf

Comments

Thanks for the response alexmm. Ultimately I need to solve fluid flow problems using the Stokes equations. Eventually I plan on solving more complicated problems than the one that I posted (i.e., with different domain geometries and a nonlinear fluid viscosity), and that problem I posted was more or less a benchmark. The main goal right now is to implement the Nitsche method properly. Maybe I need to find a better benchmark problem...

---

In that case you can benchmark using the Stokes solution for flow around a cylinder:

http://www.lmm.jussieu.fr/~lagree/COURS/M2MHP/petitRe.pdf

You can get a quantitative measure of the accuracy of your solution by computing the drag on the cylinder and checking that against the analytic solution for the drag

---

It seems this paper covers the solution for the no-slip boundary condition, however I am trying to model free-slip. Along the boundary of the circular obstacle, I want to enforce $\boldsymbol{u} \cdot \boldsymbol{\hat{n}}=0$ rather than $\boldsymbol{u} = \boldsymbol{0}$.

---

May I ask what the application is? Usually Stokes problems use a no-slip boundary condition

---

The application is glacier dynamics - I'm using the Stokes equations to model the creeping flow of glaciers. Up till now, I've been modeling glacier flow over a flat, frictionless ground line. So basically, the bottom of the domain is treated as a free-slip surface (similar to the "rollers" boundary condition in an elasticity problem). But in the near future I want to model free slip over arbitrarily shaped (e.g., curved or slanted) ground lines.