[Jabo]

Hi,

PengRobinsonGas crashes with the following error:

Program received signal SIGFPE, Arithmetic exception.
0x00007ffff34bba95 in __ieee754_pow_sse2 (x=-1.2895509926384718e-08, y=0.33333333333333331) at ../sysdeps/ieee754/dbl-64/e_pow.c:157
157     ../sysdeps/ieee754/dbl-64/e_pow.c: No such file or directory.
(gdb) bt
#0  0x00007ffff34bba95 in __ieee754_pow_sse2 (x=-1.2895509926384718e-08, y=0.33333333333333331) at ../sysdeps/ieee754/dbl-64/e_pow.c:157
#1  0x00007ffff34cc9ac in __pow (x=-1.2895509926384718e-08, y=0.33333333333333331) at w_pow.c:27
#2  0x00007ffff7aa36aa in Foam::PengRobinsonGas<Foam::specie>::Z(double, double) const () from /home/james/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/libfluidThermophysicalModels.so
#3  0x00007ffff7ab3596 in Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::PengRobinsonGas<Foam::specie>, 8>, Foam::sensibleEnthalpy>, 8> > >::calculate() () from /home/james/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/libfluidThermophysicalModels.so
#4  0x00007ffff7ab9b0a in Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::PengRobinsonGas<Foam::specie>, 8>, Foam::sensibleEnthalpy>, 8> > >::correct() () from /home/james/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/libfluidThermophysicalModels.so
#5  0x000000000041e75b in main ()

My case simulates nitrogen flow above critical pressure and temperature. Values used are 3.97 MPa and 126.9 K.

The case runs fine with a very coarse mesh, 8000 hex cells. The geometry I need to use has around 4.2 million hex cells. Both meshes pass checks when running checkMesh -allGeometry -allTopology. The case always crashes like this with the fine mesh, and I've tried cases in sprayFoam and in rhoPimpleFoam.

I've tried many different fvSolution tolerances, residual control, correctors. All crash.

My case is attached. Please, does anybody have any advice on what to try?

[dkxls]

Originally Posted by Jabo

I've tried many different fvSolution tolerances, residual control, correctors. All crash.

Some general hints:
I have had just a very quick look at your case, but what sticks out are your solver tolerances. I would generally tighten them a bit for h, U, p, rho; specifically your final tolerances should be the same or lower than the intermediate ones.

The pitzDaily tutorial gives you a pretty good starting point:

tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSolution

Just add the enthalpy to the velocity/energy parameters, i.e.:

"(U|h|e|k|nuTilda)"
"(U|h|e|k|nuTilda)Final"

Starting from there, you can then lower the solver tolerances to 1e-8 or 1e-9, this depends a lot on your case and the 1e-6 might be already enough. Note that relTol should always be 0 for the final iteration, while 0.01 is OK for the intermediate ones.

Regarding PISO vs. PIMPLE, I would probably start with the settings as given in the pitzDaily tutorial, i.e. 3 PIMPLE iterations, with each 1 PISO corrector.

Also the relaxationFactors from the pitzDaily tutorial are a good starting point, i.e. no under-relaxation for the fields and 1 for all equations.

How are the non-orthogonality values of your mesh?



Some more specific points regarding the real gas equation:
1. Have you tried to run the case with a ideal gas model? The results might be off, but it could help in debugging, especially considering that the PengRobinson implementation is rather new and I guess not yet widely tested.
2. The current implementation of the real gas equations in OpenFOAM 2.3.x doesn't consider departure functions for heat capacity, enthalpy or entropy. This will affect your results as well, but I don't think it will lead to the crashes you have.
The two real gas implementations in the "extend" version have departure functions implemented AFAIK: http://sourceforge.net/projects/open...oam-extend-3.1

[Jabo]

Hi Armin, thanks for your help.

I started with the pitzDaily tutorial and modified it for this case.

The Mesh non-orthogonality Max: 35.6935 average: 3.63311

I've tried lowering the tolerances in the way you suggest.

I think I will try using the extend version next, thanks for that suggestion.

[dkxls]

Originally Posted by Jabo

The Mesh non-orthogonality Max: 35.6935 average: 3.63311

OK, so the non-orthogonality shouldn't be a problem, though you could try to set the nNonOrthogonalCorrectors to 1 or 2. But I don't think that fixes the problem.

I also just noted that you don't solve for the momentum equation, try to run with the case with lower tolerances and the momentum predictor step, i.e.

momentumPredictor yes;

Originally Posted by Jabo

I think I will try using the extend version next, thanks for that suggestion.

I would be surprised if using the extend version fixes the problem.

[Jabo]

I tried the momentumPredictor already and nNonOrthogonalCorrectors 1.

The geometries I compared are a bit different, in that the fine one is a cylinder with an inlet and a faceplate, the coarse one is a square pipe with just inlet and outlet. I have been considering that the zeroGradient boundary conditions on the faceplate have been the issue here but so far changing them to other types hasn't work.

The crash typically happens at ~1.5e-5 seconds so maybe it's the boundary conditions.

[Jabo]

Any suggestions other than fvSolution modifications? Getting rid of the zeroGradient conditions on the faceplate didn't work.

[dkxls]

I just saw that you are actually using linear as div scheme for pretty much everything. That will most likely not work.

Again, for LES with an explicit SGS model I would recommend you the settings from the pitzDaily tutorial, i.e. use a TVD/NVD scheme!

In the pitzDaily tutorial they use a mix of filteredLinear2 and limitedLinear, which should stabilise your simulation significantly!

There are probably better choices for div schemes, but that again depends on your case and the pitzDaily tutorial is a really good starting point. There are numerous discussions on that topic here in the forum, just search for fvSchemes, FCT, TVD or the like. The topic is so old and widely discussed that you schould find information on this in pretty much every book on CFD. Here is a pretty good collection of the original references and book recommendations, this should get you started: http://en.wikipedia.org/wiki/Flux_limiter