Hello,
I am solving the Stokes equations on the unit square. Along one of the boundaries, I would like to prescribe the normal component of velocity and tangential component of traction, i.e.,
$$ \vec u\cdot\vec n=g $$
and
$$ \vec n\cdot\sigma(\vec u, p)\cdot\vec t = f. $$
This exact question was asked previously [1] in the context of complex geometries, but I got lost in the rather general implementation that was provided in the solution. [2] asked for independent control of the stress components. [3] asked an ill-posed version of the same problem. I got the first line of code from [4].
[1] https://fenicsproject.org/qa/9926/boundary-conditions-in-normal-and-tangent-directions
[2] https://fenicsproject.org/qa/8335/normal-tangential-stress-boundary-conditions-elasticity
[3] https://fenicsproject.org/qa/8543/tangential-gradient-boundary-conditions
[4] https://fenicsproject.org/qa/5236/setting-boundary-condition-vector-component-dirichletbc
Below is a minimal working example of the environment where I would like to implement this type of boundary condition. I would like to implement the "normal velocity, tangential traction" boundary condition on "bed".
from dolfin import *
mesh = UnitSquareMesh(32,32)
n = FacetNormal(mesh)
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
def left(x, on_boundary): return near(x[0],0.0)
inflow = Expression(("-x[1]*(x[1]-2)", "0.0"), degree=2)
bc1 = DirichletBC(W.sub(0), inflow, left)
class SlidingBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0)
bed= SlidingBoundary()
bc0 = DirichletBC(W.sub(0).sub(1), Constant(0), bed)
bcs = [bc1,bc0]
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
bed.mark(boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
g = Constant((1, 0))
a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx
L = inner(f, v)*dx + inner(g,v)*ds(1)
w = Function(W)
solve(a == L, w, bcs)
u, p = w.split()