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.