Rocket nozzle simulation with rhoCentralFoam
 

Dear Foamers

I have experience with CFX and want to switch over to OpenFOAM for my calculations. I keep running into the same problem and I am running out of idea's to fix it. (since am quite new to openfoam please call me out on any stupid mistakes I might have made).


I am currently trying to simulate a rocket nozzle with rhocentralfoam. For the velocity I have a fixedvalue inlet condition and an inletOutlet condition on the outlet. For the pressure I have a fixed value on the inlet and outlet. I am currently just running a laminar case to get the whole thing to work. I have made my mesh with netgen.


At normally the same time the solver crashed with the following message:


--> FOAM FATAL ERROR:
Maximum number of iterations exceeded


From function thermo<Thermo, Type>::T(scalar f, scalar T0, scalar (thermo<Thermo, Type>::*F)(const scalar) const, scalar (thermo<Thermo, Type>::*dFdT)(const scalar) const, scalar (thermo<Thermo, Type>::*limit)(const scalar) const) const
in file /home/hilbert/OpenFOAM/OpenFOAM-2.4.x/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.


FOAM aborting


#0 Foam::error:: printStack(Foam::Ostream&) at ??:?
#1 Foam::error::abort() at ??:?
#2 Foam::species::thermo<Foam::hConstThermo<Foam:: perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>::TEs(double, double, double) const at ??:?
#3 Foam::hePsiThermo<Foam: :psiThermo, Foam:: pureMixture<Foam::constTransport<Foam::species::th ermo<Foam::hConstThermo<Foam:; perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ??:?
#4 Foam::hePsiThermo<Foam::psiThermo, Foam:: pureMixture<Foam::constTransport<Foam::species::th ermo<Foam::hConstThermo<Foam:: perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ??:?
#5 ? at ??:?
#6 __libc_start_main in "/lib64/libc.so.6"
#7 ? at /home/abuild/rpmbuild/BUILD/glibc-2.18/csu/../sysdeps/x86_64/start.S:125


I have already tried to put the maximum number of iterations up, but that doesn't help.
From what I have read on this forum this error comes up when there are mesh issues.


To check this I have taken one of my meshes from my CFX calculations(which runs in CFX), and loaded that one into Openfoam. (slightly different geometry, same boundary conditions).


Now I get a very similar type of error.
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUx, Initial residual = 0, Final residual = 0, No Iterations 0
diagonal: Solving for rhoUy, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver: Solving for Ux, Initial residual = 4.90935e-10, Final residual = 4.90935e-10, No Iterations 0
smoothSolver: Solving for Uy, Initial residual = 1.10525e-10, Final residual = 1.10525e-10, No Iterations 0
diagonal: Solving for rhoE, Initial residual = 0, Final residual = 0, No Iterations 0
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib64/libc.so.6"
#3 Foam::hePsiThermo<Foam::psiThermo, Foam:pureMixture<Foam::sutherlandTransport<Foam::s pecies::thermo<Foam::hConstThermo<Foam::perfectGas <Foam::specie> >, Foam::sensibleEnthalpy> > > >::calculate() at ??:?
#4 Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::sutherlandTransport<Foam:: species::thermo<Foam::hConstThermo<Foam::perfectGa s<Foam::specie> >, Foam::sensibleEnthalpy> > > >::correct() at ??:?
#5 ? at ??:?
#6 __libc_start_main in "/lib64/libc.so.6"
#7 ? at /home/abuild/rpmbuild/BUILD/glibc-2.18/csu/../sysdeps/x86_64/start.S:125
Floating point exception


As a final check on my mesh I ran the following command on my mesh that ran in cfx(mesh generated with ICEMCFD):
checkMesh -allGeometry -allTopology


Checking geometry...
Overall domain bounding box (-11.505 -15 0) (11.505 36.6 1)
Mesh (non-empty, non-wedge) directions (1 1 0)
Mesh (non-empty) directions (1 1 0)
All edges aligned with or perpendicular to non-empty directions.
Boundary openness (-3.33549e-17 -4.40429e-17 1.23435e-14) OK.
Max cell openness = 2.63552e-14 OK.
Max aspect ratio = 368.724 OK.
Minimum face area = 3.82905e-09. Maximum face area = 0.125174. Face area magnitudes OK.
Min volume = 3.82905e-09. Max volume = 0.00618363. Total volume = 638.575. Cell volumes OK.
Mesh non-orthogonality Max: 19.9021 average: 8.25572
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.90249 OK.
Coupled point location match (average 0) OK.
Face tets OK.
*Edges too small, min/max edge length = 4.1539e-05 1, number too small: 2040
<<Writing 2040 points on short edges to set shortEdges
All angles in faces OK.
Face flatness (1 = flat, 0 = butterfly) : min = 1 average = 1
All face flatness OK.
Cell determinant (wellposedness) : minimum: 1.4384e-05 average: 2.44698
***Cells with small determinant (< 0.001) found, number of cells: 3546
<<Writing 3546 under-determined cells to set underdeterminedCells
Concave cell check OK.
Face interpolation weight : minimum: 0.0947609 average: 0.491795
Face interpolation weight check OK.
Face volume ratio : minimum: 0.104671 average: 0.971557
Face volume ratio check OK.


Failed 1 mesh checks.

Checkmesh seems to suggest a mesh issue, but if I look at my mesh statistics in icem is says that all element have a determinant between 0.95 and 1. Which would lead to the conclusion that there is an issue in the conversion script. (From icem I saved my mesh as a fluent file and used the fluentMeshToFoam utility to get my mesh)


This seems a bit unrealistic to me but I am a bit out of idea's. So if anybody have an other idea, that would be greatly appreciated.

This could be related to mesh issues, but does not have to (other examples: time stepping, wrong boundary conditions, ...). How does the solution look like just before the crash? Often, this helps in finding out where the error comes from.

Regards,

Martin

Hi Martin,

Thanks for your reply. That is were I am a bit puzzled, I can't find anything weird in the pressure, temp or velocity.

Checked time stepping, but I also get the same error with small time steps (1e-8).

Looking at the solution there is the most change at my inlet. At the my inlet I have a fixed value b.c for my velocity pressure and temperature, which looks oke to me.

Do you know of any other root causes for that error?

regards,
Hilbert

Well indeed if checkMesh fails then it means that there is a problem with the exporting process of your mesh. The error doesnt sound good.

Moreover, I would use a zeroGradient for your inlet p BC if you are planning to prescribe the inlet velocity

Is your case wedge?

You are using a compressible solver I suppose, so is the temperature 0 in somewhere of your domain?

Temperature looks fine. Also tried to put a total pressure bc on the inlet, but that was also to no avail. After some more experimentation on my mesh from icem I found the issue. The flow is accelerating in the first part of my domain. (first couple of mm) to very high speeds (were physically the flow shouldn't) and then the solver fails. I am still not really sure on what the root cause is though.

Is the flow subsonic at the inlet? In that case, I would set pressure and velocity at the same patch. I think this could be the source of error.

yes the flow is subsonic there. (I am initializing the flow subsonic through the entire domain).

Thanks for your tip, but i don't follow completely what you mean.

In my 2D domain my inlet and outlet are patches, my symmetry walls are empty, and the rest of my walls or of wall type.

On my Inlet I set both a fixedValue B.C on the pressure and velocity. Is this what you mean?

Hi, in my experience, rhoCentralFoam is not suitable for subsonic flows.

To overcome this limitation, i wrote hybrid solver, which you can download here

http://www.cfd-online.com/Forums/ope...ntralfoam.html

I'm interested in real life tests

mkraposhin, thanks for your reply, I will look into it. The problem turned out to be the mesh in the end and not the solver. I though I would post how I worked through my issues for anybody else who has the same issue.

The first point was the boundary conditions. I supplied a fixed value boundary condition on the inlet for the pressure velocity and temperature. This overspecifies the problem, and therefore one should only specify a fixed value B.C for 2 of the 3 and the 3the one should be a zeroGradient.

The problem that I was left with was that all the inlet moment was not convected downsteam. Therefore very high velocities were obtained directly at the inlet. The same problem has also been posted on this forum elsewhere: http://www.cfd-online.com/Forums/ope...outlet-bc.html (no solution here though) . I have now rebuild most of my mesh with blockmesh and now the solver runs like a charm. So somehow my mesh make with icem, that runs well in CFX does not run with openfoam. I have to note here that I did switch to sonicfoam in the process of problem solving.

I havn't figured out yet what the underlying issue is, but it is a mesh issue.

You have several cells with a determinant smaller than 0.001 in your mesh and you get a floating point error when using your cfx mesh.
I had a similar problem, although a different solver. Getting rid of the small determinants helped.

Hi TobM,

I gathered that from the checkmesh output. But I still find it a bit weird since I checked the determinanats in icem and icem told me that all determinants where greater that 0.9. Which in my view then means that the fluentMeshToFoam utility screws up the mesh.