Instability solving low-density plume expansion with rhoCentralFoam
 

Dear all,

I'm working on the simulation of an expanding nozzle plume using rhoCentralFoam. I am currently trying the replicate the results obtained by Ivanov et al. in http://arc.aiaa.org/doi/abs/10.2514/...ournalCode=jpp

They model the nozzle and the core flow with Navier-Stokes and then use DSMC to calculate the outer flow. I am only interested in the Navier-Stokes part. I've successfully created a case which computes the conditions inside the nozzle by using totalPressure and totalTemperature conditions for the inlet. For the outlet, I use a waveTransmissive boundary condition on the pressure, and zeroGradient on the temperature. Velocity is set to zeroGradient for both cases, as the nozzle flow is pressure-driven. The results of the simulation in OpenFOAM match the results in the paper almost exactly.

However, when I try to also model the expansion of the flow into vacuum, I run into problems. The domain I use is similar to the one used in the paper for continuum solutions. I initialize the pressure with a modified Lagrangian solver to obtain a smoothed pressure field (from 17 bars at the inlet to 1e-4 Pa at the outlet patches). When I run the case with the second order discretizations, it blows immediately. Using the following gradScheme, I am at least able to get it to run for 1e-4 seconds. This solution is relatively close to convergence, but there is severe numerical dissipation in it. This can be seen from the attached velocity field, with a very thick boundary.

No matter which schemes I use, the solution always diverges in the top left corner between the obstacle (left boundary) and the outlet (top boundary). At some point, the temperature becomes negative (or the density, I don't know), and the solver stops. Can anyone shed some light on why this happens and how I can avoid this behaviour? I have noticed that it happens faster on a fine grid with low numerical or physical dissipation, whereas coarse grids with high dissipation are more stable (but less accurate overall).

I am fully aware that a large part of this domain (especially the top left part near the nozzle lip) is too rarefied for Navier-Stokes to be accurate. However, this should not mean that the simulation should blow up, only that the solution cannot be trusted in this part of the domain. I have found out that the size of the domain strongly influences the behaviour of the core flow, so I do not want to simply make the domain smaller to avoid this behaviour.

If anyone could give me some hints for where to look, I would be very happy.

Hey all,

I've digged into this problem a bit deeper. It turns out that the authors of the original paper on rhoCentralFoam (Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows by Greenshields et al.) already recognized and mentioned this problem.

The issue is that temperature is not solved for directly, instead rhoCentralFoam solves for the energy per grid point. From this energy, the kinetic energy is subtracted to find the internal energy and thus the temperature:

The authors of the paper state that in one case, they also obtained negative temperatures and modified one of the interpolation schemes (or flux limiters, whatever you want to call them) to solve this problem. I also tried switching schemes and found that when using vanAlbada, the negative temperature occurs later than when using vanLeer. With upwind, it starts even later. Interestingly, the negative temperature occurs in a different location depending on the chosen limiter.

However, the problem at hand was not solved using this approach. In the end, it became clear to me that it was to be expected that negative temperature will almost unavoidably occur somewhere in the domain with this solver setup. As we are dealing with very low temperature and high velocity flows, the smallest error can cause a negative internal energy.

Therefore, I cheated somewhat. I implemented a minimum temperature by changing the line above to this:

I also write a small file that reads in a minTemperature from the controlDict file and translates it to a minimum energy value (minE), depending on the value of Cv (I use constant Cv and Cp).

I recognize that this is in fact cheating, but preliminary results indicate that this indeed makes it possible for the solver to converge to my validation case. I hope that it will be possible after initial transients to set the minimum temperature to zero without obtaining a negative temperature anywhere in the domain.

Hi LukaD,
Thanks for letting all know about the modification which wrks!
I to tried to do something similar to avoid this problem. I limited p to minimum value of 0.00001 pa. Can u please have a look? I have uploaded the files
http://www.cfd-online.com/Forums/ope...ntralfoam.html
Thanks,
Chirag

Hi lukaD,
can you pls share your modified solver ? my id is chiragkhalde@gmail.com
I will try to check whether it work for me too.
Thanks in advance!!
Chirag

Hi LukaD,

How do you pass the minE to the function that you have created?

Cheers,

Chris

Hi Chris,

I don't remember exactly (it was a long time ago) and I cannot access the code anymore as the work was for a previous employer. However, the lines above should give a clue. I made a new entry in controlDict, called "minTemperature" and somehow loaded that into the solver... But it was definitely not the hardest part for me, so you should be able to work it out I guess.

Good luck!