[gallon]

Hello everyone, here's a problem I've been trying to solve for a long time.
I have an stl model of a bullet and I want to calculate its supersonic movement under normal conditions.
However at some moment the solver yields negative pressure in the trail => negative temperature => fails.


Now for the details.
1. Case setup.
compressible/sonicFoam/ras/nacaAirfoil was used as an example/template for BC.
For time step I define ~1e-9 by default, but I tried steps down to 1e-12, the result is the same. With these steps, inlet velocity 850 m/s and cell size down to 1e-4 m the Courant number is about ~0.001 which should be more than sufficient.

2. Mesh
I tried sHM, cfMesh and ICEM CFD. Layers or no layers, the general behaviour is the same: the larger the cell size and the smaller the velocity, the farther the computation goes before failing. At velocity ~400 m/s and cell size ~1-2 mm the computations manages to come to a steady solution but it is wrong because the shock wave is not resolved.

3. The error
Here's the output:

[185] --> FOAM FATAL ERROR:
[185] Maximum number of iterations exceeded
[185]
[185]     From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar)const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar)const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar)const) const [with Thermo = Foam::hConstThermo<Foam::perfectGas<Foam::specie> >; Type = Foam::sensibleInternalEnergy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy>]
[185]     in file /common/software/OpenFOAM/OpenFOAM-4.0/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 69.

I opened thermoI.H, added some printfs, recompiled and rerun. The T is negative at that point, so it is not the root of the problem.

4. The search for success
What else I tried?.. Switching schemes in fvSchemes, a hundred different mesh setups, adjusting timestep

Yet I took `nacaairfoil` example, made it 3D by changing empty->zeroGradient and refineMesh command and directed the inlet normally to the wing plane to test. Well, it worked fine! Even at 7 Machs...

Here's the full case setup attached. The bullet.stl is changed to a simple handmade model for understandable reasons. I checked, the result is the same. PS controlDict is for OF 4. You can disable postprocessing for compatibility.
To run:

blockMesh
surfaceFeatureExtract
mesh
sonicFoam

Decomposition is strongly advised for perfomance.
bullet.zip

PPS I tried this on another model but there's a good chance that if you flip the wind (exchange inlet with outlet and swap sign in 0/defines/U), everything will go well.

[KarenRei]

I'm having the same problem with sonicFoam I'm not as far along as you at finding solutions as it's taken me days to compile the latest version of OpenFOAM (I'm on Fedora, so building it is always a pain, the makefiles are missing libs, some libs get put in the wrong place, etc)

Temperature = -0.22K on my last run for some reason... low courant number... no luck.

[KarenRei]

Found the solution in my case I dropped the p tolerance in fvsolutions by two orders of magnitude, that totally stabilized the run (and let me dramatically increase my timestep!)

[gallon]

Originally Posted by KarenRei

Found the solution in my case I dropped the p tolerance in fvsolutions by two orders of magnitude, that totally stabilized the run (and let me dramatically increase my timestep!)

Thanks for the info. Didn't work for me, unfortunately.

[KarenRei]

Originally Posted by gallon

Thanks for the info. Didn't work for me, unfortunately.

Unfortunate. I've been doing a lot of different runs with different solvers, and randomly running into such negative temperature excursions... and each time it really seems like it's just guesswork to find things that work. It's almost always by some combination of tweaking fvSchemes, fvSolutions, and increasing the mesh resolution in sensitive areas where oscillatory or turbulent behavior spirals out of control. But I've not yet found some sort of easy, always-works solution. As you noted, shrinking the timestep rarely seem to be the solution if the courant number is already sufficiently low.

What's your fvSchemes that you're using right now? And what does "checkMesh -constant -allGeometry -allTopology" say? No problems with the mesh? Orthogonality? Skewness?

Oh, one more thing: one "excursion" case I had was due to me setting up the geometry wrong; I thought I had linked some faces up proper boundary conditions but it had been the wrong faces, leaving some faces with boundaries that they shouldn't have had and some faces with no boundary conditions set whatsoever. And it led to some really weird behavior. So probably good to double check your geometry - sometimes geometry can "pass" all checks but still be wrong.

[gallon]

Originally Posted by KarenRei

It's almost always by some combination of tweaking fvSchemes, fvSolutions, and increasing the mesh resolution in sensitive areas where oscillatory or turbulent behavior spirals out of control. But I've not yet found some sort of easy, always-works solution. As you noted, shrinking the timestep rarely seem to be the solution if the courant number is already sufficiently low.

The weird thing is that the finer I make the mesh the quicker negative p/T is achieved, even if timestep is reduced accordingly... while the solution is generally expected to converge better on a finer grid.

Originally Posted by KarenRei

What's your fvSchemes that you're using right now? And what does "checkMesh -constant -allGeometry -allTopology" say? No problems with the mesh? Orthogonality? Skewness?
Oh, one more thing: one "excursion" case I had was due to me setting up the geometry wrong; I thought I had linked some faces up proper boundary conditions but it had been the wrong faces, leaving some faces with boundaries that they shouldn't have had and some faces with no boundary conditions set whatsoever. And it led to some really weird behavior. So probably good to double check your geometry - sometimes geometry can "pass" all checks but still be wrong.

I guess I tried almost every combination of numeric schemes
checkMesh usually tells everything is ok. I also tried both real STL from CAD and a simple self-made "projectile" formed as a half of an ellipsoid.
The geometry looks rather simple: a cylinder (one circle for inlet, other parts for outlet) with a wall inside...

If only someone reported the reason of this weird behavior!

[Vyssion]

One thing that I have used successfully in the past is to set up a temperature limit within an fvOptions file:

TemperatureLimit1
{
	type		temperatureLimitsContraint;
	active		true;
	selectionMode	all;

	temperatureLimitsConstraintCoeffs
	{
		Tmin	 200;
		Tmax	 2000;
	}
}

[KarenRei]

Ooh neat - didn't know that was possible!

[KarenRei]

Originally Posted by Vyssion

One thing that I have used successfully in the past is to set up a temperature limit within an fvOptions file:

TemperatureLimit1
{
	type		temperatureLimitsContraint;
	active		true;
	selectionMode	all;

	temperatureLimitsConstraintCoeffs
	{
		Tmin	 200;
		Tmax	 2000;
	}
}

Hmm. This isn't working for me.

First off, apparently in newer versions of OpenFoam it's now "limitTemperature", not "temperatureLimitsConstraint". To get it to start running, I ended up having to go with (in system/fvOptions):

temperatureLimit
{
type limitTemperature;
active true;

limitTemperatureCoeffs
{
selectionMode all;
min 100;
max 5000;
}
}

(I don't know what's supposed to go on that "temperatureLimit line, I think it's just a name?)

It's clearly parsing this file because if I have any sort of typo it doesn't start. But the limits don't actually seem to limit anything:

--> FOAM Warning :
From function Foam::scalar Foam::janafThermo<EquationOfState>::limit(Foam::sc alar) const [with EquationOfState = Foam:erfectGas<Foam::specie>; Foam::scalar = double]
in file /opt/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/specie/lnInclude/janafThermoI.H at line 105
attempt to use janafThermo<EquationOfState> out of temperature range 1e-06 -> 6000; T = -1.3418404
--> FOAM Warning :
From function Foam::scalar Foam::janafThermo<EquationOfState>::limit(Foam::sc alar) const [with EquationOfState = Foam:erfectGas<Foam::specie>; Foam::scalar = double]
in file /opt/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/specie/lnInclude/janafThermoI.H at line 105
attempt to use janafThermo<EquationOfState> out of temperature range 1e-06 -> 6000; T = -1.3418404
--> FOAM Warning :
From function Foam::scalar Foam::janafThermo<EquationOfState>::limit(Foam::sc alar) const [with EquationOfState = Foam:erfectGas<Foam::specie>; Foam::scalar = double]
in file /opt/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/specie/lnInclude/janafThermoI.H at line 105
attempt to use janafThermo<EquationOfState> out of temperature range 1e-06 -> 6000; T = -1.3418404


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

From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::janafThermo<Foam:erfectGas<Foam::specie> >; Type = Foam::sensibleInternalEnergy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::janafThermo<Foam:erf ectGas<Foam::specie> >, Foam::sensibleInternalEnergy>]
in file /opt/OpenFOAM/OpenFOAM-dev/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.

FOAM aborting

Am I doing something wrong?

[Vyssion]

Originally Posted by KarenRei

First off, apparently in newer versions of OpenFoam it's now "limitTemperature", not "temperatureLimitsConstraint"

I don't use the latest version of OpenFOAM so it may be different yes.

Originally Posted by KarenRei

--> FOAM Warning :
From function Foam::scalar Foam::janafThermo<EquationOfState>

This janafThermo error is a nightmare to deal with sometimes - your boundary conditions and/or thermalProperties script may not have the correct thermal properties for your setup.

There isn't too much help that I can offer with this issue tbh... I did come across this thread which someone was whinging about it aimlessly - but it did have some links to other threads which apparently helped other people, so I wish you luck... sorry

[cfdsolver1]

I am having kind of same problem with chtMultiRegionSimpleFoam. Negative temperature is a problem in my case.

[gallon]

Originally Posted by KarenRei

Hmm. This isn't working for me.

True, for OF 4.0 syntax has to be adjusted, but it's done easily.

I track min/max p/U/T using postprocessing, so it worked for me like this (p/T/U respectively):


Looks like the constraint tries to work but is applied too seldom or something.

You can try to add this to controlDict to get minmax postprocessing:

functions
{
    minMax
    {
        type            fieldMinMax;
        libs ("libfieldFunctionObjects.so");
        enabled         true;
        log             true;
        write           true;
        // Fields to be monitored - runTime modifiable
        fields(T U p);
    }
}

This is the script I wrote to view data with gnuplot in realtime:

#!/usr/bin/bash
set -e

if [ $# != 1 ]; then
  echo "Usage: gp.sh <FILE>"
  exit -1
fi

glaunch () {
cat << EOF > gpscript."$2"
  plot "$1" every 3::$2 using 1:3 title 'x' with points, "$1" every 3::$2 using 1:8 title 'x' with points
  pause 1
  reread
EOF

gnuplot ./gpscript."$2"
}

glaunch $1 0 &
glaunch $1 1 &
glaunch $1 2 &

Simply issue something like:

$ ./this_script.sh ./postProcessing/minMax/0/fieldMinMax.dat

having launched the solver or after it finishes.

[gallon]

Good news everybody. I tweaked system/fvSchemes so that the solver doesn't fail.

gradSchemes
{
    default cellLimited Gauss linear 1; //default: Gauss linear;
}

divSchemes
{
    div(phi,U) Gauss linearUpwind cellLimited Gauss linear 1; // default: Gauss limitedLinearV 1;
}

laplacianSchemes
{
    default Gauss linear limited 0.333; // default: Gauss linear limited corrected 0.5;
}

snGradSchemes
{
    default limited 0.333; // default: corrected;
}

Honestly, I have no idea how and why this works but it works. Probably some of these changes are redundant but right now I am not ready to try it. Anyway, it works as is.

Perhaps fellow members could explain why.

[Icemaan]

Originally Posted by gallon

Good news everybody. I tweaked system/fvSchemes so that the solver doesn't fail.

gradSchemes
{
    default cellLimited Gauss linear 1; //default: Gauss linear;
}

divSchemes
{
    div(phi,U) Gauss linearUpwind cellLimited Gauss linear 1; // default: Gauss limitedLinearV 1;
}

laplacianSchemes
{
    default Gauss linear limited 0.333; // default: Gauss linear limited corrected 0.5;
}

snGradSchemes
{
    default limited 0.333; // default: corrected;
}

Honestly, I have no idea how and why this works but it works. Probably some of these changes are redundant but right now I am not ready to try it. Anyway, it works as is.

Perhaps fellow members could explain why.

Hi,

I had the same negative temperature issue while solving for Mach 2 external flow using sonicFoam. I did a similar thing. I guess one of the reason this tweak works is because of the order of the divergence schemes. I only changed the divergence scheme from 'llimitedLinear V 1' (which is second order) to 'Gauss Upwind' (which is first order). This took care of this problem for my simulation. Also, after obtaining a good enough solution with first order scheme, I was able to change it to second order and continue, although at a smaller timestep. I did get the same negative temperature problem again with the second order simulation, but this time it was only due to trying to up the timestep too fast.

[Paul Caicedo]

Originally Posted by gallon

True, for OF 4.0 syntax has to be adjusted, but it's done easily.


Looks like the constraint tries to work but is applied too seldom or something.

You can try to add this to controlDict to get minmax postprocessing:

functions
{
    minMax
    {
        type            fieldMinMax;
        libs ("libfieldFunctionObjects.so");
        enabled         true;
        log             true;
        write           true;
        // Fields to be monitored - runTime modifiable
        fields(T U p);
    }
}

Hello Vladislav,

Perhaps, do you know how to set this in version 5.x? As an example, in previous versions to use MRF you needed to include the data in the system/fvOptions file, but now it is located in constant/MRFProperties file.

Regards,
Paul

[foadsf]

Adding a fvOptions file in system folder with :

temperatureLimit
{
    type limitTemperature;
    active true;

    limitTemperatureCoeffs
    {
        selectionMode all;
        min <Tmin>;
        max <Tmax>;
    }
}

seems to resolve the issue.

[thanosniotis]

Originally Posted by foadsf

Adding a fvOptions file in system folder with :

temperatureLimit
{
    type limitTemperature;
    active true;

    limitTemperatureCoeffs
    {
        selectionMode all;
        min <Tmin>;
        max <Tmax>;
    }
}

seems to resolve the issue.

That is good! I will implement it and we will see what happens.
But I suspect that the problem with the negative temperature exists when I use turbulence flow (Kepsilon), in case of laminar flow everything works fine. I dont know if any fellow foamer has noticed something relevant.

[Chintan-21]

Hello guys


I am running multi phase simulation inside the helical pipe, and I am getting the error as follow for the volume fraction.

PIMPLE: Iteration 2
MULES: Solving for alpha.gas
MULES: Solving for alpha.liquid
gas fraction, min, max = 0.002992778214 -2.60889491e-06 1
liquid fraction, min, max = 0.997007311 -1.776623502e-36 1.000002609
Phase-sum volume fraction, min, max = 1.000000089 -7.771561172e-16 1.000721991
MULES: Solving for alpha.gas
MULES: Solving for alpha.liquid
gas fraction, min, max = 0.002992170717 -2.646514245e-06 1
liquid fraction, min, max = 0.9970079185 -1.792762148e-36 1.000002647
Phase-sum volume fraction, min, max = 1.000000089 -7.771561172e-16 1.000721991
MULES: Solving for alpha.gas
MULES: Solving for alpha.liquid
gas fraction, min, max = 0.002991565366 -2.629352244e-06 1
liquid fraction, min, max = 0.9970085238 -1.78660011e-36 1.000002629
Phase-sum volume fraction, min, max = 1.000000089 -7.771561172e-16 1.000721991
Constructing momentum equations
Tf.gas_liquid: min = 289.0432295, mean = 292.9999477, max = 293.0876964
DILUPBiCGStab: Solving for O2.gas, Initial residual = 0.001402672873, Final residual = 1.716052678e-13, No Iterations 1
DILUPBiCGStab: Solving for water.gas, Initial residual = 0.002668600287, Final residual = 1.226195265e-16, No Iterations 1
DILUPBiCGStab: Solving for N2.liquid, Initial residual = 2.366996742e-05, Final residual = 1.318267e-16, No Iterations 1
DILUPBiCGStab: Solving for O2.liquid, Initial residual = 4.701264042e-05, Final residual = 4.83641327e-16, No Iterations 1
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:perator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#5 Foam::TwoResistanceHeatTransferPhaseSystem<Foam::M omentumTransferPhaseSystem<Foam:haseSystem> >::heatTransfer() const at ??:?
#6 Foam::InterfaceCompositionPhaseChangePhaseSystem<F oam::PhaseTransferPhaseSystem<Foam::TwoResistanceH eatTransferPhaseSystem<Foam::MomentumTransferPhase System<Foam:haseSystem> > > >::heatTransfer() const at ??:?
#7 Foam::solvers::multiphaseEuler::energyPredictor() at ??:?
#8 Foam::solvers::multiphaseEuler::thermophysicalPred ictor() at ??:?
#9 ? in "/opt/openfoam11/platforms/linux64GccDPInt32Opt/bin/foamRun"
#10 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#11 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12 ? in "/opt/openfoam11/platforms/linux64GccDPInt32Opt/bin/foamRun"
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:perator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#5 Foam::TwoResistanceHeatTransferPhaseSystem<Foam::M omentumTransferPhaseSystem<Foam:haseSystem> >::heatTransfer() const at ??:?
#6 Foam::InterfaceCompositionPhaseChangePhaseSystem<F oam::PhaseTransferPhaseSystem<Foam::TwoResistanceH eatTransferPhaseSystem<Foam::MomentumTransferPhase System<Foam:haseSystem> > > >::heatTransfer() const at ??:?
#7 Foam::solvers::multiphaseEuler::energyPredictor() at ??:?
#8 Foam::solvers::multiphaseEuler::thermophysicalPred ictor() at ??:?
#9 ? in "/opt/openfoam11/platforms/linux64GccDPInt32Opt/bin/foamRun"
#10 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#11 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12 ? in "/opt/openfoam11/platforms/linux64GccDPInt32Opt/bin


even I dropped the pressure tolerance in fvSolution, but still I am getting the error

any help would be greatly appreciated.