Strange Behavior of Pressure Solution to Stationary NS-equations

Question

Hi all,

I'm attempting to compute solutions to the 2D stationary NS equations in a unit square. I prescribe no-slip BC on the top and bottom walls, and periodic BC on the left/right walls. To drive the fluid, I prescribe a pressure drop from the left to right. (p = 1 on the left side, p = 0 on the right). When I plot the computed solution, I get strange results for the computed pressure.

For UnitSquareMesh(16,16), I get pressure values that are O(10^3). Mesh refinement brings the value down; for instance, on UnitSquareMesh(64,64) I get ~O(10^2). But, further mesh refinement (UnitSquareMesh(128,128)) computes pressure to be O(10^5). What is strange to me is that:

(i) I'm in a low Re regime--I have viscosity set to 1 .
(ii) No matter how fine or coarse the mesh, I still get the expected "Poiseuille flow" i.e. a parabolic wavefront.
(iii) With the exact same problem setup, solving the Stokes equation (i.e. dropping the convective term and solving a linear system) always gives me (essentially) the same pressure values, no matter how fine the mesh.

Can anyone help me explain this behavior? Attached below is my code (modulo defining domain subboundaries). Any help is much appreciated.

#Create no-slip BC for velocity and pressure BCs
no_slip = Constant((0.0, 0.0))
p_right = Constant(0.0)
p_left = Constant(1.0)
bc0 = DirichletBC(W.sub(0), noslip, dbc)
bcs = [bc0]

# Define new measures for exterior boundaries
n = FacetNormal(mesh)
ds = Measure("ds")[boundaries]

# set parameter values -- viscosity
nu = 1.0

# Define variational problem
(v, q) = TestFunctions(W)
w = Function(W)
u, p = split(w)

F = nu*inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx + dot(grad(u)*u, v)*dx \
    + p_left*dot(v,n)*ds(1) + p_right*dot(v,n)*ds(2)

solve(F == 0, w, bcs)

u, p = split(w)
plot(u)
plot(p)
plot(mesh)

interactive()

Answer 1

Hi,

Since you are imposing Dirichlet boundary conditions on the top and the bottom walls and periodic boundary conditions on the right and left walls, the pressure is defined up to an arbitrary constant (Only $\nabla p$ and not $p$ appears in the N-S equations).

To fix the constant value for the pressure you can:

• Use a lagragian multiplier to enforce $\int_\Omega p = 0$.

• Use a pointwise dirichlet condition to enforce p = 0 at a particular point at the domain.

  def origin(x,on_boundary):
  return x[0] < dl.DOLFIN_EPS and x[1] < dl.DOLFIN_EPS
  bc_pressure_at_point = dl.DirichletBC(W.sub(1), dl.Constant(0), origin, 'pointwise')

Comments

This is true--the pressure that appears in the N-S equations is only defined up to a constant. However, by enforcing that pressure on the left wall must be 1.0, and pressure on the right wall must be 0.0, does that not have the effect of your second suggestion? Enforcing the pressure on an entire boundary certainly enforces it at a point.

Using a Lagrange multiplier to ensure mean zero pressure is something I will try though. Thanks for the input.

Regards,
-Sean

---

Hi,

Since you imposed periodic boundary condition, you don't have boundary integral

$$ \int_\partial\Omega -p \cdot \mathbf{n} + \boldsymbol{\sigma}_n. $$

Therefore the integrals:

+ p_left*dot(v,n)*ds(1) + p_right*dot(v,n)*ds(2)

are not forcing p to be equal to 1 on the left boundary and 0 on the right boundary in the momentum equation.

You are simply imposing a concentrated (dirach delta) source on the periodic boundary.

---

Well I want to impose periodicity on the velocity, not the pressure. From a mathematical point of view, prescribing $$p = p_{left}$$ on $$(x=0) \times [0,1]$$ and $$p = p_{right}$$ on $$ (x=1) \times [0,1]$$ means that when you integrate by parts the term

$$ \int_{\Omega} \nabla p \cdot v dx$$

you get

$$ - \int_{\Omega} p \nabla \cdot v dx + \int_{\Gamma_{in}} p_{left}v\cdot n \,\,\,ds + \int_{\Gamma_{out}} p_{right} v\cdot n\,\,\, ds $$

where $$\partial \Omega = \Gamma_{rigid} \cup \Gamma_{in} \cup \Gamma_{out}$$, and the boundary integrals over $$\Gamma_{rigid}$$ vanish as a consequence of the no slip condition.

Does this make sense? I have a feeling the problem with my code is that when I define my function spaces:

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

I get an error if I don't include the constrained_domain = PeriodicBoundary() on the pressure space, Q. Is there any way to define my mixed function space without making both of them periodic?

Best,
-Sean

---

Hi Sean,

The problem is that after you define a periodic boundary you have

$$ \partial \Omega = \Gamma_{rigid}, $$

And therefore:

$$ \int_\Omega \nabla p \cdot v = - \int_\Omega p \nabla \cdot v + \int_{\Gamma_{rigid}} p n \cdot v. $$

This shows why the two boundary integral you wrote play the role of concentrated forces in the momentum equation, but do not set the value for the pressure at $\Gamma_{in}$ or $\Gamma_{out}$.

Hope this explains better my previous point.

Best,

Umbe