This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Prescribed normal velocity and tangential traction

+2 votes

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()
asked Jan 19, 2017 by radbiv_kylops FEniCS Novice (680 points)
edited Jan 19, 2017 by radbiv_kylops

HI,

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

Best Regards!

Hamilton

4 Answers

0 votes

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

answered Jan 19, 2017 by SM FEniCS Novice (630 points)

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

+1 vote

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$$

answered Jan 19, 2017 by radbiv_kylops FEniCS Novice (680 points)
0 votes

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.

answered Jan 19, 2017 by radbiv_kylops FEniCS Novice (680 points)
edited Jan 19, 2017 by radbiv_kylops
0 votes

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

answered Mar 6, 2017 by dajuno FEniCS User (4,140 points)
edited Mar 6, 2017 by dajuno
...