Termination in parallel, but not serial calculation

[hennis722]

Hello,

I have successfully implemented the solver pimpleCentralFoam and chtPimpleCentralFoam in Openfoam version 18.12. Then I tested the pimpleCentralFoam solver using the attached forwardStep tutorial. The simulation does not crash with a serial calculation. Strangely, the simulations in parallel calculations crash due to negative temperatures. I've attached the tutorial case and the solver. It can be seen in the log file that within one time step the temperature becomes negative. Do you know what the problem is?

ForwardStep.tar.gz

[unicfdlab]

Hi, the problem is with stability of the numerical scheme -- you are using Adams-Bashforth method for derivatives (which is called "backward" in OpenFoam). This method has smaller stability region. When you run simulation in serial mode, the random disturbances to the solution are different to what you may encounter in parallel execution. Thus, the model which works in serial mode could diverge in parallel mode.

Try to play with numerical scheme settings: 1 ddtSchemes: Euler/backward 2 divSchemes: try to set defaultAdvScheme and defaultAdvSchemeV to Minmod(MinmodV)/vanAlbada(vanAlbadaV) 3 maximum Courant: try to set maxCo to 0.15 or 0.1 4 number of outer iterations: try to set nOuterCorrectors to 3 or 4 5 I always recommend to use PBiCGStab instead of PBiCG asymmetric solvers

[hennis722]

Hi,

thanks for the recommendations, the change to the matrix solver PBiCGStab and backward helped. I still have open questions about the pimpleCentralFoam and chtPimpleCentralFoam solvers because it looks very interesting for my problem.

Questions:

1. The paper under the link https://doi.org/10.1016/j.procs.2015.11.007 is for the solver pisoCentralFoam what exactly is the difference to the pimpleCentralFoam solver? Is there another paper on that? And is there a paper showing the validation of the chtPimpleFoam solver?

2. Do you have any recommendations for the solver settings (fvSolution, fvScheme, Courant number, ...) for the simulation of gas-dynamic processes to be resolved over time (many reflections of shock waves or like the popular sod problem) or has something like this been calculated using the pimpleCentralFoam / rhoPimpleCentralFoam solver? Since I want to calculate a transient gas dynamic effect in an unstructured mesh and the temperature behavior of a thermocouple at Helium Mach 5.

3. The difference between pimpleCentralFoam and rhoPimpleCentralFoam lies in the use of ideal and real gases. But wouldn't it just be possible to specify a real gas model in the thermphysicalProperties file of the pimpleCentralFoam solver?

4. Has the solver already calculated impact reflections so that a significant increase in pressure was measured on the wall or on a pressure sensor?

5. Have large eddy simulations of supersonic flows been performed with the solver and showed good results? Since the Kurganov Scheme is known for a high numerical dissipation and occurring eddies are not well resolved in comparison to RANS calculations. Maybe also in a paper?

[unicfdlab]

1. pisoCentralFoam uses PISO algorithm for pressure-velocity coupling (i.e., without outer iterations), while pimpleCentralFoam could use outer (SIMPLE) and inner (PISO) iterations. Inner iterations do not conserve total energy in compressible, high speed flows. However, I tried 2 and 3 outer iterations (SIMPLE) for some Riemann problems (see E. Toro) and it worked well for Courant numbers up to 0.9. For pimpleCentralFoam, it is better to refer to https://doi.org/10.1002/fld.4512 article.

2. I can give you only some general advice: start with the most diffusive and stable configuration (upwind, Euler, Co=0.1, 2 outer correctors) and then continue by improving the accuracy (upwind->Minmod->vanAlbada->vanLeer, Euler->backward) and the computational costs (Co=0.1->0.2->0.3). One note: the higher is Mach number, the less is the scheme stability. I would recommend to use upwind for Ma > 6, or blended schemes which can switch between upwind/vanLeer depending on local Mach number.

3. No, it will not work. The difference between EoS leads to the difference in a pressure equation.

4. Please, clarify this question.

5. See this project: https://www.researchgate.net/project/Development-and-implementation-of-hybrid-Density-Pressure-scheme-for-compressible-flows-simulation-in-OpenFOAM and also works from Prof. M. Pfitzner and his student Matthias Banholzer: they've done LES of supercritical liquids with hybrid approach. We are working on the implementation of HLLC fluxes for hybrid approach, but I don't know when it will be finished. As an alternative, you can try QGDFoam -- it is more computationally expensive, but gives better results in transient simulations due to absence of Riemann solvers and flux limiters

[hennis722]

First of all, thank you for your feedback. Do you think there is a way to increase the CFL number? At the moment I am looking at a similar problem to the one from your paper https://doi.org/10.1109/IVMEM.2018.00019 on the airbag gas generator and calculate with a CFL number of 0.2. With backward and Minmod / MinmodV, it is important that the second order is still calculated.

[unicfdlab]

You can try CFL=0.3. If it doesn't work, you can implement blended flux limiter, which blends between Minmod (or vanLeer) and upwind depending on the local (face) Mach number. And also one important note -- sometime you have to limit temperature and other variables -- see for example https://github.com/unicfdlab/realLifeExamples/blob/master/PlasmaExpansion/OpenFOAM-3.0.0/100PaOutlet/system/fvOptions . For my paper I made calculations with CFL=0.5.

[hennis722]

Okay, but isn't the Minmod or vanLeer process already a blending function to switch between 1st and 2nd order, but depending on the velocity? It is perhaps important to know that my problem has to be resolved in time and that it is about emptying a gas tank. Yes, I have limited the temperature to 1 Kelvin using the janaf table.

[unicfdlab]

1) No, Minmod and vanLeer blend between 1st and 2nd order depending on the gradient of the function (pressure, velocity, density, etc). However, since the computational stencil of OpenFOAM is compact, the implementation of Minmod and vanLeer doesn't correspond to the classical definition of these limiters. And sometimes for large Mach number it is better to use pure upwind without limiters. 2) I would recommend to try fvOption limitTemperature

[hennis722]

Okay strange, then how do I know if density, pressure or speed is selected for the gradient of the function? Have you ever programmed the limiter depending on the Mach number? Unfortunately, I have no experience with it yet.

[unicfdlab]

I mean, the interpolation of "p" to face is limited according to gradient of "p". For the blending of Minmod and upwind, please refer to localBlended

[hennis722]

I have already tested many settings in the last few days, but so far without success. I calculated a supersonic jet similar to the Ladenburg supersonic jet. The mesh and showing time are the same. The initial condition is a pressure driven flow with on the left side 150 bar and right side 1 bar.

[unicfdlab]

Try to set ddtSchemes to Euler, maxCo to 0.2

[unicfdlab]

And set defaultAdvScheme to Minmod

[hennis722]

With Euler, Minmod and maxCo 0.2 the Mach number looks more physical. However, the deviations from rhoCentralFoam are still large. Normally should pimpleCentralFoam only use the Kurganov Scheme because of the supersonic flow? And is a calculation with the 2nd order of the temporal term not possible?

[unicfdlab]

Generally, the 2nd order Adams-Bashforth time derivative approximation scheme (aka 'backward') is not intended for usage with Kurganov/Tadmor schemes.

I know that pimpleCentralFoam solver has problems with Mach > 5-6 on triangular girds and limited schemes. Set defaultAdvScheme(V) to upwind.

If it works, then you can switch to upwind/vanLeer scheme.

[hennis722]

Which mesh would you recommend? I use polyeders normally.

Where would you implement the localBlendend in the fvScheme?

Did you run the airbag gas generator from the paper https://doi.org/10.1109/IVMEM.2018.00019 with Euler, upwind, maxCo 0.5?

[unicfdlab]

Here are settings from Airbag case controlDict.txt fvOptions.txt fvSchemes.txt fvSolution.txt

[unicfdlab]

And I've used mesh from snappyHexMesh

[unicfdlab]

Where would you implement the localBlendend in the fvScheme?

To implement localBlended , you need small changes in the source code of pimpleCentralFoam. We can do it later, when it will be clear why your solution diverges

[hennis722]

Thank you ! The limitations in the fvOptions have helped to approach the rhoCentralFoam solver. I used backward, minmod, celllimited Gauss linear 1 and maxCo 0.45.

If I want to use the Peng Robinson equation in the rhoPimpleCentralFoam solver, is only an adjustment necessary in the thermophysicalProperties Dict from "covolumeGas" to "PengRobinson"? Or also in the "additionalRhoThermos.C" file you defined in the pimple library in the function "makeThermo(...)"? additionalRhoThermos.C.txt

[unicfdlab]

Thank you ! The limitations in the fvOptions have helped to approach the rhoCentralFoam solver. I used backward, minmod, celllimited Gauss linear 1 and maxCo 0.45.

I've recently remembered, that usage of non-vector version of limiter (e.g., Minmod instead of MinmodV) could also improve the situation

[unicfdlab]

If I want to use the Peng Robinson equation in the rhoPimpleCentralFoam solver, is only an adjustment necessary in the thermophysicalProperties Dict from "covolumeGas" to "PengRobinson"? Or also in the "additionalRhoThermos.C" file you defined in the pimple library in the function "makeThermo(...)"? additionalRhoThermos.C.txt

If the combination of properties you specified in thermophysicalProperties doesn't exist, then you have to create a new one in additionalRhoThermos.C and then you must recompile the library (compressibleTools?)

[hennis722]

I used the pisoCentral library for this, but I just saw that libcompressibleTools also contains the file.

I think I have one last question: Is there a paper for the solver reactingPimpleCentralFoam in which it has already been validated? I will need it for the mixing of two gases without reaction.

[unicfdlab]

I used the pisoCentral library for this, but I just saw that libcompressibleTools also contains the file.

pisoCentral is also OK. I forgot, that it also contains such classes. However, libcompressibleTools contains very usefule BC when you can have both sub- and supersonic flow at an outlet

[unicfdlab]

I used the pisoCentral library for this, but I just saw that libcompressibleTools also contains the file.

pisoCentral is also OK. I forget, that it also contains such classes. However, libcompressibleTools contains very usefule BC when you can have both sub- and supersonic flow at an outlet

I made some comparisons, but I don't have a paper on this solver. On the other hand, many people used it for their research successfully, check Scopus database for citations of main articles.

You can also refer to this PhD thesis:

B. E. Schmidt, On the stability of supersonic boundary layers with injection, Ph.D. thesis, California Institute of Technology (2016).

[unicfdlab]

I found good article for the verifcation:

J.R.Edwards, M. Ling, Low-Diffusion Flux-Splitting Methods for Flows at All Speeds.AIAA Journal 1998

[hennis722]

Okay thanks, I'll take a look at the paper. I am currently investigating a conjugate heat transfer simulation in the supersonic / hypersonic range with the k - Omega SST turbulence model, my mesh has a y+ value <1. Do you have any recommendations from your experience which wall function I should use there for alphat, k, nut and omega?

[unicfdlab]

Our solver QGDFoam could work for Mach numbers 100 (and maybe more) without limiters. Of course, there are other issues with this algorithm, but in some cases it could show advantages