[wildfire230]

Hi all,


I have done a lot of searching to find an answer to this question, but unfortunately I haven't been very successful.


It seems to be very common in the OpenFOAM community to use the zeroGradient boundary condition for pressure on all solid walls in most simulations. However, it is obvious from the Navier-Stokes equations that this boundary condition cannot be correct in general. I realize it is a reasonable approximation for certain types of flows, but it is certainly not general.


In fact, for some flows, such as low Reynolds number flows past obstacles, the errors due to zeroGradient pressure boundary conditions could be the same order of magnitude as the solution itself.


My question is, what should we do? What is the analytically correct pressure boundary condition on channel walls, and has it been successfully implemented into OpenFOAM?


I am specifically interested in the zero Reynolds number, Stokes flow limit, but I am also interested in the general case. For Stokes flow I have tried using groovyBC as follows:

type groovyBC;
valueExpression "0";
value uniform 0;
gradientExpression "normal() & lapU";
fractionExpression "0";

where lapU is calculated in my solver as fvc::laplacian(U).
I'm not 100% sure if this is working correctly or not, but it seems to make the solver less stable compared to simply using zeroGradient.


I would greatly appreciate any thoughts or advice on this issue. It seems to be a key question regarding the use of OpenFOAM which I have not seen a satisfactory answer for.


Best regards

[wildfire230]

Any ideas? I've been trying to just write a boundary element method solver from scratch, but it would be much more convenient to figure this out with OpenFOAM.

[wildfire230]

Sorry for the bump, I still haven't had any luck finding an answer.

[RobertHB]

Originally Posted by wildfire230

Sorry for the bump, I still haven't had any luck finding an answer.

Hey,
i'm dealing with Stokes flows myself and i gotta admit that i did not think about differences in the pressure variable or BC treatment until now. What i found is that the velocity field corresponds perfectly to analytical profiles even with the p zeroGradient BC. That was/is good enough for me.

[wildfire230]

Thanks for the response. Out of curiosity, does your flow have any fine geometric features like flow through a narrow gap or something like that? My suspicion is that these errors might be a problem when the velocity gradients are large.

[RobertHB]

My inlet and outlet are cylinders with a diameter of 0.8 and 0.7mm and the domain in 0.5mm heigh. Nothing narrower in my case.

[christiank_ucl]

Hi,
I was wondering if you have looked further into this problem?

I am running similar simulations- Stokes flow past confined cylinders. I have run the boundary conditions dp/dn=0 and the boundary condition you mention dp/dn=lap(u).n and the velocity fields are exactly the same but dp/dr is different on the cylinder. When setting the boundary condition dp/dn=0, and then extracting the normal gradient of the pressure on the cylinder I found this to be a non-zero value - I am not sure if this is because of how the value is calculated or is some other numerical artefact.
Best regards,
Christian

[wildfire230]

To be honest, I have not thought about this problem much in a while. If you are setting dp/dn=0, and then you calculate dp/dn and it is not zero, that must be some numerical artifact though I am sure.

[Santiago]

Originally Posted by wildfire230

Hi all,

It seems to be very common in the OpenFOAM community to use the zeroGradient boundary condition for pressure on all solid walls in most simulations. However, it is obvious from the Navier-Stokes equations that this boundary condition cannot be correct in general. I realize it is a reasonable approximation for certain types of flows, but it is certainly not general.

This statement is not true. Mathematically speaking, the projection method for the solution of the incompressible N-S equations REQUIRES NO BOUNDARY CONDITIONS for the pressure Poisson equation, since all weak formulations only require p to be L2 continuous (not H1, or sobolev, hence no boundary data). SEE HERE FOR DETAILS. You only need that the initial velocity field is divergence-free EVERYWHERE, including the boundaries.

However, when implementing SEGREGATED solver in real life, well, you get boundary conditions for p. Using a Natural BC (zeroGradient) or a projection of the momentum eqn on the wall (not like you proposed) for the pressure just affects the regularity (smoothness) of the solution. In low-order methods this shouldn't produce any problems, and the effects on turbulence quantities are negligible, so long you keep your deltaT around the kolmogorov time scale and you use a good turbulence model, etc...

[christiank_ucl]

Hi Wildfire,

Thank you for your quick reply. I was wondering when you were working on it if you tried different ways of implementing the boundary condition as I've tried the way you describe:

{
type groovyBC;
valueExpression "0";
value uniform 0;
gradientExpression "normal() & Lap_U";
fractionExpression "0";

}

and I am seeing quite odd results. For example, when I refine the mesh, dp/dr on the cylinder does not converge to the analytical solution.



I also tried the boundary condition



{

type groovyBC;
valueExpression "0";
value uniform 0;
variables "grad_temp=-snGrad(Lap_U);"
gradientExpression "grad_temp;"
fractionExpression "0";
}


but the above is not giving results close to the other implementation or the analytical solution.


If you had any thoughts or experiences to share about these points, that'd be great!
Thanks,
c

[Santiago]

Originally Posted by christiank_ucl

Hi Wildfire,

Thank you for your quick reply. I was wondering when you were working on it if you tried different ways of implementing the boundary condition as I've tried the way you describe:

{
type groovyBC;
valueExpression "0";
value uniform 0;
gradientExpression "normal() & Lap_U";
fractionExpression "0";

}

and I am seeing quite odd results. For example, when I refine the mesh, dp/dr on the cylinder does not converge to the analytical solution.



I also tried the boundary condition



{

type groovyBC;
valueExpression "0";
value uniform 0;
variables "grad_temp=-snGrad(Lap_U);"
gradientExpression "grad_temp;"
fractionExpression "0";
}


but the above is not giving results close to the other implementation or the analytical solution.


If you had any thoughts or experiences to share about these points, that'd be great!
Thanks,
c

Have you tried using zeroGradient?

[Santiago]

Originally Posted by christiank_ucl

Hi,
I was wondering if you have looked further into this problem?

I am running similar simulations- Stokes flow past confined cylinders. I have run the boundary conditions dp/dn=0 and the boundary condition you mention dp/dn=lap(u).n and the velocity fields are exactly the same but dp/dr is different on the cylinder. When setting the boundary condition dp/dn=0, and then extracting the normal gradient of the pressure on the cylinder I found this to be a non-zero value - I am not sure if this is because of how the value is calculated or is some other numerical artefact.
Best regards,
Christian

Where are you post-processing your results/calculating grad(p)?

[christiank_ucl]

Hi Santiago,
I am looking at Stokes flow past a confined cylinder.

Yes, Ive tried dp/dn=0 and also the gradient conditions I described above. For all the pressure gradient boundary conditions, there is no difference in the velocity field or the pressure on the cylinder surface but dp/dr on the cylinder is different.

I have tried two ways of calculating the pressure gradient on the cylinder:
1. Using the foamCalcEx utility to get dp/dx and dp/dy and then taking these values out on the cylinder surface and then calculating dp/dn=dp/dx cos theta + dp/dy sin theta.
2. Using a utility similar to dT/dn for heat flux but changing the variable to pressure.



I have not found consistent values for the two methods. Also, I have not found the values converge when I carried out a mesh independence study.


The main reason for using the full boundary condition (consistent with NS equations) is that even if it is possible (as in there is no effect on the velocity field) to use dp/dn=0 for Stokes flow, I am not sure if dp/dn=0 is valid when the Reynolds number is increased.


Thank you for any advice/thoughts/experience,
Christian

[Santiago]

Originally Posted by christiank_ucl

I am not sure if dp/dn=0 is valid when the Reynolds number is increased.

As I explained in earlier posts, setting up either natural boundary condition, or projecting the momentum equation onto walls (or, more precisely, where Dirichlet BCs are used for velocity) for p, it only affects the regularity of the field p. Why? Because the weak formulation of the non-incremental (and incremental) pressure projection method of Chorin (Yanenko, and others) does not require a BC for pressure. The need for BCs for p comes out in some FEM/FVM/FDM because you need to solve algebraic equations derived from the mathematical statements! Since these equations (pressure poisson equation and momentum equations) are derived from one another, if you set dirichlet BCs for velocity, you need to "unconstrain" p (set it to zeroGradient), and viceversa.

On the contrary, as Re increases the zero gradient assumption becomes more similar to the momentum projection onto the BC since (1/Re) laplacian(U)*n --> 0 as Re --> \infty, hence you get more regular p fields using the natural BC for high Re wall-bounded flows. Besides, in wall-bounded turbulent flows velocity behaves linearly (in the mean) close to the wall (in the viscous region), that is: U+ = Y+, hence your laplacian is near zero anyways...

... Long story short:there is no problem with the solution running piso/pimple/ico/Foam using zeroGradient for p, where dirichlet is used for velocity...

I cannot speculate what gone wrong with your 2 attempts, so I will not comment about that.... My suggestion is to produce grad(p) during runtime either by creating a volVectorField using codeStreams, or to use postProcess to produce grad(p) directly from your solution. Like this, at least you are consistent on how the gradients were actually determined.

[christiank_ucl]

Hi Santiago,
Thank you for the information. Could you please share the references that are relevant for this?
Thanks,
Christian

[Santiago]

Originally Posted by christiank_ucl

Hi Santiago,
Thank you for the information. Could you please share the references that are relevant for this?
Thanks,
Christian

References:

LINK 1

LINK 2

LINK 3

LINK 4