CFD Online Logo CFD Online URL
Home > Forums > OpenFOAM Running, Solving & CFD

Instability solving low-density plume expansion with rhoCentralFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Display Modes
Old   November 5, 2014, 10:05
Default Instability solving low-density plume expansion with rhoCentralFoam
New Member
Luka Denies
Join Date: Oct 2014
Posts: 23
Rep Power: 3
LukaD is on a distinguished road
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

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.

    default         cellLimited Gauss linear 1; 
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.
Attached Images
File Type: jpg nozzleExpansion_p_t136e-6.jpg (13.5 KB, 23 views)
File Type: jpg nozzleExpansion_rho_t136e-6.jpg (13.6 KB, 20 views)
File Type: jpg nozzleExpansionT_t130e-6.jpg (14.4 KB, 16 views)
File Type: jpg nozzleExpansionT_t136e-6.jpg (15.3 KB, 17 views)
File Type: jpg nozzleExpansionU_t130e-6.jpg (13.4 KB, 18 views)
LukaD is offline   Reply With Quote

Old   November 10, 2014, 04:13
New Member
Luka Denies
Join Date: Oct 2014
Posts: 23
Rep Power: 3
LukaD is on a distinguished road
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:

e = rhoE/rho - 0.5*magSqr(U);
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:

e = max(rhoE/rho - 0.5*magSqr(U),minE);
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.
LukaD is offline   Reply With Quote


negative density, nozzle, plume, rhocentralfoam

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
buoyantSimpleFoam and watertank Tobi OpenFOAM Running, Solving & CFD 48 December 26, 2014 09:49
Moving mesh Niklas Wikstrom (Wikstrom) OpenFOAM Running, Solving & CFD 122 June 15, 2014 06:20
rhoSimplecFoam Mach0.8 no pressure values CFDnewbie147 OpenFOAM Running, Solving & CFD 16 November 23, 2013 06:58
Why RNGkepsilon model gives floating error shipman OpenFOAM Running, Solving & CFD 3 September 7, 2013 08:00
High Courant Number @ icoFoam Artex85 OpenFOAM Running, Solving & CFD 9 January 3, 2012 09:06

All times are GMT -4. The time now is 09:59.