[Winston Virtaus]

Originally Posted by FMDenaro

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.

Yes, I have confidence in the method now. I was planning to apply this later to external flows also and see how well it works e.g. in nozzle type flows where a jet exits into theoretically zero pressure environment.

[FMDenaro]

Originally Posted by Winston Virtaus

Yes, I have confidence in the method now. I was planning to apply this later to external flows also and see how well it works e.g. in nozzle type flows where a jet exits into theoretically zero pressure environment.

Just remember to consider always the Hodge decomposition. On an outflow, if you apply a normal derivative of the velocity to be zero, then you have a relation for the second derivative of the pressure along the normal direction to use in the Poisson equation. If you do all correctly, on a staggered grid you will be able to drive the div V = 0 constraint to machine precision.

[Winston Virtaus]

Originally Posted by FMDenaro

Just remember to consider always the Hodge decomposition. On an outflow, if you apply a normal derivative of the velocity to be zero, then you have a relation for the second derivative of the pressure along the normal direction to use in the Poisson equation. If you do all correctly, on a staggered grid you will be able to drive the div V = 0 constraint to machine precision.

The decomposition seems indeed to be the preferred method to think about valid boundary conditions, I'll keep this in mind!

Usually I have been using the standard do-nothing boundary condition with zero normal velocity gradient and zero pressure at outlet. Previously I have also tried the traction-free boundary condition where pressure is set by the normal velocity gradient but this didn't provide much benefit as an outlet condition over the standard boundary condition. I didn't pay much attention if these were consistent with the decomposition.

[Winston Virtaus]

I was able to some more testing on the Neumann pressure boundary conditions on cell centered grids. Hopefully soon proceeding to more complex 3D cases from the backwards facing step case. No scaling tests yet perfomed though.

Originally Posted by FMDenaro

Just remember to consider always the Hodge decomposition. On an outflow, if you apply a normal derivative of the velocity to be zero, then you have a relation for the second derivative of the pressure along the normal direction to use in the Poisson equation.

It was mentioned a while back that a zero gradient condition for velocity would lead to a zero second gradient of pressure via the Hodge decomposition. Is it a bad idea to assume a constant
pressure gradient across the boundary cell even if there is a nonzero velocity gradient present e.g. in wall boundaries?

I mean that if intermediate velocity gradients at boundary are set to be equal to physical velocity gradients before intermediate velocity field is computed (not implying same boundary values, but same gradients), wouldnt zero second gradient of pressure be a direct consequence of the Hodge decomposition?

[FMDenaro]

Originally Posted by Winston Virtaus

I was able to some more testing on the Neumann pressure boundary conditions on cell centered grids. Hopefully soon proceeding to more complex 3D cases from the backwards facing step case. No scaling tests yet perfomed though.



It was mentioned a while back that a zero gradient condition for velocity would lead to a zero second gradient of pressure via the Hodge decomposition. Is it a bad idea to assume a constant
pressure gradient across the boundary cell even if there is a nonzero velocity gradient present e.g. in wall boundaries?

I mean that if intermediate velocity gradients at boundary are set to be equal to physical velocity gradients before intermediate velocity field is computed (not implying same boundary values, but same gradients), wouldnt zero second gradient of pressure be a direct consequence of the Hodge decomposition?

The Hodge decomposition is a relation between continuous vector fields, therefore you can apply any derivative on it.

From a rigorous mathematical aspect, that is no longer true on confinement boundary where derivative is defined only from the interior.
However, let assume any derivative can be done also on the boundary, the issue would be when you prescribe some approximation for the derivative.

Thus, yes, if you assume the velocities derivative are the same the second derivative of the pressure will vanish on the boundary. A different issue is what that means in the physics of that region.

[Winston Virtaus]

Originally Posted by FMDenaro

The Hodge decomposition is a relation between continuous vector fields, therefore you can apply any derivative on it.

From a rigorous mathematical aspect, that is no longer true on confinement boundary where derivative is defined only from the interior.
However, let assume any derivative can be done also on the boundary, the issue would be when you prescribe some approximation for the derivative.

Thus, yes, if you assume the velocities derivative are the same the second derivative of the pressure will vanish on the boundary. A different issue is what that means in the physics of that region.

Its good to hear zero second derivative of pressure has some validity behind it. Its a convenient option to use numerically so i decided to go with that route.

Currently I'm using v* = v(n+1) +deltaT*gradP(n) for intermediate velocity step and the dP/dn = (v*-v).n at the Poisson step.
Are there any rules of thumb to make an assessment if the resulting v is actually divergence free enough after the projection step? Average Div(v) is approximately O(1e-4) at Co=1 in the backwards facing step case.

[FMDenaro]

Originally Posted by Winston Virtaus

Its good to hear zero second derivative of pressure has some validity behind it. Its a convenient option to use numerically so i decided to go with that route.

Currently I'm using v* = v(n+1) +deltaT*gradP(n) for intermediate velocity step and the dP/dn = (v*-v).n at the Poisson step.
Are there any rules of thumb to make an assessment if the resulting v is actually divergence free enough after the projection step? Average Div(v) is approximately O(1e-4) at Co=1 in the backwards facing step case.

In a staggered arrangement of the variables you must have in each cell the divergence-free constraint verified at the same error of the stop there should in an iterative solver. That can be also zero at the roundoff level. On non-staggered arrangement, the constraint is satisfied up to the magnitude of the local truncation error. At least in a general case, some specific method can produce the same error of the staggered case.
Check the max error of the div V condition and the position.

[Winston Virtaus]

Originally Posted by FMDenaro

In a staggered arrangement of the variables you must have in each cell the divergence-free constraint verified at the same error of the stop there should in an iterative solver. That can be also zero at the roundoff level. On non-staggered arrangement, the constraint is satisfied up to the magnitude of the local truncation error. At least in a general case, some specific method can produce the same error of the staggered case.
Check the max error of the div V condition and the position.

The largest divergence seems to be at the boundary cells close to the 90 degree corner where flow separates due to abrupt expansion of the channel height. There div(V) is around O(100) which is many orders of magnitude higher than the average would suggest. Both positive and negative values of divergence occur in the boundary and the high divergence is limited only to a couple of cells. In the internal cells I can see some divergence errors O(10) occuring further away from the corner at the shearing layer that is formed after the backward step. The errors approach the average levels O(1e-3 to 1e-4) fairly quickly in the streamwise direction.

The only physical interpretations of div(V) I currently have is that its related to the rate of change of volume per unit volume and its related to fluid normal stresses via the Newtonian fluid stress tensor but nothing conclusive yet.
So this in a way leads to the question that are the divergence numbers indicative of the solution quality on itself or should they be evaluated in relation to some other physical quantity? Perhaps get div(V) so low that it would constitute to less than 0.1% of the normal stresses?

The face fluxes have been in general adjustable via solver tolerance but div(V) at cell center has been harder to get down to arbitrary levels which is probably due to the inevitable local truncation error like you pointed out.

[FMDenaro]

Originally Posted by Winston Virtaus

The largest divergence seems to be at the boundary cells close to the 90 degree corner where flow separates due to abrupt expansion of the channel height. There div(V) is around O(100) which is many orders of magnitude higher than the average would suggest. Both positive and negative values of divergence occur in the boundary and the high divergence is limited only to a couple of cells. In the internal cells I can see some divergence errors O(10) occuring further away from the corner at the shearing layer that is formed after the backward step. The errors approach the average levels O(1e-3 to 1e-4) fairly quickly in the streamwise direction.

The only physical interpretations of div(V) I currently have is that its related to the rate of change of volume per unit volume and its related to fluid normal stresses via the Newtonian fluid stress tensor but nothing conclusive yet.
So this in a way leads to the question that are the divergence numbers indicative of the solution quality on itself or should they be evaluated in relation to some other physical quantity? Perhaps get div(V) so low that it would constitute to less than 0.1% of the normal stresses?

The face fluxes have been in general adjustable via solver tolerance but div(V) at cell center has been harder to get down to arbitrary levels which is probably due to the inevitable local truncation error like you pointed out.

But what about the grid arrangement ? The error seems due to some bug.
Have you used the same discrete operators for Div V and Div Grad p ?

[Winston Virtaus]

Originally Posted by FMDenaro

But what about the grid arrangement ? The error seems due to some bug.
Have you used the same discrete operators for Div V and Div Grad p ?

Block-structured grid with co-located arrangement.

Linear interpolation is used for the source term and compact laplacian for div(grad(p)). The schemes should be consistent.

I'll try to double check everything.

Its possible my post processing step has the bug. I'll try to post process the source term of the Poisson equation.

[FMDenaro]

Originally Posted by Winston Virtaus

Block-structured grid with co-located arrangement.

Linear interpolation is used for the source term and compact laplacian for div(grad(p)). The schemes should be consistent.

I'll try to double check everything.

Its possible my post processing step has the bug. I'll try to post process the source term of the Poisson equation.

Let me remark that a projection method on co-located grid and linear flux reconstruction produces a div v error of the order of the local truncation error in case of the compact DivGrad operator (Approximate Projection Method), as I've shown here (see Eq.(27)):


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


and also here



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





Now, if your solution is performed using non-dimensional variables, the high error in the divergence-free constraint is due to some error in prescribing the BCs.


I strongly suggest to use a standard test-case, single block structured, such as a simple Poiseuille solution to test your code.

[Winston Virtaus]

Originally Posted by FMDenaro

Let me remark that a projection method on co-located grid and linear flux reconstruction produces a div v error of the order of the local truncation error in case of the compact DivGrad operator (Approximate Projection Method), as I've shown here (see Eq.(27)):


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


and also here



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





Now, if your solution is performed using non-dimensional variables, the high error in the divergence-free constraint is due to some error in prescribing the BCs.


I strongly suggest to use a standard test-case, single block structured, such as a simple Poiseuille solution to test your code.

Thanks, I'll take a look!

Just to clarify, there are three options that I can think of currently to estimate the divergence at postprocessing. All options have 1/s unit if normalised via cell volume or m3/s if not normalised.

Option 1: Interpolate cell-centre velocity to cell faces, perform dot product with face area vector and calculate div(V) using 1/Volume*sum(Vf*Sf).

Option 2: Calculate div(V) as trace(grad(V)). Values are usually equal to option 1 within some roundoff error and slightly dependent on grad(V) discretisation choice.

Option 3: Calculate div(V) directly from the face fluxes which are decomposed as Vf = Vf* - deltaT*dP/dn and perform dot product with face area vector and calculate div(V) using 1/Volume*sum(Vf*Sf).

The values that were reported earlier were obtained from option 1 (unit 1/s) so it seems that they were per unit volume. I referred to the values given by options 1-2 as divergence of cell centre velocity and option 3 as divergence of fluxes. Which was probably abit inaccurate. The level of divergence of the decomposed face fluxes can usually be reduced via tightening Poisson solver tolerance.

If option 1 is used without the cell volume normalisation (unit m3/s), the largest divergence is in the same area as before but the values are now substantially lower. Roughly O(1e-8) at flow separation point and O(1e-9) at the shear layer. Which correspond to about 0.01% and 0.001% of flowrate through the domain. The flowrates are of course dependent on the width of the domain now that my case is actually 2D but it shouldn't affect the percentages.

If option 3 is used, the divergence is around O(1e-12-1e-13) and its scattered fairly well throughout the domain without any visible problem areas.

Should I be mainly looking at the values obtained from option 1 in m3/s units?

[FMDenaro]

You have no option! You must use exactly the same discrete operators you used to solve the pressure equation, no other choice is correct.
See the papers I addressed.

[Winston Virtaus]

The first paper hopefully clarified some concepts.

So to estimate the divergence we can take a look at two different sets of face velocities:

Set 1: "MAC-like" face velocities that live in the cell faces and are obtained through Vf =V* - deltaT*dP/dn, where V* is an interpolated quantity. Divergence should be obtained down to machine precision.

Set 2: Obtained from linear interpolation of cell-centre variables. Vf = V +deltaT*grad(P). Divergence of this quantity should be monitored due to its sensitivity to truncation errors.

Am I back on the right track? It looks like this is exactly the issue that Rhie-Chow interpolation is trying to reduce. The discrepancy between large and compact stencils.

[FMDenaro]

Originally Posted by Winston Virtaus

The first paper hopefully clarified some concepts.

So to estimate the divergence we can take a look at two different sets of face velocities:

Set 1: "MAC-like" face velocities that live in the cell faces and are obtained through Vf =V* - deltaT*dP/dn, where V* is an interpolated quantity. Divergence should be obtained down to machine precision.

Set 2: Obtained from linear interpolation of cell-centre variables. Vf = V +deltaT*grad(P). Divergence of this quantity should be monitored due to its sensitivity to truncation errors.

Am I back on the right track? It looks like this is exactly the issue that Rhie-Chow interpolation is trying to reduce. The discrepancy between large and compact stencils.

The MAC strategy is effective on staggered grid where you can get machine precision, but you have the APM on your colocated arrangement.
You have to check the divergence-free constraint in each cell, but using the congruent discretization for the div operator. The result should be an error of the magnitude of the local truncation error.
Of course, also the grad operator must be congruent.

[Winston Virtaus]

Thanks for the input, much appreciated!

Actually I wasn't aware that the mac-like face velocities werent as effective on co-located grids as I had previously thought.

Should the APM face velocity set Vf = V +deltaT*grad(P) be then used at correction step too and not just in postprocessing?

[FMDenaro]

Originally Posted by Winston Virtaus

Thanks for the input, much appreciated!

Actually I wasn't aware that the mac-like face velocities werent as effective on co-located grids as I had previously thought.

Should the APM face velocity set Vf = V +deltaT*grad(P) be then used at correction step too and not just in postprocessing?

The correction step must indeed be performed using the same discretization. In my paper you will find the full procedure for the APM (do not consider the DPM part). Just observe that the check for the divergence-free error is a part of the solver.

[Dr Youssef Hafez]

Sorry I have not followed all the comments here but just wanted to share my view. There is an alternative formulation for the pressure in which no need to use the Poisson's equation for the pressure. In this way the pressure is solved on the element or cell level to enforce the continuity equation which is not needed thus reducing the numebr of equations to only the momentum equations as follows.
The pressure has a mathematical, not a thermodynamic meaning. According to RANS the pressure is not uniquely defined as any pressure function could satisfy the RANS up to an additive constant so a family of parallel pressure surface fields would satisfy the RANS mathematically speaking. However, knowing one pressure value would tie the solution. In fact, from a CFD or mathematical prospective the pressure is just a constraint to enforce the incompressibility constraint or the velocity divergence free condition (Div. V =0).

One possible treatment of the pressure is through expressing the pressure using the penalty approach as:
- p/k = (Div. V) ; or p = - K (Div. V)
where K is a penalty parameter that is supposed to have a very large values (ideally infinite value) so at such very large values Div. (V) =0.0 .
A better formulation is the iterative penalty method which I used in my codes which is:
p (at iteration i+1) = p (at iteration i) - K (Div. V)
In the iterative penalty method smaller values could be taken for k, e.g. k maybe 100 to avoid using very large numbers.
When using the pressure penalty approach, no need to solve the continuity equation on the global scale and in fact the pressure is substituted for using the above mentioned equation in RANS and the pressure is no longer an unknown which reduces the number of unknowns to only the velocities. This is a big computational advantage saving significant memory storage. However, the pressure can be recovered by solving on the element or cell scale using the penalty equation and the obtained values are true within an additive constant. Knowing one value of pressure, the pressure values obtained could be easily scaled.

[Winston Virtaus]

Originally Posted by FMDenaro

The correction step must indeed be performed using the same discretization. In my paper you will find the full procedure for the APM (do not consider the DPM part). Just observe that the check for the divergence-free error is a part of the solver.

I was able to get hopefully some better divergence error estimate numbers by monitoring directly the divergence of the mac-like face velocities during runtime. However, I havent yet switched to the APM procedure which would be the more correct choice.

The solver was tried on a 3D channel flow @ Retau = 395 to see if it produced any different behaviour to the backward facing step case.

In the channel flow case the maximum divergence error is around O(1e-7) and its located roughly in the buffer layer region. Usually third cell from the wall.

In the backward facing step case the maximum divergence is now around O(1e-5) and its located far in the downstream where the grid is stretched in the stream-wise direction. Usually second or third cell near the wall.

Any comments on these values now? Is the incompressibility constraint satisfied atleast somewhat acceptably?

[Winston Virtaus]

Originally Posted by Dr Youssef Hafez

Sorry I have not followed all the comments here but just wanted to share my view. There is an alternative formulation for the pressure in which no need to use the Poisson's equation for the pressure. In this way the pressure is solved on the element or cell level to enforce the continuity equation which is not needed thus reducing the numebr of equations to only the momentum equations as follows.
The pressure has a mathematical, not a thermodynamic meaning. According to RANS the pressure is not uniquely defined as any pressure function could satisfy the RANS up to an additive constant so a family of parallel pressure surface fields would satisfy the RANS mathematically speaking. However, knowing one pressure value would tie the solution. In fact, from a CFD or mathematical prospective the pressure is just a constraint to enforce the incompressibility constraint or the velocity divergence free condition (Div. V =0).

One possible treatment of the pressure is through expressing the pressure using the penalty approach as:
- p/k = (Div. V) ; or p = - K (Div. V)
where K is a penalty parameter that is supposed to have a very large values (ideally infinite value) so at such very large values Div. (V) =0.0 .
A better formulation is the iterative penalty method which I used in my codes which is:
p (at iteration i+1) = p (at iteration i) - K (Div. V)
In the iterative penalty method smaller values could be taken for k, e.g. k maybe 100 to avoid using very large numbers.
When using the pressure penalty approach, no need to solve the continuity equation on the global scale and in fact the pressure is substituted for using the above mentioned equation in RANS and the pressure is no longer an unknown which reduces the number of unknowns to only the velocities. This is a big computational advantage saving significant memory storage. However, the pressure can be recovered by solving on the element or cell scale using the penalty equation and the obtained values are true within an additive constant. Knowing one value of pressure, the pressure values obtained could be easily scaled.

Interesting concept!

Is this mainly used in the FEM-community? I've usually seen penalty approach used for displacement boundary conditions and contact constraints in solid mechanics.