[sebastian]

Hi all!


I am trying to perform a steady state simulation of a compressible flow through a nozzle with a development of a jet behind it, using the k-epsilon model.

I was playing around now for a couple of days and have tested a lot of setups but so far I did not get any solution. Hope any of you experts can help me with that.

MY CASE:
The flow of my current case is expected to be subsonic (Ma=0.5). But later on I want to observe also transonic flow (Ma=1.3), which I have not tested yet.
I am using a total pressure and a total temperature inlet at the inlets.
I want to use the standard-k-epsilon model, as well as the realizable k-epsilon model.
The case is steady state.

THE PROBLEM:
After a few hundred iterations, density and pressure start oscillating and then my solution blows up.
It seems that the disturbance occurs first at the outlet of the nozzle and than propagates back to the inlets...

WHAT I HAVE TRIED SO FAR:
First I have tested the incompressible case with simpleFoam. The solution looked pretty okay.
Then I switched to the compressible case using rhoSimpleFoam.
I have tested first order as well as second order discretization. I also played with the relaxation factors and the relative solver tolarance.
I also have tested several meshes.

I have calculated the case in Fluent, using the compressible SIMPLE algorithm without any problems. This should be the same as rhoSimpleFoam right? So why the problems in OF?

Well, okay maybe OF has a problem with that and rhoSimpleFoam not the adequate solver for that.. I am a little bit lost to be honest.

If anybody is interested, I have uploaded my case file.


Would be glad to hear some ideas


Thanks in advance!

Sebastian

[FelixL]

Hi, Sebastian,


is checkMesh reporting any kind of troubles regarding your mesh? A picture of the mesh for your case would be helpful. You said, the disturbance seems to first occur at the outlet, so it might be logical to search for the error source at the outlet first. What kind of boundary conditions have you set there?

Another thing: have you tried running the case with turbulence switched off? If the simulation still blows up, we may know that it's not (only) the turbulence model, which is causing the problems.


Greetings,
Felix

[sebastian]

Hello Felix!


Thanks for the quick answer!

Yes, the problem is the same without turbulence. checkmesh says everything is okay.

Sorry, maybe my explanation was a bit misleading. I didnt mean that the problem occurs at the boundary outlet of my domain. It occurs at some state of the calculation at the exit of the nozzle when the flow enters into the farfield behind the nozzle.

I have attached some pics of the whole domain and the nozzle configuration.
I also attached of the distribution of the static pressure, a few iterations before the simulation stops. The pressure distribution is not physical. It should be highest at the inlets. The pressure at the nozzle exit should be the same as ambient...


Hope this helps to understand the case better!


Sebastian

[FelixL]

Hi, Sebastian,


at first sight I can't really tell, what's causing the problem. The mesh looks great, I don't see any troublemakers there. Is it the whole domain you're showing in the first picture? If so, it probably might be better if you move the outlet boundary further downstream where the jet has dissipated more.

Could you post your BCs, please? (all of them, including velocity) And a velocity or Mach-number plot of the nozzle would be great, as it's a more "natural" way for interpreting the flow.


Greetings,
Felix.

[sebastian]

Hello Felix!


Yes, its the whole domain. I know its small. But for now this should be okay. I think this doesnt cause my problems.

Here my BCs. And a plot of the velocity-magnitude. Think for this time of progress of the simulation, the velocity looks good. Max velocity here is ca 110 m/s, maximum Mach number is 0.35.
Do you think the rhoSimpleFoam solver should handle this in general?

Best regards,

Sebastian

[FelixL]

Hey, there,


I had a look at your files and images: to my understanding everything seems to be setup fine. I have to say that I'm not too familiar with compressible simulations in OpenFOAM, but the totalPressure BCs along with zeroGradient velocity seem to be reasonable to me. Also you don't seem to have any backflow issues at the boundaries, so that shouldn't be a problem here.

Another thought here: I had a similiar issue with that "choked" pressure distribution one can see in the picture you attached at post #3. That issue was in an incompressible case, though, and was a result of a local mesh issue, if I recall correctly. I was able to solve the problem by setting

gradSchemes
{
    default         faceMDLimited Gauss linear 0.5;
}

inside fvSchemes. Have you already tried that?

Hope that helps,

Greetings,
Felix.

[sebastian]

Hi Felix!

Thanks! No, I have not tried that. What exactly happens here if I use that grad scheme?

I will give it a try


Best regards,

Sebastian

[FelixL]

Hey,


it's a limited version of the Gauss linear gradScheme. It assures that the extrapolated face values are within the limits of the neighboring cell values:

Description
    faceMDLimitedGrad gradient scheme applied to a runTime selected
    base gradient scheme.

    The scalar limiter based on limiting the extrapolated face values
    between the face-neighbour cell values and is applied to the gradient
    in each face direction separately.

Good luck with that,

Greetings,
Felix.

[Chris Lucas]

Hi,

first of all, could you please post the full case (without mesh).

Are you using the transonic setting? If not, your pressure equation (at your high velocities) is not correct (as I understand this solver).

Regards,
Christian

[sebastian]

Hi!

I have tried the limitation of the gradient now. It seems that the calculation is a bit more stable and it runs a bit longer. But the result is the same - same problem occurs at the exit of my nozzle...


@Christian: what do you mean by transonic setting?

I have uploaded now the whole case again. So anybody feel free to tell me their ideas about the problem! Thanks!


Best regards,
Sebastian

[sebastian]

Hi!

I think I found the option transonic. So I switched it on by changing fvSolution:

SIMPLE
{
nNonOrthogonalCorrectors 0;
pMin	pMin [1 -1 -2 0 0 0 0] 50000;	
transonic true;
}

The calculation with that option seems to me highly unstable. In every configuration, it doesnt last more than a few iterations before it blows up:



Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.00064076898, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.0011129843, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.0014193024, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99999977, Final residual = 1.2279612e-05, No Iterations 1
GAMG: Solving for p, Initial residual = 1, Final residual = 0.00010163833, No Iterations 1
time step continuity errors : sum local = 5249.3272, global = 5249.2913, cumulative = 5249.2913
rho max/min : 1.4641599 1.2200726
ExecutionTime = 44.24 s ClockTime = 45 s

Time = 2

DILUPBiCG: Solving for Ux, Initial residual = 0.063224505, Final residual = 0.0014127685, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.77219605, Final residual = 0.018765348, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.38431967, Final residual = 0.010036484, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.9805994, Final residual = 0.00013796488, No Iterations 1
GAMG: Solving for p, Initial residual = 0.61791917, Final residual = 2.2252849e-05, No Iterations 1
time step continuity errors : sum local = 35648.227, global = 35647.963, cumulative = 40897.254
rho max/min : 1.4641542 1.2045404
ExecutionTime = 74.99 s ClockTime = 76 s

Time = 3

DILUPBiCG: Solving for Ux, Initial residual = 0.19958378, Final residual = 0.0038654242, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.524207, Final residual = 0.009996843, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.17581458, Final residual = 0.0036239134, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.98209612, Final residual = 0.00020747175, No Iterations 1
GAMG: Solving for p, Initial residual = 0.60087116, Final residual = 1.5788899e-05, No Iterations 1
time step continuity errors : sum local = 43254.307, global = 43254.274, cumulative = 84151.528
rho max/min : 1.4641489 1.1446924
ExecutionTime = 105.68 s ClockTime = 106 s

Time = 4

DILUPBiCG: Solving for Ux, Initial residual = 0.25990016, Final residual = 0.0059634889, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.32796848, Final residual = 0.0048554188, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.095194627, Final residual = 0.0014654287, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99835796, Final residual = 0.00022092011, No Iterations 1
GAMG: Solving for p, Initial residual = 0.57627831, Final residual = 3.4923119e-06, No Iterations 1
time step continuity errors : sum local = 126889.55, global = 126889.52, cumulative = 211041.05
rho max/min : 1.4641441 1.0874578
ExecutionTime = 136.46 s ClockTime = 137 s

Time = 5

DILUPBiCG: Solving for Ux, Initial residual = 0.10431322, Final residual = 0.00092621411, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.28540798, Final residual = 0.0043205614, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.049779571, Final residual = 0.00081063858, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99972076, Final residual = 0.00016023754, No Iterations 1
GAMG: Solving for p, Initial residual = 0.4916998, Final residual = 2.2014033e-07, No Iterations 1
time step continuity errors : sum local = 1955369.2, global = 1955369.2, cumulative = 2166410.2
rho max/min : 1.4641398 1.033085
ExecutionTime = 167.13 s ClockTime = 168 s

Time = 6

DILUPBiCG: Solving for Ux, Initial residual = 0.0086557065, Final residual = 4.4554317e-05, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.046355992, Final residual = 0.0004336927, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.01445842, Final residual = 0.00019740974, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99998175, Final residual = 0.0001601321, No Iterations 1
GAMG: Solving for p, Initial residual = 0.41766058, Final residual = 2.0001359e-08, No Iterations 1
time step continuity errors : sum local = 22979085, global = 22979085, cumulative = 25145496
rho max/min : 1.4641358 0.98143071
ExecutionTime = 197.77 s ClockTime = 198 s

Time = 7

DILUPBiCG: Solving for Ux, Initial residual = 0.0010534608, Final residual = 2.3964112e-06, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.010440758, Final residual = 4.4895448e-05, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0029302907, Final residual = 1.7381004e-05, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99999912, Final residual = 0.00015988003, No Iterations 1
GAMG: Solving for p, Initial residual = 0.45276891, Final residual = 2.7486885e-08, No Iterations 1
time step continuity errors : sum local = 58809550, global = 58809550, cumulative = 83955046
rho max/min : 1.4641322 0.93235917
ExecutionTime = 228.46 s ClockTime = 229 s

Time = 8

DILUPBiCG: Solving for Ux, Initial residual = 0.00038112419, Final residual = 5.4734558e-06, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0069200946, Final residual = 7.5209829e-05, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0021328287, Final residual = 2.5482588e-05, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99999963, Final residual = 0.00015815714, No Iterations 1
GAMG: Solving for p, Initial residual = 0.52127838, Final residual = 4.6696771e-09, No Iterations 1
time step continuity errors : sum local = 65978178, global = 65978178, cumulative = 1.4993322e+08
rho max/min : 1.4641289 0.88574121
ExecutionTime = 259.08 s ClockTime = 260 s

Time = 9

DILUPBiCG: Solving for Ux, Initial residual = 0.00023764328, Final residual = 3.953079e-07, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 0.0054798111, Final residual = 0.00020048139, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0025317365, Final residual = 7.6321079e-05, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99999965, Final residual = 0.00015632946, No Iterations 1
GAMG: Solving for p, Initial residual = 0.51953995, Final residual = 2.1894755e-09, No Iterations 1
time step continuity errors : sum local = 67810556, global = 67810556, cumulative = 2.1774378e+08
rho max/min : 1.464126 0.84145415
ExecutionTime = 290.13 s ClockTime = 291 s

Time = 10

DILUPBiCG: Solving for Ux, Initial residual = 0.00021875401, Final residual = 1.5713824e-05, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0042287586, Final residual = 1.0771838e-05, No Iterations 2
DILUPBiCG: Solving for Uz, Initial residual = 0.002683303, Final residual = 0.00025108918, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99999967, Final residual = 0.00015442018, No Iterations 1
GAMG: Solving for p, Initial residual = 0.48863879, Final residual = 1.6088557e-09, No Iterations 1
time step continuity errors : sum local = 68560890, global = 68560890, cumulative = 2.8630467e+08
rho max/min : 1.4641233 0.79938145
ExecutionTime = 321.16 s ClockTime = 322 s

Time = 11

DILUPBiCG: Solving for Ux, Initial residual = 0.00017913139, Final residual = 2.244739e-06, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0041400616, Final residual = 8.179217e-05, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0027678658, Final residual = 4.9276336e-05, No Iterations 1
DILUPBiCG: Solving for h, Initial residual = 0.99999968, Final residual = 0.00015256293, No Iterations 1
[0]
[0]
[0] Maximum number of iterations exceeded#0 Foam::error:rintStack(Foam::Ostream&) in "/home/sebastian/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/sebastian/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Foam::hPsiThermo<Foam:ureMixture<Foam::constTran sport<Foam::specieThermo<Foam::hConstThermo<Foam:: perfectGas> > > > >::calculate() in "/home/sebastian/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#3 Foam::hPsiThermo<Foam:ureMixture<Foam::constTran sport<Foam::specieThermo<Foam::hConstThermo<Foam:: perfectGas> > > > >::correct() in "/home/sebastian/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#4 main in "/home/sebastian/OpenFOAM/OpenFOAM-1.6/applications/bin/linuxGccDPOpt/rhoSimpleFoam"
#5 __libc_start_main in "/lib/tls/libc.so.6"
#6 _start at ../sysdeps/i386/elf/start.S:122
[0]
[0]
[0] From function specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const
[0] in file /home/dm2/henry/OpenFOAM/OpenFOAM-1.6/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 68.
[0]
FOAM parallel run aborting
[0]
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun has exited due to process rank 0 with PID 18347 on
node node02 exiting without calling "finalize". This may
have caused other processes in the application to be
terminated by signals sent by mpirun (as reported here).
--------------------------------------------------------------------------

[FelixL]

Hi, Sebastian,


well, due to my limited experience with the rhoSimpleFoam-solver this case is beginning to give me a headache. I can only guess now and I think that's not the help you're looking for - nevertheless: two more thoughts that came to my mind:

- What relaxation factors did you try? Maybe this thread gives you a new kind of inspiration: http://www.cfd-online.com/Forums/ope...am-solver.html

- Your initial field is a uniform field, as far as I can see. You might want to try to map the simpleFoam results to your compressible case (which is pretty tough because of the normalized pressure values) or you can try getting a starting field with a short transient calculation and have rhoSimpleFoam start a calculation with the resulting initial values.


Sorry about the guesswork, but I'm pretty confident someone else here knows a fine answer.


Greetings,
Felix.

[Chris Lucas]

Hi,

the reason why the simulation crashes is that either your pressure or your enthalpy is negative, so the Newton solver can't find the correct T for your H.

rhoSimpleFoam is not the best solver for high speed steady state simulations (as far as I now). For this kind of simulations, a coupled solver is needed.

I would use a transient solver (rhoPisoFoam). This will take some time, but you will get result.

Regards,
Christian

[sebastian]

Hi Chris!

Thanks for the reply!

Whish there was installed a coupled solver here in OpenFoam. Of course I would have used that one instead.
I am currently running a transient calculation with the piso solver and it takes me ages. But why is Fluent able to handle the prob with the simple algorithm? Does anybody know that? Or do they just use some fancy numerical tricks?


Best regards,

Sebastian

[cleoo]

Originally Posted by Chris Lucas

I would use a transient solver (rhoPisoFoam). This will take some time, but you will get result.

Hi Chris,

since there is no more rhoPisoFoam do you think sonicFoam or rhoCentralFoam might be ok with this same case?

Also, Sebastian did the transient case work out for you?