[bastil]

Hello all,

I have a case that is running fine in chtMultiRegionSimpleFoam in steady-state. Now I am trying to run the same case (with steady-state BCs) in transient using chtMultiRegionFoam (I will switch to transient BCs later on). So far this is a desaster. What I did:

1. Use the case without any changes but run "chtMultiRegionFoam" instead of "chtMultiRegionSimpleFoam". This gives me mass conservation for time 1. However, as expected courant numbers are way to high for PISO and the simulation is crashing after a few time steps.

2. Next step was to reduce time-stepping to 1e-4. This gives me max. courant of 0.43. However, this calculation gives significant mass inconsistency between my in- and outlet and massively to high pressure drop after the first time step. It ends up with unphysical high velocities after a few time-steps. I realized that reducing time step further (1e-6) makes everything even worse which I don't really understand. However, larger time steps (1e-2) gives better but not perfect behavior.

Has anybody observed similar effects? I have no explanation for this nor an idea what can be done to overcome this to get a stable solution.

Regards bastil

[Bloerb]

Please post your fvSchemes and fvSolution files.

[bastil]

See attached for both fluid (it's a liquid) and solid.

Inbetween I increased nOuterCorrectors from 1 to 10 to switch from PISO to pimple.

This improves the continuity error (seems to improve due to available under-relaxation) but gives me divergence typically when solving velocity.

[Bloerb]

While your files are all over the place I can't spot anything obviously wrong in them. There are things in there like a momentumPredictor in a solid, which is not read since there is no flow being solved, but nothing that makes a crash inevitable. Your schemes are fine. You might want to try changing your div(phi,U) scheme to upwind. The celllimited version you are using is fine though. Since you did not change anything in 0 and constant as I understood it is only the start up phase.

How long did it take your steady state solution to converge to reasonable levels? Did you try initializing your transient simulation with the steady state one? When and how does it crash? Can you post a log file of your run?

[bastil]

Yes, it is the start-up phase that is causing issues.

The steady-state solution took about 2000 Iters to converge. However, it was only the solid temperature that was quite slow. Flow field convergence was notably quicker.

I can post a log-file but this will take some time, I am out for a few days now. Mostly the crash happens at the beginning of a new timestep before solving velocities. I have a quite verbose log-file that gives me min and max field values for pressures, temperatures and velocities. After increasing nOuterCorrectors to 10 there is nothing obviously wrong about them before the crash occurs. Before that, mass conservation was not ok.

I tried initializing from steady-state without any success. I suspect this has something to do with different pressure levels: I use pressure "0" at the outlet for steady-state which I changed to "1e5" for transient. From my understanding the incompressible solver is pressure-normalized ("0") whereas the transient one requires absolute pressures.

[Bloerb]

There is likely your error. You forgot to change the pressure in both files, p and p_rgh to 1e5 in the initial field. You could use the zero in the transient solver as well as in the steady state one 1e5 though. It depends on your settings in thermoPhysicalProperties. If you set rhoConst there you are not solving for density and can hence use the normalized 0 pressure.

[bastil]

No I did that. I have p and p_rgh "1e5" for the transient (chtMultiRegionFoam) case. I had "0" for steady-state (chtMultiRegionSimpleFoam). And that is - I guess - why an initialization is causing more troubles than starting from scratch.

[bastil]

Attached is a log-File for a "pimple-like" run. I am currently re-running using 1e-4 timestep and piso and will post this log as well.

[Bloerb]

The obvious thing in your log file is the drastic temperature drop from iteration 1 onwards. A common error with these solvers. Usually due to poor turbulence modelling, here i am not sure though. You can use automatic time step selection in your controlDict.

maxCo           0.1;
maxDi           1e-3;
adjustTimeStep  true;

The fishy thing from my perspective about your log files is, that the momentum Equation is not solved. Try enabling momentumPredictor in your PIMPLE subdict of fvSolution. Since no turbulence model is chosen I presume you disabled it due to the low reynoldsnumber.

[bastil]

maxC0 0.1 only makes sence with PISO-like settings otherwise it may take years to run. However, it is also my first aim to get it running PISO-like (nOuterCorrectors=1).
momentumPredictor is off as this is the case in the tutorial es well. I am not really sure what it does but will turn it on.

[Tobi]

Hi guys,

just a few hints:

• SIMPLE vs PIMPLE is different and even in PISO
• That means, the mesh can make trouble or the boundary too (even if it is working for SIMPLE, it can be a problem - I guess not but you never know) - but based on your extreme temperature increase/drop after the first iteration I would expect some problem in a cell.
• You can use PIMPLE in a correct way with under-relaxation etc. Just increasing the nOuterCorrections (for non-linear and inter-coupled terms) without under-relaxing is meaningless - you will just go faster to a diverged solution (more information about that can be checked out in my book)
• momentumPredictor will just turn on or off the estimation of the velocity field; in a lot of cases this makes trouble and is not really needed if you take care about the numerical solution procedure.

Good luck.

[bastil]

Originally Posted by Tobi

SIMPLE vs PIMPLE is different and even in PISO
You can use PIMPLE in a correct way with under-relaxation etc. Just increasing the nOuterCorrections (for non-linear and inter-coupled terms) without under-relaxing is meaningless - you will just go faster to a diverged solution (more information about that can be checked out in my book)

I know. Of course I have under-relaxation turned on when I have nOuterCorrectors > 1.

momentumPredictor will just turn on or off the estimation of the velocity field; in a lot of cases this makes trouble and is not really needed if you take care about the numerical solution procedure.

Ok. I also did not realize a large difference with or without momentumPredictor.

Attached is a log-file where I did the following changes compared to the previous one:

• time step: 1e-4 instead of 1e-2
• nOuterCorrectors = 1 instead of 10

I have much more experience with PISO out of my history from other codes and therefore prefer to get it done with PISO first. As it can be seen in the log the temperature range is fine. However, velocity is getting far to high quickly. I am expecting velocities below 5 m/s. I looked into the location where this happens and it is in a relatively small hole inside my geometry (mesh there looks ok) where far to much flow is passing through at the beginning. Also it can be seen that mass continuity is not fulfilled at the beginning(inflow_wasser vs. outflow_wasser should be the same). What can be done to get a better velocity field solution?

Edit: so far I tried the following without success:
1. Reduce time-step from 1e-4 to 1e-5: Same behavior, calculation runs 5e-5s instead of 5e-4s
2. Increase nOutercorrectors from 1 to 2 to make use of under-relaxation (factor 0.5): Helps a little but not to much.

[Tobi]

Hi Bastil,

its a bit difficult to get any insight to your case and the physics behind that. The region named - kopf_al (I think its some aluminum block - engine?) and your water has a very low temperature. I guess it is fine right? The error you get is based on the Temperature evaluateion - recalculate T from h. This is based on the fact that some cell makes trouble. You can see it that your Co explodes to > 500 and thus, I think the velocity gets extreme high, the pressure drops and the evaluation in the thermodynamic model fails.

What you could do:

• Shift to a newer FOAM version (the 2.3.x is already old)
• Checking the mesh
• Reducing the time step to 0.00005 and save each time step to enable that the abort does not happen in one iteration but more like diverging within 2 - 6 time steps. If this works, you will see which cell or location will cause the problem
• PIMPLE with two outerCorrectors does not make too much sense because the first iteration will use the under-relaxing factor and the last one will use 1 (only if you explicitly tell FOAM you want to under-relax the last one too).
• PIMPLE should be set-up with residual control, in order to leave the pimple loop after your convergence criteria are fulfilled; however I know that you cannot converge always but this is case depended. That's why I set-up my nOuterCorrector to 100 and normally leave the loop after 6-7 outer corrs.
  
  But you should ensure that the solid regions should converge fast (very small relaxation for h like 0.95) plus nNonOrthoCorrections (for converging non-linear terms a bit faster) ...
• The last point could be based on the numerical schemes. I would first try to solve it with upwind. If this fails, then there should be definitely a mesh problem.

Just some hints but as always, I am not an expert and its just my experience. Maybe you can extract some helpful information out of it.

[bastil]

Originally Posted by Tobi

• Shift to a newer FOAM version (the 2.3.x is already old)

I am currently trying to complie 4.x

• Checking the mesh

Create time

Create polyMesh for time = 0
Time = 0
Mesh stats
points: 511755
faces: 4343955
internal faces: 4106445
cells: 2027153
faces per cell: 4.1686
boundary patches: 9
point zones: 0
face zones: 3
cell zones: 3
Overall number of cells of each type:
hexahedra: 0
prisms: 341769
wedges: 0
pyramids: 19
tet wedges: 0
tetrahedra: 1685365
polyhedra: 0
Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology 
...
Checking geometry...
Overall domain bounding box (-0.2795 -0.159141 -0.264175) (0.261 0.238676 0.099)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (1.43513e-17 5.98081e-18 -2.31787e-16) OK.
Max cell openness = 6.67657e-16 OK.
Max aspect ratio = 43.2037 OK.
Minimum face area = 2.22653e-08. Maximum face area = 9.52577e-05. Face area magnitudes OK.
Min volume = 2.47071e-12. Max volume = 2.8834e-07. Total volume = 0.0176281. Cell volumes OK.
Mesh non-orthogonality Max: 86.7125 average: 20.3952
*Number of severely non-orthogonal (> 70 degrees) faces: 2309.
Non-orthogonality check OK.
<<Writing 2309 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
***Max skewness = 4.6416, 1 highly skew faces detected which may impair the quality of the results
<<Writing 1 skew faces to set skewFaces
Coupled point location match (average 0) OK.
Failed 1 mesh checks.
End

• Reducing the time step to 0.00005 and save each time step to enable that the abort does not happen in one iteration but more like diverging within 2 - 6 time steps. If this works, you will see which cell or location will cause the problem

Already done. I know where the location is. It is a relatively small hole. I have several of theses holes where the flow needs to distribute through but nearly the whole mass is passing through this single hole. Therefore high velocities in there ... No idea how to overcome this, seems to be cause by the insufficient velocity solution.

(only if you explicitly tell FOAM you want to under-relax the last one too).

How?

• The last point could be based on the numerical schemes. I would first try to solve it with upwind. If this fails, then there should be definitely a mesh problem.

Also not running with upwind. However, mesh in the region with high velocity looks fine for me.

[Bloerb]

I thought checking the mesh was to obvious of an answer. I am proven wrong once again. A nonOrthogonality of 90 is invalid. Getting a stable solution with a value of 87 is extremely challenging. checkMesh -allTopology -allGeometry will likely fail additional checks. I'd highly recommend creating a mesh that it is worth computing on. Because this is definitely a mesh problem. Changing the version won't do you any good in this case.

Otherwise you might want to increase nOuterCorr to something big like Tobias mentioned and use the same relaxation as in your simple solution. You can afterwards lower it. ResidualControl for automatically ending this loop however was only recently added by him.
As he stated PIMPLE uses a Final flag in fvSolution. So relaxationFactors for the last loop would be hFinal for example.

[bastil]

Ok. I will try to improve the mesh. What are the recommended limited for a transient chtMultiRegionFoam case?