How many boundary conditions do we need for the pressure?
 

Hello guys,
I am still somewhat new to CFD (just starting my masters degree) so please bear with me if my question is dumb.

I am simulating fluid flow in a channel for an assignment. The aim of the assignment is to simulate 2 dimensional incompressible Poiseulle flow, moving from left to right. It is a very simple simulation. There are no obstacles and it is a simple rectangular domain. For the velocity, I enforce a velocity inlet condition at the inlet and a homogeneous Neumann condition at the outlet. I enforce a homogeneous Dirichlet condition (no-slip) at the top and bottom wall. This all makes sense to me as the INSE are second order in u .

But, when I was looking over some example tutorial cases in OpenFOAM, I noticed that their simulations were enforcing a boundary condition for the pressure on both the inlet AND the outlet, and correspondingly, on both the bottom AND the top wall. They typically used a homogeneous Neumann condition on the top and the bottom wall, a homogeneous Neumann condition on the inlet, and a Dirichlet condition on the outlet.

This does not make sense to me from a mathematical point of view, because the INSE are first order in p in both the x and y coordinate, so specifying boundary conditions at only, for example, the bottom and left wall should already be enough to fully specify the solution.

I know that the INSE are invariant under constant pressure shifts p --> p+p0 which gives the system a little extra freedom. I also know that typically the pressure is simulated using the pressure Poisson equation, which is of course second order. But this is a consequence of the INSE, not a governing equation in its own right.

I would appreciate some elucidation on this topic. Thanks.

The number and types of boundary condition needed at each boundary depend on the nature of the characteristics at the boundary.

When you learn that a 2nd order PDE needs "only 2 boundary conditions" you are:
1) dealing with an equation that is an independent pde in one dynamic variable (i.e. temperature) that is not coupled with anybody else
2) learning how to solve laplace's equation, where all extrema occur on domain extremes
3) solving on a 1D line or 2D rectangle where there is a clear distinction between x and y coordinates and so you have an inlet and outlet (i.e. left and right boundary) in the x direction and inlet and outlet in the y direction (i.e. top and bottom). And you associate this with "two boundary conditions"

Consider for example that you are solving navier-stokes on a hexagonal domain... clearly "two boundary conditions" is not enough to solve this problem. But hmm, that's strange... suddenly the number of BCs needed depends on the number of boundaries....? Consider also that you can't arbitrarily specify higher 4th order or n-th order derivative on all boundaries to solve this problem. Keep in mind that to a mathematician, an inlet and outlet doesn't mean anything in terms of topology, all boundaries are compact surfaces bounded in . They don't even count how many of them there are.

Thank you for the response. I am familiar with the traditional polygonal domain example. Allow me to be more precise.

Let be a connected, open, bounded, convex subset of , where =2 or 3. (Think of some sort of ball.) Now draw a line through E. Because of the assumptions we made on , this line intersects exactly twice.

Now let's say we are solving some partial differential equation in . As I understand it, for a second order problem or higher, we need to specify two (EXACTLY two) pieces of information, either both at the same point, or one at either point at which this line intersects , in order to uniquely specify a solution (assuming a unique solution exists). The most common way is to use a Dirichlet condition at each point.

However, for a first order problem, (again, as I understand it) we require one, and ONLY one piece of information on this line, so the most we can do is to specify some condition at only one of these points of intersection. So, in your hexagonal domain example, boundary data on only one of the sides is enough to fully specify a solution (matching our intuition).

That is essentially what I am asking about. When we specify two boundary conditions along this line for the pressure, are we just hoping that some solution exists that meets the criteria, and praying that our numerical scheme doesn't diverge?

I hope this frames my question better.

Numerical scheme diverging or not has nothing to do with the structure of the underlying continuous pde. So no need to discuss what happens numerically when you want a general answer to how many BCs you need for a given PDE.

Again you are using an example of an independent equation in E, Navier stokes is a coupled system of non-linear pde's

And if this line enters and exits the same boundary.... what then? You still can have a posedness problem. Oh but you say it can't intersect itself, well a different point on the same "surface" is also not intersecting itself and can be partitioned to be a different surface via "compact surface in R3." I'm not cheating, that's what differential means. I can always carve out your surface to two smol surfaces (i.e. two points) and reduce this argument to specifying conditions at two points and not two surfaces.

Furthermore, the "two" pieces of information you need to provide also aren't arbitrary for posedness. You can overconstrain a problem by specifying E everywhere depending on what the characteristics at each boundary are (the infamous static vs total properties). Consider what happens at a boundary if you have 0 velocity, + infinite velocity, and - infinite velocity. Really what matters is the Mach number, but I hope that drives the point.

I will add that it is straightforward to prove uniqueness in a linear ODE/PDE and there you can easily show that you only need "two boundary conditions" for a 2nd-order ode/pde. But Navier-Stokes is not Laplace's equation. Furthermore, you can use a simple substitution to convert an n-order linear ode/pde into a system of n 1st order ode/pde's, i.e. starting with

define u=y' and then you can rewrite it as




These properties make the analysis easier to generalize to higher orders for linear ODE's and PDE's because if you can show uniqueness for the 1st order case, you can then show it for the 2nd order case and so on by using this trick. But Navier-Stokes is not Laplace's equation.

Let me summarize this by saying all linear 2nd order PDE's look very much alike. Coupled sets of non-linear ones don't necessarily have the same simple rules.

Since you are new in CFD, I suggest to stop a moment using OpenFOAM and just think about a more general statement of the problem.
You have a system of equations, the divergence-free velocity constraint and the momentum equation, they should be solved as a coupled system. The hyperbolic divergence-free constraint is transformed into the elliptic Poisson equation for the "pressure" (that has nothing to do with the thermodynamic pressure) and in order for the problem to be well posed, you have to prescribe the BC along the whole frontier. The BC for the pressure is a direct derivation of the Hodge decomposition, for this reason Neumann BCs for the pressure are prescribed. However, the singularity of the matrix requires that you fulfill with your BC the compatibility condition, that is you fulfill the continuity equation.
Therefore, you are wrong, the pressure equation is a governing equation substituting the continuity.

I don't understand everything you have said, but I gather that the point you are trying to convey is "It's not so simple". I suppose all we can do is assign boundary conditions that match the physics of the problem, and hope for the best. Is that correct?

-----

ADDENDUM: Why this concerns me.

The momentum equation can be rewritten as



Where the right hand side side is curl-free. Let's suppose we somehow knew the solution for ahead of time, and all we have to do is find the pressure that fulfills our prescribed BCs. So, in general, is it possible to find ? Essentially, we focus on the following:

Given a connected, open, bounded subset of (n=2,3) and a curl-free vector field , is it always possible to find a scalar field that satisfies



The answer is NO. For example, consider a two-dimensional case and take . It is easily shown that all of the that satisfy are of the form . Therefore taking for example there is no solution.

Interesting, and thanks for your remarks. Are you essentially saying that since the INSE are second order in u, they are basically second order in p as well by proxy? I am happy enough with this simple explanation if it is true.

In some sense you are right...the formulation for the incompressible flows is somehow both not simple and not physical ... you are assuming that the pressure wave moves at infinite speed, a fact that is not physically meaningful. But it works and provides a good approximation for low Mach flows.



I think you have to start from the theory of incompressible flows to understand that the "pressure" is nothing but a lagrangian multiplier.


I just try to clarify the point in a different way:



Assume dv/dt=a to be the eulerian acceleration and a* to be convective+diffusive terms.


The system writes as


Div a =0 (Continuity)


a + Grad p = a* (Momentum)


Now you have to consider the Helmholtz-Hodge decomposition theorem to understand that



Div Grad p = Div a*



is nothing but the continuity equation and the BCs are deduced by evaluating the momentum equation on the frontier



dp/dn = (a*-a).n



This way you have that the compatibility equation



Int[V] Div Grad p dV = Int[V] Div a* dV

is satisfied.


That allows to have a solution apart an arbitrary function of time. In case the integral relation is not satisfied, your iterative solver will never converge. Fixing a Dirichlet condition is only a way to produce a wrong convergence, it is not necessary at all.

about me
 

FYI everyone, I have a pretty extensive background in mathematics, so please feel free to communicate your ideas to me using any and all mathematical language you wish. I actually will appreciate it far better than an argument using words alone.

You can read my old paper about the topic
https://www.researchgate.net/publica...ary_conditions

I'm not sure I understand your issue but, from the original NSE one can derive a pressure equation which is a Helmholtz equation. Taking the incompressible limit, this equation collapses into a Poisson equation. I won't go into the details of which boundary conditions these equations need.

Incompressible pressure based methods do just the same, solve for a set of equations which has been artificially forced to give an elliptic nature to the pressure equation.

As an additional note, we also have preconditioned density based methods, which work exactly in the opposite direction, by artificially forcing the equations to always have the compressible nature, even in the full incompressible limit, trough an artificial speed of sound. When using such methods, boundary conditions are indeed assigned following the artificially hyperbolic nature of the equations.

Would the compatibility equation

Int[V] Div Grad p dV = Int[V] Div a* dV

get satisfied if the boundary condition for p is given as dp/dn=0 and at the boundary a*.n is simply set equal to a.n?

I'm assuming a*.n would need to be otherwise somehow interpolated to the boundary if cell centered FVM was used.

Yes, the Poisson problem is well posted if, equivalently, you prescribe homogeneous Neumann BC and modify the source term close to the boundaries.

Very interesting! If this kind of a boundary condition is applied to the conventional

u = u*-deltaT*gradP

decomposition, would this boundary condition allow to have nonzero tangential velocities prescribed at boundary aswell?

Or just a nonzero normal component?

Only the normal component can be prescribed. The tangential component is a result of the numerical solution.
This is a known issue in the fractional time step method.

Thanks for the info!

I was also wondering if this kind of a boundary condition would also work with impulsively started flows?

For example if flow in a pipe was initially at rest and and then inlet velocity (say 1 m/s) was given at the the first time step, would solution of the poisson equation still enforce continuity?

Using the incompressible flow model you cannot simulate a case where the flow enters and propagates in the interior. Pressure wave are assumed to travel at infinite speed.
In other words, you need to prescribe an outlet BC to ensure that the mass entering and the mass leaving are impulsively equal. In such a case the Poisson problem is well posed.

So in theory a badly posed Poisson problem at the start would continue being badly posed through subsequent time steps and eventually even result in a solver crash?

Would a potential flow solution then be a reasonable starting point for the solver?

The solver does not converge since from the first time step. A potential initial solution can be used.

So to elaborate more on the proper use of tangential components, does the fundamental issue restrict us to use tangential components from the solution during the momentum predictor step or at the Poisson equation step or perhaps in both?

I'm trying to figure out if and/or when i should be interpolating the tangential components from the interior of the solution in cell centered FVM.

After the projection step you have the correction step where the pressure gradient is added to the intermediate velocity. That could be done also for the tangential component on the boundary and you could check the resulting error, it does not produce exactly the physical tangential velocity. The practice is simply to set the value of the tangential component to the correct physical value.
If the error propagates or not depends on the orthogonality of the decomposition. Generally, a numerical boundary layer is produced.

This is a brief note, a full paper is addressed in reference.
https://www.researchgate.net/publica...ressible_flows

Thanks, I'll check out the paper!

I had a look at that paper which explained many of the questions I had. The treatment of boundary conditions was a bit more nuanced than I expected.

Here's what's my current understanding of the method:
1. Compute the intermediate velocity field v* subject to the boundary conditions v* = v(n+1) +deltaT*gradP(n)
2. After solving equation for the intermediate velocity, the Poisson equation is solved with boundary conditions gradP=0 and at the boundary v*.n = v(n+1).n is set.
3. Velocities are subsequently corrected to be divergence free as v(n+1) = v* - deltaT*gradP(n+1)

So if exact physical velocity with nonzero normal and tangential components is used at the boundary, there is no guarantee that after corrector step 3 the tangential components are satisfied? Does this also imply that there is a possibility that no-slip conditions aren't satisfied in the tangential direction?

That’s correctl, of course you do not use the computed tangential components that could be wrong but you set the physical one for the next time step.

I had some more time to think about the boundary conditions :D

If the time continuous idea of v* = v*(t) would be used in a method where pressure gradient is retained in the intermediate velocity equation, would the correct boundary conditions then be v* = v(n+1) for such a case?

Also wouldn't that mean that the pressure equation with v*.n = v(n+1).n and subsequent projection would still produce a divergence free solution?

You can retain only the old pressure gradient but not the total integral in time of the pressure gradient. Therefore you have some terms to add in the BC for the intermediate velocity.

The Poisson pressure does not need the value of V* on the boundary.

Thanks for the info, I think I see the universality behind the Poisson equation now.

One side note, if div(V)=0 is enforced, would it also require that the exact physical boundary velocity field would fulfill an additional compatibility condition div(ω)=0, when ω = curl(V)?

There is an issue because the vorticity cannot be specified in general as Dirichlet BC. But you can see in my paper how the vorticity can be introduced in the Hodge decomposition along with the potential vector function.

Hello! It's great to hear that you're starting your masters degree in CFD and exploring fluid flow simulations. Your question is not dumb at all, and it's important to seek clarification to deepen your understanding.

In the context of simulating the 2D incompressible Poiseuille flow in a channel, it is indeed true that the incompressible Navier-Stokes equations (INSE) are second order in velocity components (u, v), while they are first order in pressure (p). This distinction has implications for the boundary conditions that need to be specified.

The no-slip condition (Dirichlet condition) you mentioned on the top and bottom walls is appropriate since it ensures zero velocity at these solid walls, reflecting the physical behavior of the fluid. Similarly, the velocity inlet condition (Dirichlet condition) at the inlet defines the desired velocity profile for the Poiseuille flow.

Regarding the pressure boundary conditions, it is common to enforce a homogeneous Neumann condition (also known as a zero-gradient condition) at the outlet. This means that the pressure gradient in the flow direction is assumed to be zero, which is a reasonable assumption for fully developed flow. This condition allows the pressure to adjust naturally throughout the domain while preserving the overall pressure drop imposed by the velocity inlet condition.

As for the pressure boundary condition at the top and bottom walls, applying a homogeneous Neumann condition is a common approach. This condition assumes that the pressure gradients in the normal direction (y-coordinate in your case) are also zero, meaning there are no pressure variations along these walls. This assumption is often reasonable for channels with a sufficiently large aspect ratio, where the flow is primarily influenced by the velocity profile and not the pressure distribution.

You correctly noted that the pressure can be shifted by a constant (p --> p + p0) without affecting the velocity field, as the INSE are invariant under such constant pressure shifts. This is why only specifying boundary conditions for pressure on the bottom and left walls is mathematically sufficient to determine the solution. However, by also enforcing the homogeneous Neumann condition on the top and outlet boundaries, you provide additional information to the solver, which can help stabilize the pressure calculation and ensure better numerical behavior in the simulation.

Overall, the choice of pressure boundary conditions depends on the specific problem and its boundary characteristics. The approach you described, where the pressure boundary conditions are only specified on the bottom and left walls, can indeed be valid in certain cases. However, enforcing homogeneous Neumann conditions on the top and outlet boundaries is a common practice that can help achieve more stable and accurate results, especially in complex flow scenarios.

I hope this explanation clarifies the rationale behind the pressure boundary conditions in the OpenFOAM tutorial cases and provides you with a better understanding. Feel free to ask further questions if you need more information

Just to be clear, this is pure non sense which could only be written by some AI.

Thanks for the tip, I'll see what I can find. Non-primitive formulations can sometimes be a bit unintuitive at first.

If i understood the paper correctly, the vorticity can be solved explicitly from the intial value problem

d(ω)/dt = curl(a*)

where a* is the intermediate acceleration. As result, one ends up with the correct values of vorticity for time step n+1.

Is it then permissible to deduce the exact physical tangential velocities for time step n+1 from the now available vorticity and carry on with the solution strategy as described in earlier posts?

What is your final goal?



First, in order to be well posed the problem, the physical velocity BC.s must be known for any type of approach, they are not deduced from a formulation.


The problem appears because in the Hodge decomposition only one BC is required. Using the vorticity, you cannot prescribe the normal velocity component. That can be done only by prescribing the elliptic problem for the potential vector.

Note that the Hodge decomposition applies at any time but has never to do with a time evolution from tn to tn+1.

The initial condition at t0 allows to do one step for the vorticity but then you have to approximate the vorticity at the walls.


Maybe I can help you more if you give more details about the numerical method you want to implement.

My idea was to get good estimates for the tangential velocity components at outlet boundaries and evolve these through the time somehow using d(ω)/dt = curl(a*). There is usually very little information available beforehand to what these tangential velocity components should actually be at flow outlets. Tangential velocity components at inlets and walls have usually been well known beforehand.

However, I think I had mixed some key concepts between the different decompositions and their limitations. :D

So am I safe to assume that if I'm using the decomposition

v = v*-deltaT*gradP

with gradP=0 and v*.n = v.n on the boundary at Poisson step, that the only correct way to get a proper decomposition would be to set v to have zero tangential velocity components at boundaries? Therefore making the search for tangential velocities reduntant as they must be set to zero anyhow?

I did some small testing using the above decomposition but with nonzero tangential velocities and had the impression that vorticity didn't pass through the outlet boundary, like it did with the usual p=0 and grad(v).n = 0 conditions. That led me to suspect that I had improperly set my tangential velocity components at the outlet thus affecting the vorticity.

Be careful, the Hodge decomposition in this case does not prescribe anything about the tangential velocity component.

After having performed the decomposition with the normal velocity component, from a theoretical point of view you have to compute the tangential component



v.t = v*.t-deltaT*(gradP).t

on the all the boundary. Of course, no one will do that on the boundaries when you have a known physical value for the tangential velocity, for example on a wall. You set the known physical value and go to the next time step.
Have a look to page 59 here:


https://www.researchgate.net/publica...ary_conditions

As far as the specific issue on the outlet, what normal condition you prescribed? Did you get convergence in the Poisson solver without problems?

Thanks for the elaboration, I think I have a more complete understanding now. So in practice the tangential components are a byproduct of the projection method and these can be nonzero but are often discarded and physical ones are described instead.

I'll see if incrementing the tangential components using the method described in the paper has any effect.


I used a zero normal gradient condition for the normal velocity component. If tangential components were determined using second order upwinding (e.g. V = V+distance*gradV) the results were slightly better than just using zero normal gradient for the tangential components.

The Poisson solver didn't have any noticeable problems.

Actually now that i think of it, the boundaries that seem to work properly have one thing in common, they are all steady. The outlet is the only unsteady one. I think i need to double check my solver.

I finally had some success with the Poisson equation with Neumann conditions in a basic 2D backward-facing step case. The case converged nicely and appeared to be physically acceptable. :)

Turns out, my velocity boundary condition wasn't behaving correctly and it let some inflow at the outlet boundary resulting in vorticity.
I was manually modifying the outlet velocities to satisfy global continuity and did not properly handle the possible backflow situation. This resulted in small short circuiting of momentum.

My initial impression is that Neumann conditions result in lower div(V) values in general than with the Dirichlet condition case.

I need to see how well the method scales with respect to time step and test more in 3D when Im back in office after the summer break.

Thanks for the tips!

I had no problem with Neumann conditions, tested both on analytical solutions and on the 3D backward facing step.
On outflow I adopted fully developed flow, that is zero normal derivative for the velocity. Consequently, at the outflow you have the Poisson equation only in a plane.