Stokes flow simulation by imposing pressure bc does not give sensible results!

Question

I did some Stokes flow simulation with and without stabilization by imposing pressure bc for the inlet and outlet but the results does not seem to be reasonable. It gives me high velocity close to the boundary but very small values for other parts of the mesh. I tested different codes in fenics demo for stokes flow and they gave me similar results. But if I use inlet velocity and some pressure value for the outlet then the results are fine. Here also you can find my code:

from future import print_function
from dolfin import *

Load mesh and subdomains

mesh = UnitSquareMesh(50, 50)
set_log_level(PROGRESS)

Define function spaces

V = VectorFunctionSpace(mesh, "CG", 2)
P = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, P]) #Taylor-Hood

No-slip boundary condition for velocity

x1 = 0, x1 = 1 and around the dolphin

noslip = Constant((0, 0))
inflow = Expression(("100", "0"))
zero = Constant(0)

Inflow boundary condition for velocity

x0 = 1

def inlet_boundary(x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[0]) < tol

def outlet_boundary(x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[0]-1) < tol

def wall_boundary(x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]) < tol or abs(x[1]-1) < tol

Collect boundary conditions

bc_1 = DirichletBC(W.sub(0), noslip, wall_boundary)

Inflow boundary condition for velocity

bc_2 = DirichletBC(W.sub(1), Constant(0.0), inlet_boundary)
bc_3 = DirichletBC(W.sub(1), Constant(10.0), outlet_boundary)

Define variational problem

(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
a = (inner(grad(u), grad(v)) - div(v)p + qdiv(u))dx
L = inner(f, v)dx

Compute solution

A, b = assemble_system(a, L)
bc_2.apply(A,b)
bc_1.apply(A,b)
bc_3.apply(A,b)

w = Function(W)
solve(A, w.vector(), b,'lu')
u, p = w.split()

# Split the mixed solution using a shallow copy

Save solution in VTK format

ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p

Answer 1

Hi,
Pressure conditions have to be imposed weakly as a Neumann condition.
The stokes demos show you how to do this.

Comments

Hi,

Thanks for your response. I know how to impose Neumann bc, but for stokes flow Neumann bc is "grad(u).n + p.n =g" which is a combination of velocity gradient and pressure. Is there any way to impose only pressure as a bc?

---

It is common to say that you set the pressure, but actually you set the stress.
this is the only way. However, pressure is usually the dominating part of
the stress (or normal stress).

---

Yes, you are right. But in some commercial software such as COMSOL it has one option for imposing pressure and another option for imposing normal stress.

---

If you have a operator splitting scheme for a flow solver you may assign also pressure.
It is however not correct from a mathematical point of view. But for a Stokes
solver like the one you propose here, you should do as already stated. Unless you
want to dwell into deeper mathematical issues with Lagrange multipliers etc which
may potentially give you more freedom. I am not sure.

---

Thanks, got it.