Originally Posted by FMDenaro (Post 852763)

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.

Originally Posted by Winston Virtaus (Post 852771)

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.

Originally Posted by FMDenaro (Post 852772)

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.

Originally Posted by FMDenaro (Post 852772)

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.

Originally Posted by Winston Virtaus (Post 855851)

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?

Originally Posted by FMDenaro (Post 855864)

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.

Originally Posted by Winston Virtaus (Post 855903)

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.

Originally Posted by FMDenaro (Post 855907)

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.

Originally Posted by Winston Virtaus (Post 855935)

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? :rolleyes:

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.

Originally Posted by FMDenaro (Post 855937)

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 ?

Originally Posted by Winston Virtaus (Post 855939)

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.

Originally Posted by FMDenaro (Post 855944)

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.

Originally Posted by Winston Virtaus (Post 856026)

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.

Originally Posted by Winston Virtaus (Post 856032)

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?

Originally Posted by FMDenaro (Post 856033)

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.

Originally Posted by Dr Youssef Hafez (Post 856034)

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.