Poiseuille flow

Question

I have probably a very stupid problem. I can't solve a simple Poiseuille flow in a straight 2-D channel driven by a pressure drop. Results are complete nonsense. (Setting zero pressure on the outflow gives trivial solution/no flow. Taking pressure on the outflow non-zero (but still having some pressure drop) gives a solution flowing in to the channel from both ends). Also satisfying of the incompressibility is very vague. Could you, please, look at my code? I have no idea what I am doing wrong.

Marek

from dolfin import *
mesh = Mesh("P.xml")

V = VectorFunctionSpace(mesh,"CG",2)
Q = FunctionSpace(mesh,"CG",1)
W = V*Q

w = Function(W)

nu = 100
L = 20

noslip = DirichletBC(W.sub(0),(0,0),"on_boundary && (x[1] > 4.0 - DOLFIN_EPS | x[1] < DOLFIN_EPS)")
inflow = DirichletBC(W.sub(1),10,"on_boundary && x[0] < DOLFIN_EPS")
outflow = DirichletBC(W.sub(1),2,"on_boundary && x[0] > L - DOLFIN_EPS")

bcs = [noslip,inflow,outflow]

(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)

f = Constant((0,0))
a = nu*inner(grad(u),grad(v))*dx + p*div(v)*dx + div(u)*q*dx
L = inner(f,v)*dx

problem = LinearVariationalProblem(a,L,w,bcs=bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
(u,p) = split(w)

plot(p)
plot(div(u))
plot(u)
interactive()

Answer 1

This actually works out fine if you don't integrate the pressure term by parts and write your variational problem as follows:

a = nu*inner(grad(u),grad(v))*dx + dot(grad(p), v)*dx + div(u)*q*dx

You also also need write the actual numerical value (not L) in your definition of the outflow boundary. Your variable L will otherwise not be used. I'm surprised you didn't get an error message that L is undefined.

Answer 2

Part of your problem is your choice of boundary conditions. This paper (which I wrote) develops the criteria for imposing a pressure boundary condition on the incompressible Navier-Stokes equations. The short version is that you need to impose the normal component of the normal traction on the inlet and outlet to be the pressures you want on each surface, and also force the flow to be normal to those surfaces. The first condition is achieved by through a natural boundary condition, and the latter can be done through a Dirichlet condition if the surfaces are axis-aligned and flat (otherwise you'll have to use a Largrange multiplier, penalty, or constraint).

In a typical Taylor-Hood FEM formulation you should not impose Dirichlet conditions on the pressure. It looks to me like you're trying to use these elements (piecewise quadratic velocities and piecewise linear pressures), but I'm no expert in what I'm guessing is FENICs (and so I can't help you with the code itself).