Slip condition on flat boundary, Stokes with constant viscosity

variational formulation
1

Hello,
I would like to impose a slip condition for my problem, and I am aware that a good option is Nitsche’s method, together with dolfin_dg. As I am just starting out, I wanted to run the simplest possible case of my problem, which is Stokes equations with a constant viscosity. My boundary conditions are a bit unusual though. I have a square domain, with outflow at the bottom, a moving wall on the right and a symmetry condition at the centre. However, at the top, my domain is split into two halves. I have an inlet velocity on the left half, and I would like a slip condition on the right half. The slip condition should be that the normal velocity is zero and the tangential stress is zero. I believe the latter is satisfied automatically, and I used W.sub(0).sub(1) to impose zero normal velocity, as with a flat boundary it would correspond to the second component. However, it does not converge. Is there anything I am missing? I am worried that there is something fundamental that I don’t understand which in turn will make it much harder to implement Nitsche’s for my ‘real’ problem. Here is a minimal example of the code

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry

mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)

# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))

# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))

# Define the viscosity and bcs

mu = Constant(1.0)
u_in = Constant(-2.0) # inlet velocity
u_c = Constant(-1.0) # wall velocity, and outlet velocity

inflow = 'near(x[1], 1.0) && x[0]<=0.5' # left half of top boundary
wall = 'near(x[0], 1.0)' # right boundary
centre = 'near(x[0], 0.0)' # left boundary
outflow = 'near(x[1], 0.0)' # bottom boundary
weird = 'near(x[1], 1.0) && x[0]>=0.5' # right half of top boundary
bcu_inflow = DirichletBC(W.sub(0), (0.0, u_in), inflow)
bcu_wall = DirichletBC(W.sub(0), (0.0, u_c), wall)
bcu_outflow = DirichletBC(W.sub(0), (0.0, u_c), outflow)
bcu_symmetry = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre)
bcu_slip = DirichletBC(W.sub(0).sub(1), Constant(0.0), weird) # slip condition that should be right but is somehow wrong
bcs = [bcu_inflow, bcu_wall, bcu_outflow, bcu_symmetry, bcu_slip]

x = SpatialCoordinate(mesh)

# Define stress tensor
def epsilon(v):
   return sym(as_tensor([[v[0].dx(0), v[0].dx(1), 0],
                         [v[1].dx(0), v[1].dx(1), 0],
                         [0, 0, 0]]))

# stress tensor
def sigma(v, p):
    return 2*mu*epsilon(v)-Id(p)

def Id(p):
    return as_tensor([[p, 0, 0],
                      [0, p, 0],
                     [0, 0, p]])

def cond(v):
    return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
                          [v[1].dx(0), v[1].dx(1)]]))

def sigmabc(v, p):
    return 2*mu*cond(v) - p*Identity(2)


# Define the variational form
f = Constant((0, -1)) # forcing term ( I have also tried with this being zero, still no convergence)

colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0)  # default to zero
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 5)
CompiledSubDomain("near(x[1], 1.0) && x[0]<=0.5").mark(colors, 1) # top left
CompiledSubDomain("near(x[1], 1.0) && x[0]>=0.5").mark(colors, 2) # top right
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 3)  # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 4)  # outflow

# Create the measure
ds = Measure("ds", subdomain_data=colors)

a1 = (inner(sigma(u, p), epsilon(v))) * dx
a2 = (- div(u) * q - dot(f, v)) * dx

b1 = - dot(dot(sigmabc(u, p), v), n) * ds(1)
b2 = -dot(dot(sigmabc(u, p), v), n) * ds(2)
b3 = - dot(dot(sigmabc(u, p), v), n) * ds(3)
b4 = - dot(dot(sigmabc(u, p), v), n) * ds(4)

F = a1 + a2 + b1 + b2 + b3 + b4

# Solve problem
solve(F == 0, w, bcs)

# Plot solutions
(u, p) = w.split()

File("Results/velocityIsothermalFreeSlip.pvd") << u
c = plot(u, title='Velocity Isothermal')
plt.colorbar(c)
plt.show()

I really hope it is something simple I have overlooked! I appreciate your time.

2

What’s the purpose of b1 through b4? I don’t see how they can appear in the variational formulation.

After removing these it looks like you’ve over prescribed the BCs and therefore your problem (as written) is ilposed.

Intuitively you’re forcing mass into an incompressible system from which it cannot escape.

A simple demonstration which will yield a result:
Your “outflow” BCs are not actually outflow BCs. E.g. just removing your “bcu_outflow” from the bcs list will enforce a natural outflow condition and mass may escape the system.

1 Like
3

Thank you! That’s right, it now matches the results I get with Nitsche’s method. I thought that if I don’t impose b1 to b4 I am setting zero stress on all the boundaries, but I don’t want to impose zero stress on all the boundaries … Did I misunderstand?