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.0) && x<=0.5' # left half of top boundary wall = 'near(x, 1.0)' # right boundary centre = 'near(x, 0.0)' # left boundary outflow = 'near(x, 0.0)' # bottom boundary weird = 'near(x, 1.0) && x>=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.dx(0), v.dx(1), 0], [v.dx(0), v.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.dx(0), v.dx(1)], [v.dx(0), v.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)").mark(colors, 5) CompiledSubDomain("near(x, 1.0) && x<=0.5").mark(colors, 1) # top left CompiledSubDomain("near(x, 1.0) && x>=0.5").mark(colors, 2) # top right CompiledSubDomain("near(x, 1.0)").mark(colors, 3) # wall CompiledSubDomain("near(x, 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.