Is it possible to enforce Non-slip boundary condition with BDM elements?

Hello everyone,

I have a problem to solve Stokes problem using BDM elements. I am trying a very simple approach to validate my mode, a steady-state flow of a fluid between two parallel infinite plates. I believe that the error is relate to the non-slip boundary conditions on the wall for DBM elements. I use this approach:

bc2 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundaries, 2)
bc4 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundaries, 4)

Can I apply the non-slip boundary conditions on the wall for BDM elements? How Can I do it?

Here is my strong and weak formulation for this problem:

- \mu \nabla ^{2} \textbf u + \nabla p = 0

\nabla \cdot \mathbf{u} = 0 ,

\mathbf{u} = 0 , \partial \Omega _{wall}

where \mu is the dynamic viscosity, \textbf{u} and p are the velocity and pressure respectively.

The weak form is:

\int_{\Omega}^{}\nabla \mathbf{u} \cdot q dx= 0

\int_{\Omega}^{} \mu \nabla \textbf u \cdot \nabla v \: dx - \int_{\Omega}^{}p \cdot \nabla v \:dx + \int_{\Gamma}^{}n \cdot vp \: ds + Stab= 0
Stab = \mu \left( \sum_{E \:\epsilon \: \varepsilon }^{} \frac{\alpha }{h}\left \langle \left [ \mathbf{u} \right],\left [ v \right ] \right \rangle_{E} - \left \langle \left \{ \nabla\mathbf{u} \right \},\left [ v \right ] \right \rangle_{E} - \left \langle \left \{ \nabla v \right \},\left [ \mathbf{u} \right ] \right \rangle_{E}\right)

where \left [ \: \right ] and \left \{ \right \} are the jump and average operators, h is the cell size and \alpha is a positive constant on each edge.

My code is:

from dolfin import *
import ufl


Nx = 50
Ny = 50

mesh = UnitSquareMesh(Nx, Ny)
n = FacetNormal(mesh)

order = 1
V = FiniteElement("BDM", mesh.ufl_cell(), order)
Q = FiniteElement("DG", mesh.ufl_cell(), order - 1)
Element = MixedElement([V, Q])
W = FunctionSpace(mesh, Element)


boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))


left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

ds = Measure("ds", domain=mesh, subdomain_data=boundaries)


bc2 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundaries, 2)
bc4 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundaries, 4)

bcs = [bc2, bc4]

(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
mu = 1
pin = 10  # 1.0
pout = 0.0

alpha = 35
h = CellDiameter(mesh)
h2 = ufl.Min(h("+"), h("-"))


def tensor_jump(v, n):
    return ufl.outer(v, n)("+") + ufl.outer(v, n)("-")


stab = (
    mu * (alpha / h2) * inner(tensor_jump(u, n), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(u)), tensor_jump(v, n)) * dS
    - mu * inner(avg(grad(v)), tensor_jump(u, n)) * dS
)

a = mu * inner(grad(u), grad(v)) * dx - div(v) * p * dx + div(u) * q * dx + stab
L = inner(f, v) * dx - pin * dot(v, n) * ds(1) - pout * dot(v, n) * ds(3)

w = Function(W)
solve(a == L, w, bcs)

(u, p) = w.split(True)

Any suggestions would be very helpful.

You could use Nitsche’s method. Assuming you still enforce the no-penetration constraint in the normal direction strongly (i.e., with a DirichletBC), Nitsche’s method would be analogous to how you’re already weakly enforcing tangential velocity continuity between elements, but with the following replacements:

\lbrack\mathbf{v}\rbrack~\to~\mathbf{v}\otimes\mathbf{n}~~,~~\lbrack\mathbf{u}\rbrack~\to~(\mathbf{u} - \mathbf{g})\otimes\mathbf{n}~~,~~\{\nabla\mathbf{u}\}~\to~\nabla\mathbf{u}~~,~~\text{and}~~\{\nabla\mathbf{v}\}~\to~\nabla\mathbf{v}~~,

where \mathbf{g} is the prescribed boundary velocity (e.g., \mathbf{g} = \mathbf{0} for no-slip). This can be thought of as “DG at the boundary”.

2 Likes

Dear @kamensky, this solution works. Thank you

My only concern is that in the boundary the velocity is not zero. Values are between 10^{-4} and 10^{-6} and in some points negative values ( -4 \cdot 10 ^ {- 6} ) depending on the mesh refinement like in the figure. Is it
expected ?

The piece of code add:

nitsche = (
    alpha / h * inner(outer(v, n), outer(u, n)) * ds(2)
    - inner(grad(u), outer(v, n)) * ds(2)
    - inner(grad(v), outer(u, n)) * ds(2)
    + alpha / h * inner(outer(v, n), outer(u, n)) * ds(4)
    - inner(grad(u), outer(v, n)) * ds(4)
    - inner(grad(v), outer(u, n)) * ds(4)
)

a = (
    mu * inner(grad(u), grad(v)) * dx
    - div(v) * p * dx
    + div(u) * q * dx
    + stab
    + nitche
)

Weak enforcement of boundary conditions (e.g., Nitsche’s method) doesn’t satisfy the boundary conditions exactly on a finite mesh, but the error should converge to zero as h\to 0.

As an aside, I’ll mention that this inexactness is actually an advantage in many situations. It’s not so important in diffusion-dominated problems like Stokes flow, but when the exact solution contains sharp layers at the boundary (as in convection-dominated flow problems), relaxing the Dirichlet boundary conditions can actually give much more accurate solutions on the interior of the domain. The analogy to explicit wall models for turbulent flow is discussed in detail here.

1 Like

Thanks for the reference @kamensky