Prescribed normal velocity and tangential traction

Question

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()

Comments

HI,

Do you find a way to solve this problem ?
Hoping for your recent update!

Best Regards!

Hamilton

Answer 1

May you please post the code until the expression for L?

Comments

I've edited the original post to include a minimal working code. It makes no attempt to implement the type of boundary condition that I'd like to implement because, well, I don't know how to do it :)

Answer 2

I thought through the variational form for this boundary condition. I would just like to ask for the community's help to translate this form into FENICS.

Performing a standard integration by parts on the momentum balance equations (as in this demo) results in a vector-valued boundary integral. I'll focus only on the boundary y=0, x=[0,1]. The individual components of the boundary integral are,
$$\int_0^1 \phi_1 (\nabla u_1 \cdot n) dx,$$
$$\int_0^1 \phi_2 (\nabla u_2 \cdot n) dx,$$
or, expanding the vector notation,
$$\int_0^1 \phi_1 \left(n_1 \frac{\partial u_1}{\partial x} + n_2 \frac{\partial u_1}{\partial y} \right) dx$$
$$\int_0^1 \phi_2 \left(n_1 \frac{\partial u_2}{\partial x} + n_2 \frac{\partial u_2}{\partial y} \right) dx$$

I want to enforce $u_2 = 0$ and $ \frac{\partial u_1}{\partial y} = g$ on this boundary. Furthermore, the unit normal on this boundary is $n = (0,-1)'$. Using these two pieces of info gives,

$$\int_0^1 -\phi_1 \frac{\partial u_1}{\partial y} dx = \int_0^1 -\phi_1 g dx$$
$$\int_0^1 -\phi_2 \frac{\partial u_2}{\partial y} dx = 0$$

Answer 3

I've updated the code in the main question. There's no longer an error, but it doesn't seem to implement the boundary condition as intended.

Answer 4

e: sorry for the edits, made a mistake.
What's the question exactly? And what is g in your code, the stress vector?

In general, you can use: $\vec T = \vec n\cdot \sigma = (\vec T\cdot\vec n)\vec n + (\vec T\cdot \vec t)\vec t$ in the natural boundary integral.
In your particular case, where the boundaries are aligned with the coordinate
axis, you simply have at the bottom boundary $n = -e_y$ and $t=e_x$. Because of
the Dirichlet BC on the normal component, $\vec v\cdot\vec n = v_y = 0$.
So the boundary integral simplifies to
$$\int\vec n\cdot\sigma = \int f\vec v\cdot\vec t $$, i.e., f*v[0]*ds(1).