[steampunc]

Thanks in advance for the potential help. I'm relatively new to OpenFoam, so I'm really sorry if this is a stupid mistake. I can't seem to find the same problem online, but I've been searching for a couple days.

I've been trying to apply some diagnostic/validation tests to incompressible pimpleFoam configurations and have been finding some anomalous results.

For any incompressible fluid, div(U) should be 0 or very near it. From what I understand, the pimpleFoam solver uses the div(U) -> 0 condition as it solves the pressure field, and it certainly seems to be doing that, because in each step the time continuity errors are very small (order 10^-17 or so). However, when I use funkySetFields/postProcess to compute the div(U) field at a given logged time, the values are much larger (on order 1 or even larger). I've tried it on my own configurations, along with several of the pimpleFoam/pisoFoam tutorials, and all of them result in large div(U) fields.

To compute the field, I've been using these commands:

 funkySetFields -time 0: -create -field divU -expression "div(U)"

and

postProcess -func "div(U)"

From what I understand, postProcessing (and funkySetFields?) uses the same divergence schemes as OpenFoam itself, so I'm confused why it's so large. I'm almost completely certain it's due to user error, but it seems like such a simple test that there's not much room for it to go wrong.

I've been using a Gaussian divergence scheme with linear interpolation in system/fvSchemes for div(U).

Attached are some images of the div(U) field in different pimpleFoam tutorial cases in paraView.









Also, here are some excerpts of logs from the pimpleFoam runner on channel395:

PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.0117778, Final residual = 5.81593e-06, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0603226, Final residual = 2.72548e-06, No Iterations 3
smoothSolver:  Solving for Uz, Initial residual = 0.0585004, Final residual = 2.7339e-06, No Iterations 3
Pressure gradient source: uncorrected Ubar = 0.1335, pressure gradient = 2.76228e-05
GAMG:  Solving for p, Initial residual = 0.146204, Final residual = 0.0126411, No Iterations 1
time step continuity errors : sum local = 1.14563e-05, global = -8.27789e-20, cumulative = -2.0927e-18
Pressure gradient source: uncorrected Ubar = 0.1335, pressure gradient = 2.76872e-05
GAMG:  Solving for p, Initial residual = 0.0163918, Final residual = 4.2532e-07, No Iterations 5
time step continuity errors : sum local = 4.5834e-10, global = -1.12784e-19, cumulative = -2.20548e-18
Pressure gradient source: uncorrected Ubar = 0.1335, pressure gradient = 2.76392e-05
ExecutionTime = 14.27 s  ClockTime = 14 s

fieldAverage fieldAverage1 write:
    Calculating averages

Courant Number mean: 0.30284 max: 0.540451
Time = 5.4

PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.0118056, Final residual = 5.80108e-06, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0603712, Final residual = 2.65463e-06, No Iterations 3
smoothSolver:  Solving for Uz, Initial residual = 0.0584964, Final residual = 2.67235e-06, No Iterations 3
Pressure gradient source: uncorrected Ubar = 0.1335, pressure gradient = 2.76884e-05
GAMG:  Solving for p, Initial residual = 0.146366, Final residual = 0.0127159, No Iterations 1
time step continuity errors : sum local = 1.15734e-05, global = -1.22909e-19, cumulative = -2.32839e-18
Pressure gradient source: uncorrected Ubar = 0.1335, pressure gradient = 2.77468e-05
GAMG:  Solving for p, Initial residual = 0.0165409, Final residual = 5.45694e-07, No Iterations 5
time step continuity errors : sum local = 5.65104e-10, global = -1.5652e-20, cumulative = -2.34404e-18
Pressure gradient source: uncorrected Ubar = 0.1335, pressure gradient = 2.76982e-05
ExecutionTime = 14.61 s  ClockTime = 15 s

fieldAverage fieldAverage1 write:
    Calculating averages

Courant Number mean: 0.302862 max: 0.539102
Time = 5.6

[HPE]

- Do you execute `div` FO on converged cases? The reason of the question is that the divergence of velocity is expected to reduce to zero, but the intermediate steps might not well resort to zero-divergence. For `steady-state` cases, `bounded` keyword can be added to the divergence scheme, please see example.
- Which versions do you use? Have you tested different versions, and obtained similar results?

[steampunc]

Thanks for the response!

Do you execute `div` FO on converged cases?

I'm pretty sure that the cases are convergent - based off of the output of pimpleFoam at each timestep, the Courant number isn't enormous, and the residual values stay relatively small as well. Is there a more rigorous way to check for convergence? I've tried adjusting the timesteps smaller, but that didn't seem to affect my calculations of divergence on the logged steps.

The reason of the question is that the divergence of velocity is expected to reduce to zero, but the intermediate steps might not well resort to zero-divergence.

What do you mean when you say that the intermediate steps might not be zero-divergence? Are you talking about the smaller steps within a given simulation timestep (i.e. computing U, p, etc)? If so, I don't know if that potential lack of zero-divergence in an intermediate step should work its way into the final timestep logs. My understanding is that at the end of any timestep, the zero-div condition should be met, but I very well could be wrong.

For `steady-state` cases, `bounded` keyword can be added to the divergence scheme,

Unfortunately these cases are not steady-state - they involve turbulence and pimpleFoam solving, so I can't use the bounded schemes which would maybe give a better result.

Which versions do you use? Have you tested different versions, and obtained similar results?

I'm doing this on OpenFoam v1906, but I haven't tried it on other versions of OpenFoam. I can try that and see what happens.

Again, thank you so much for the response! I'm really hoping we can figure this out!

[HPE]

>> My understanding is that at the end of any timestep, the zero-div condition should be met, but I very well could be wrong.

I think that might be the subtlety.

To my understanding, div(U)=0 can be reached for an incompressible flow computation at the end of a converged simulation.

Quoting from this link:

"For numerical solution of incompressible flows, div(U)=0 at convergence, at which point the third term on the right hand side is zero. Before convergence is reached, however, div(U)=0 and in some circumstances it is better to include the term within a numerical solution because it helps maintain boundedness of the solution variable and promotes better convergence."

So, I assume you measure "div(U)=0" at the end of each time-step, rather than at the end of a simulation?

[steampunc]

Sorry for not responding for a while.

I might be mixing terms, but yes - I mean "div(U)=0" at the end of each timestep, which doesn't seem to be the case when I compute it. This is a dynamic simulation, which I'm not expecting to converge to a steady state.



It's just very surprising that the divergence is nonzero in this incompressible solver. I can't wrap my head around it!

[HPE]

just numerics/floating point representation of the information - every CFD solver is prone to it. I wouldn't worry about it.