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.