Recovering pressure from velocity or streamfunction fields

Question

I am interested in 2D channel flow of an incompressible Stokes fluid (Re << 1), with periodic boundary conditions in the x-direction and no-slip at the walls in the y-direction. I have existing code that solves this using a stream function formulation which is tested and works correctly. The details of the numerical implementation shouldn't matter, but for reference I use finite differences in the y-direction, and a Fourier method in the x-direction.

In the stream-function formulation the pressure is eliminated via the incompressibility condition, however occasionally it is useful to examine it in order to gain some insight. To be clear, this pressure is an auxiliary field, and the accuracy of the simulation isn't dependent on it.

My question is how can this pressure be recovered from the stream function object? I also have the raw velocities in case they are easier to work with. This general problem appears to be called the Pressure Poisson Equation (PPE) but I'm having some difficulties understanding it.

My attempt

Begin with the Stokes equation (subject to some force $\bf f$):

$$0 = \eta \nabla^2{\bf v} - \nabla p + {\bf f} $$

Taking the divergence of this, and applying incompressibility, ${\bf 0} = {\bf \nabla}\cdot{\bf v}$, we find:

$$\nabla^2 {p} = \nabla \cdot {\bf f}, \qquad {\rm (1)}$$

i.e., a Poisson equation for the pressure.

Doubly Periodic

In a doubly periodic system (i.e., periodic in x and y) we can easily calculate the pressure (defined up to a constant), using the boundary conditions that gradients in the pressure are subject to periodic BCs: $\nabla p(x, y) = \nabla p(x+L,y)$ etc. This can easily be solved in Fourier space.

Channel flow

This is where I am unsure - in the periodic direction I know the boundary condition, but in the wall normal direction it's not clear what the boundary condition should be. Eq 16 in Ref (excluding advection, and assuming time-independent BCs) suggest that both forces and velocity must appear in the BC:

${\bf n}\cdot\nabla p = {\bf n}\cdot\left(f + \eta \nabla^2 {\bf v} \right)$

This seems consistent with Ref 2 Eq 1.9 (which sets $\bf f = 0$) where the BC reads:

$\partial p/\partial{\bf n} = {\bf n} \cdot \eta \nabla^2 {\bf v}$

Is this an appropriate way to find the pressure in channel geometries, along with Eq. (1) above? If not are there any alternative approaches, perhaps involving the streamfunction directly? Ease of implementation is preferred over efficiency as this is generally a one-off calculation.

Answer 1

So in case if you want to get pressure field there are two variants.

1) Simple one

As you have already solve the problem in terms of stream function ($\psi$), you get velocity field from the definition of stream function:

$v_{x} = \partial \psi / \partial y,\quad v_{y} = -\partial \psi / \partial x$

Then using momentum equation you can obtain field of pressure gradient, in coordinates form:

$\partial p / \partial x = \eta \left(\partial^2 v_{x}/ \partial x^2 + \partial^2 v_{x}/ \partial y^2\right) + f_{x}$

$\partial p / \partial y = \eta \left(\partial^2 v_{y}/ \partial x^2 + \partial^2 v_{y}/ \partial y^2\right) + f_{y}$

As you don't set any boundary conditions for pressure (e.g. pressure at infinity point), pressure field will not be unique. Pressure field will be unique with accuracy to some additive constant. Knowing pressure derivatives in each point of your mesh you can integrate this equations and get pressure field. Before integrating you should set pressure to some value in some particular point of mesh. Integration could be done simply by formula:

$p_{i+1, j} = p_{i, j} + g_{i, j} \Delta x_{i} \quad p_{i, j+1} = p_{i, j} + g_{i, j} \Delta y_{j}$

2) Another one

Well, you can solve elliptic problem to find out the pressure field. As you mentioned, pressure equation is obtained after taking divergence of the momentum equation. The Neumann boundary condition in case of your geometry is following:

$p_{x} = \eta \partial^2 v_{x} / \partial x^2$

on the $x = 0, x = L$ lines.

Answer 2

The answers you've provided are valid, however, if you're looking for a simpler approach that is specific to your problem, then you may use $\partial p / \partial n = 0$ on a no-slip, no-flow through wall.

The justification behind this, I believe, is that the Dirichlet conditions imply velocity is known and, if steady, means no net forces exist at the boundary. Looking at the momentum equation on the boundary normal to the wall, you'll see that all forces are zero (since the boundary is kinematic). Therefore $\partial p / \partial n = 0$ is a valid BC.