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.

Originally Posted by Winston Virtaus (Post 856140)

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?

Originally Posted by Winston Virtaus (Post 856255)

I'm still working on how the divergence error scales on different grids but I thought I'd ask a quick question about the local truncation error you mentioned.

In the Analysis of the local truncation error in the pressure-free projection method for incompressible flows: a new accurate expression of the intermediate boundary conditions paper its mentioned that the splitting error of the pressure free method can be estimated as
efs = deltaT*[grad(phi)-grap(p)]-deltaT^2/2Re*div(grad(grad(p)) + O(deltaT^3) (equation 24).

If first order implicit Euler is used for the diffusion term in the v* equation, does this then imply that I'm getting only the first order approximation of p i.e. grad(phi) = grad(p)? And if this is indeed the case, would the equation 24 splitting error reduce to efs = -deltaT^2/2Re*div(grad(grad(phi))). Of which the local truncation error part is LTE = efs/deltaT?

Originally Posted by Winston Virtaus (Post 856343)

Following similar analysis for the Implicit Euler as the Crank-Nicholson show in the paper, this was the result:

efs = deltaT*(grad(phi)-grad(p)) - 1/(2*Re)*deltaT^2*div(grad(grad(p))= 1/(2*Re)*deltaT^2*div(grad(grad(phi))) + O(deltaT^3)

leading to error proportional to O(deltaT^2). Is this correct now?

It seems that the Crank-Nicholson option would exactly cancel this leading term.

I hope the procedure was followed properly now :)

Originally Posted by Winston Virtaus (Post 856429)

I think it was due to some calculation error in the higher order terms.

However I could do a lower order approximation to the splitting error as



which can be further reduced into



assuming



but I'm not exactly sure if the above assumption is correct.

Then using



the resulting split error would be



Any thoughts? :D

Originally Posted by Winston Virtaus (Post 856482)

Thanks for the confirmation!

The paper is indeed very comprehensive and rigorous. The time-discrete analysis was abit harder to follow due to the inverse operator so thats why I followed the continuous analysis.

Its going to be interesting to see what kind of splitting error numbers I get when the error relation is postprocessed. The relation seems rather cheap to compute, is it customary to keep tracking it during runtime too?

Originally Posted by Winston Virtaus (Post 857390)

I did some more testing on the boundary conditions in a simpler test case to which a closed form solution is available:

Poisseule flow between infinite parallel plates at Re≈26.6

• Peak flow velocity of 1 m/s
• Pressure gradient of 1 m/s2
• Mesh resolution was 16x16x16.
• Timestep roughly 1e-3 seconds.

Boundary conditions

• Dirichlet for physical velocity
• Taylor series based BCs for intermediate velocity
• Neumann conditions for pressure

I decided to calculate the actual error from the analytical solution. The correct APM formula is also used for face velocities.

Is the divergence error now low enough compared to the velocity and pressure gradient errors?

Originally Posted by FMDenaro (Post 857391)

I suppose you have the height of 1m, the magnitude of the error seems quite congruent for a second order scheme. However, I would suggest to check further the BC at the inflow, there (and at the two corners), it appears some peaks. Since you are prescribing Dirichlet conditions both at inflow and outflow, I see a different distribution of the errors.

Originally Posted by Ogbency001 (Post 857427)

I want to achieve the following objectives and need your guide Develop a simplified CFD model to simulate the bio-pesticide
sprayer application process.

Analyze the spray pattern and droplet behaviour using the CFD
model.

Optimize the sprayer design and operational parameters to
enhance spray coverage and droplet size distribution.

Validate the CFD model by comparing the simulated spray
characteristics with available experimental data or previous
studies

Question1: Which solver in OpenFoam is best for this ?
Question 2: I have learn many tutorial and know how to do simulation based on the tutotial but can not find any on pesticide sprayer; what should I do

Originally Posted by Winston Virtaus (Post 858088)

Thanks for the input! :)

Two different tests were performed to see the divergence error scaling:

• In the first plot mesh size was fixed to 0.04 m and timestep was changed from roughly Co=1 to Co=0.01

• On the second plot timestep was fixed to Co=0.6 and mesh size was instead changed. The mesh size was reduced from about 0.1 m to about 0.01 m by splitting the mesh equally in all directions between runs.

All plots are on a double logarithmic scale.

It looks like there is still some inconsistency in the code.

Originally Posted by kinematic_presser (Post 837274)

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.