Calculation of pressure in a domain with no-slip boundaries

[parfenyev]

Hello everyone!

I'm trying to understand how pressure is calculated when solving the Navier-Stokes equation in a region limited by no-slip walls. In the documentation, I read that to solve the equation one uses the fractional step method. When using this method, an auxiliary quantity arises, which is designated in the documentation as $p_{non}$, see equations (7) and (8). However, the literature explains that this quantity is not pressure, which initially appears in the incompressible Navier-Stokes equation, see Sec. 3.1 in Ref. 1 and Ref. 2. I have the feeling that the quantity called non-hydrostatic pressure in model.pressures is precisely this auxiliary field $p_{non}$ and it cannot be used to analyze physical pressure near boundaries, since it satisfies wrong boundary conditions, see Eq. (19) in Ref. 1 or Eqs. (9) and (10) in Ref. 2. Could you please comment on this issue?

[glwagner]

From just a quick read, I believe our method is consistent with Ref 1 (but there are some subtleties which I discuss below). Our algorithm sets the predictor velocity wall-normal velocity on boundaries equal to the physical velocity:

Oceananigans.jl/src/Models/NonhydrostaticModels/pressure_correction.jl

Line 13 in 00006c1

fill_halo_regions!(model.velocities, model.clock, fields(model))

and then we indeed solve for pressure prescribing the wall-normal component of the gradient to zero as in Ref 1, equation 19.

The subtlety is that Ref 1 presents eq 15 and 19 as if they are written in stone, but in fact these are a choice. Eq 15 and 19 must be consistent with one another. However, one is free to choose the "boundary conditions" on the predictor velocity in eq 15 arbitrarily. (Introducing the additional predictor velocity field also introduces an additional DoF in the mathematical problem which must be determined by a choice and is not determined by the problem itself.)

It is particularly convenient though to use the choice in eq 15, because then it leads to the convenient boundary conditions on the pressure in eq 19, for which we have nice numerical methods.

I think there is an omission in ref 1 for the case that the inflow changes in time. (We also do not treat this situation correctly... something we should look into.) But I could be wrong. @simone-silvestri is probably interested in this discussion...

A long, long, long, long time ago we did convergence tests to show empirically that our pressure solve was implemented correctly. The particular case of "Forced Flow, fixed slip" is the one of interest:

https://github.com/CliMA/Oceananigans.jl/blob/main/validation/convergence_tests/src/ForcedFlowFixedSlip.jl

Unfortunately I don't think that code will run out of the box, someone would have to put the effort in to update it. But if we can resurrect it, I believe we could prove that the analytical solution (including the pressure field) converges at the expect rate with increasing grid resolution.

Are you diagnosing a problem based on the results of a simulation, or by reading the docs? It's possible the docs are wrong or misleading and we can fix that.

[parfenyev]

So, your point is that for a particular choice (15) the solution of (14)-(15) is discontinuous at the boundary and it produces the delta-function term in right-hand side of (18) that effectively changes the BC (19). I need some time to think about this and see if $\Pi$ can really be equated with pressure. Thank you!

[glwagner]

Yes, no problem! I was wondering when this might come up...

Here is a bit more detail after discussing with a few others like @navidcy.

The Navier-Stokes equation is

$$ \partial_t u = G - \nabla p $$

The predictor velocity is defined such that

$$ \partial_t u_\star = G $$

with $u_\star(t_n) = u$. (This holds up to the next time, $t_{n+1}$).

For simplicity consider impenetrable boundaries with $n \cdot u(x=x_{\Omega}) = 0$ where $n$ is the unit vector normal to the boundary. Taking $n \cdot$ the first equation gives a condition for $\nabla p$ on the boundary,

$$ n \cdot \nabla p = n \cdot G $$

Taking $n \cdot$ the predictor equation gives

$$ \partial_t \left ( n \cdot u_\star \right ) = n \cdot G = n \cdot \nabla p $$

This shows that $n \cdot u_\star$ is not, in fact, equal to $n \cdot u$ on the boundary --- at least not when there is a pressure gradient across the boundary. (This seems to be where the confusion starts for that reference 1.) This has a physical meaning --- impenetrability is in fact enforced by pressure for the physical velocity. The predictor velocity, which obeys an equation that lacks a pressure gradient, cannot be constrained to impenetrability on the boundary; it's missing pressure gradient force that would have that effect.

Now, taking the divergence of the Navier-Stokes equation yields $\nabla^2 p = \nabla \cdot G$; and the divergence of the predictor equation yields

$$ \partial_t \left ( \nabla \cdot u_\star \right ) = \nabla \cdot G = \nabla^2 p $$

Integrating this in time, and using a fully-implicit discretization for $p$,

$$ \int_{t_n}^{t_{n+1}} p , dt \approx \Delta t p^{n+1} $$

yields the pressure Poisson equation,

$$ \nabla \cdot u_{\star} = \Delta t \nabla^2 p^{n+1} $$

Note that we have used $u_{\star}(t_n) = u(t_n)$ so that $\nabla \cdot u_{\star}(t_n) = \nabla \cdot u(t_n) = 0$. To find the boundary condition for $p^{n+1}$ we do the same with the predictor equation on the boundary, using $n \cdot u_{\star}(t_n) = n \cdot u(t_n) = 0$,

$$ n \cdot u_{\star} = \Delta t n \cdot \nabla p^{n+1} $$

That's the boundary condition on the pressure Poisson equation.

This boundary condition is effectively implemented in the numerical solution of the Poisson equation by "filling halo regions" (ie spuriously enforcing impenetrability on the predictor velocity field prior to calculating the RHS of the pressure Poisson. Enforcing impenetrability moves the application of the boundary condition above into the RHS of the pressure Poisson equation.

Moving the contribution of inhomogeneous boundary conditions into the RHS of the Poisson equation is a typical technique, for example:

It is rather confusing and tricky that the way we incorporate the correct pressure Poisson boundary conditions effectively is by enforcing the "same" boundary conditions on the predictor velocity as on the physical velocity. Mathematically, the predictor velocity does not satisfy the same boudnary conditions as the physical velocity. Nevertheless, this "trick" is quite a convenient way to effectively incorporate the inhomogeneous boundary contributions into the pressure Poisson equation.

I personally believe that most authors of CFD codes do not realize that this is effectively what is happening. It puzzled me for a while. But when I worked out the detailed discrete math to understand what I needed to add to the pressure Poisson equation to correctly satisfy the boundary conditions, I found miraculously that the terms I needed to add were already there. That lead to a realization that enforcing impenetrability on the predictor field is effectively how the pressure Poisson boundary conditions are satisfied.

It's also still true that due to the crude first-order time-discretization, the numerical pressure field is only accurate to O(dt).

[parfenyev]

After the break, I returned to this topic again and realized that my misunderstanding was deeper and was not related to the numerical method at all. Let me consider the 2D Navier-Stokes equation for the incompressible fluid in a no-slip domain:

$\partial_t \vec v + (\vec v \cdot \nabla) \vec v = - \nabla p + \nu \nabla^2 \vec v,$
$\mathrm{div} \ \vec v = 0, \qquad \vec v \vert_{\partial G} = 0.$

Taking the divergence, we come to the Poisson equation for the pressure:

$\nabla^2 p = -\partial_i (v_j \partial_j) v_i, \qquad (A1)$

where we sum over the repeated indices. To solve this equation we need a boundary condition. To obtain BC you actually suggested to consider the Navier-Stokes equation in the immediate neighborhood of a boundary:

$\nabla p \vert_{\partial G} = \nu \nabla^2 \vec v \vert_{\partial G}$

where we have used $\vec v \vert_{\partial G} = 0$. Actually, it is a vector equation and so the problem is over-determined. In the numerical scheme, you project this equation on the direction normal to the boundary

$\vec n \cdot \nabla p \vert_{\partial G} = \vec n \cdot \nu \nabla^2 \vec v \vert_{\partial G}, \qquad (A2)$

and then solve the Poisson equation (A1) for the pressure with the Neumann BC (A2), am I right?

But in principle, you can project the vector equation on the other direction, e.g. on the tangential to the boundary vector:

$\vec \tau \cdot \nabla p \vert_{\partial G} = \vec \tau \cdot \nu \nabla^2 \vec v \vert_{\partial G}. \qquad (A3)$

Since the pressure is determined up to an arbitrary constant, you can integrate (A3) along the boundary and obtain the corresponding Dirichlet BC. Do you know how to prove that the solutions to equations (A1, A2) and (A1, A3) are identical?

Trying to find an answer to this question, I found a paper that states that, generally speaking, the solutions are different. This is surprising because it is then not clear which boundary condition should be used. And the question here is not so much about the numerical scheme, but about the general understanding of the Navier-Stokes equation.

[francispoulin]

I will let @glwagner respond too this in full but I will share some thoughts.

First, I'm not convinced that the problem is over determined since your equation for the gradient of pressure, between A1 and A2, is not independent of the other equations. It is the momentum equations evaluated at the boundary and simply what we need in order to know what the gradient of the pressure is on the boundary.

Second, I can't read the paper because the text comes out all weird but I'm not surprised that you can get different answers by using the different boundary conditions as they are specifying different things. A2 specifies the directional derivative of the pressure normal to the boundary and A3 specifies the directional derivative of the pressure tangent to the boundary (if I understand the notation correction).

I think that A2 is the more natural equation to use since the RHS can be written in terms of the dot product of the velocity and the normal unit vector, what the no slip boundary condition specifies.

[parfenyev]

I can't read the paper because the text comes out all weird.

You can find the normal version of the paper here, sorry for the problem.

I'm not surprised that you can get different answers by using the different boundary conditions as they are specifying different things.

But pressure is a physical quantity that can be measured experimentally. What boundary condition should I use to get an answer consistent with the laboratory experiment and why? Actually, I do not understand your point. If the answers are different, then one of them (or both) is incorrect.

[parfenyev]

I have found an answer to my question in Ref. 1. In short, it is incorrect to use the equation between (A1) and (A2) to generate BC. Also, I have found a useful Ref. 2 about the accuracy of fractional step method. In particular, it addresses the difference between pressure and pseudo-pressure, and also determines the accuracy of fulfillment of the boundary conditions for velocity. The basics of the fractional step method are well explained in Ref. 3. I wish to thank @glwagner and @francispoulin and close the discussion.

[glwagner]

Awesome!

Also note that we have performed some convergence tests that demonstrate the accuracy of the pressure solver when there is a non-zero pressure gradient across the boundary (eg comparing to analytical solutions). The tests are probably stale but can be resurrected and are in here: https://github.com/CliMA/Oceananigans.jl/tree/main/validation/convergence_tests

In particular the "forced flow, fixed slip" test implements a forced problem with a "fixed slip" boundary condition on the velocity (requiring the velocity to match a time-dependent function). I think these tests are fairly complete in that they test all directions / different velocity components, albeit always in a 2D configuration.