temperature negative occurred in rhoCentralFoam at high M

longming · August 29, 2018, 09:06

Hi,guys

When I use rhocentralFoam to solve flow over wedge or other simple shape( using RAS model), negative temperature occurred easily when M>5, even dt and maxco are set to very small value.
I followed the example in tutorial, it works when computing grid is course. I only refined the grid in boundary layer and the problem just occurred.
Is the slover unsteady at high M?
Who knows how to slove the problem? Please Help

hokhay · October 22, 2018, 22:19

Hi Longming,

I am encountering the same problem. Have you found a solution yet?

Jason

backscatter · October 23, 2018, 00:37

Quote:

Originally Posted by longming

Hi,guys

When I use rhocentralFoam to solve flow over wedge or other simple shape( using RAS model), negative temperature occurred easily when M>5, even dt and maxco are set to very small value.
I followed the example in tutorial, it works when computing grid is course. I only refined the grid in boundary layer and the problem just occurred.
Is the slover unsteady at high M?
Who knows how to slove the problem? Please Help

I've had such kind of issues in passive scalar transport in incompressible solvers though. Usual practice is to try to use a bounded scheme (div and grad) for your energy/TEqn. You can start with upwind which is more stable as it incurs numerical dissipation which stablizes the run, but it could be less accurate. Relaxing the equation can also help in some cases. In my case, the culprits were high-aspect ratio cells away from the boundary layer- When I removed those high aspect cells, my problem was solved.

Here's the link to my post:
TEqn (passive scalar) diverges all of a sudden

TheMadHungarian · October 23, 2018, 16:40

Usually setting deltaT to a very small number (1e-15 for the 1st time step) will fix this, the solver will adjust the timestep to a larger value automatically.

chc · January 21, 2024, 16:00

Quote:

Originally Posted by longming

Hi,guys

When I use rhocentralFoam to solve flow over wedge or other simple shape( using RAS model), negative temperature occurred easily when M>5, even dt and maxco are set to very small value.
I followed the example in tutorial, it works when computing grid is course. I only refined the grid in boundary layer and the problem just occurred.
Is the slover unsteady at high M?
Who knows how to slove the problem? Please Help

Hello,
I'm currently experiencing similar issues. I'm new to OpenFOAM still and was wondering if you found a solution to this and where? Thank you!.

chemak · January 22, 2024, 04:51

Hi,
I think that the cause of negative temerature in rhoCentralFoam is rather physical than numerical. Check your ground energy in "thermophysicalProperties" (Cp, Tref and Href), setting more physical values will generally prevent this problem.

chc · January 22, 2024, 21:05

My initial set up uses the settings below for that file. I've taken them from relevant articles that reference sea level values for air. I'm solving using dimensional initial values for my velocity, temperature, and pressure, so I'm not sure where to set Tref or Href. If its relevant, I'm setting the outlet to be zeroGradient for each inital value.

thermoType
type hePsiThermo;
mixture pureMixture;
transport const;
thermo hConst;
equationOfState perfectGas;
specie specie;
energy sensibleInternalEnergy;

mixture
specie

nMoles 1;
molWeight 28.96;

thermodynamics

Cp 1004.5;
Hf 0; // 2.544e+06;

transport

mu 18e-6;
Pr .71;

I appreciate your help here.

chemak · January 23, 2024, 04:30

You can add Tref and Href in this way:

thermodynamics
{
Cp 1004.5;
Hf 0;
Tref 250;
Href 250750;
}

Let me know if it solves your problem.

chc · January 23, 2024, 12:57

I see what you mean. These unfortunately don't change the overall solution. I did simplify the geometry, as I initially had a bump at the outlet location to emulate a the start of a compressor which caused the propagation of a shock upwind. It eventually spills out of my main inlet, at which point the error occurred. When I remove that bump, everything runs smoothly for now, although it takes extremely long. I need to run a type of mesh study to see what resolution is needed.

chc · February 25, 2024, 23:51

For anyone still searching for this, I still need to validate the results but what I wound up doing was using the schemes presented in this tutorial:
http://www.wolfdynamics.com/training..._tutorials.pdf
It's only for laminar flow so you still need to include your turbulent schemes if you don't already have a default to a specific one. I also found that there is a delicate balance between having a small enough maxCellSize and a large enough minCellSize so the results don't blow up or just take forever to complete and it took some time playing around with those to get a simulation to run to completion.
I can post my files here later once I figure out how to do the code format thing.

chc · February 26, 2024, 00:46

U:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];  // kg m s K mol A cd

internalField   uniform (0 0 0);

boundaryField
{
    inlet1
    {
        type            fixedValue;
        value           uniform (680 0 0); // Mach 2 at sea Level

    }

    outlet1
    {
        type            zeroGradient;

    }

    wall1
    {
        type            noSlip;
    }

    top
    {
        type            supersonicFreestream;
        pInf            101325;
        TInf            300;
        UInf            (680 0 0);
        gamma           1.4;
        value           uniform (680 0 0);
    };
    

    bottomEmptyFaces
    {
        type            empty;
    }
    topEmptyFaces
    {
        type            empty;
    }
}

// ************************************************************************* //

T:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];  // kg m s K mol A cd

internalField   uniform 300;

boundaryField
{
    inlet1
    {
        type            fixedValue;
        value           uniform 300;

    }

    outlet1
    {
        type            zeroGradient;

    }

    wall1
    {
        type            zeroGradient;
    }

    top
    {
        type            inletOutlet;
        inletValue      uniform 300;
        value           uniform 300;
    }

    bottomEmptyFaces
    {
        type            empty;
    }
    topEmptyFaces
    {
        type            empty;
    }
}


// ************************************************************************* //

p:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0]; // kg m s K mol A cd

internalField   uniform 101325; // initially atmospheric

boundaryField
{
    inlet1
    {
        type            zeroGradient;
    }

    outlet1
    {
        type            zeroGradient;

    }

    wall1
    {
        type            zeroGradient; // As a result of noSlip used for velocity
    }

    top
    {
        type            zeroGradient;
    }

    bottomEmptyFaces
    {
        type            empty;
    }
    topEmptyFaces
    {
        type            empty;
    }
}


// ************************************************************************* //

omega:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 -1 0 0 0 0]; // kg m s K mol A cd
// calculated from https://www.cfd-online.com/Tools/turbulence.php
internalField   uniform 8; // epsilon/Cmu/k 

boundaryField
{
    "(wall1)"
    {
        type            omegaWallFunction;
        value           $internalField;
    }

    "(inlet1|outlet1|top)"
    {
        type            zeroGradient;
    }

    "(bottomEmptyFaces|topEmptyFaces)"
    {
        type            empty;
    }
}


// ************************************************************************* //

nut:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    "(inlet1|outlet1|top)"
    {
        type            calculated;
        value           uniform 0;
    }

    wall1
    {
        type            nutkWallFunction;
        type            nutkWallFunction;
        value           uniform 0;
    }

    "(topEmptyFaces|bottomEmptyFaces)"
    {
        type            empty;
    }
}

// ************************************************************************* //

k:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];
// calculated from https://www.cfd-online.com/Tools/turbulence.php
internalField   uniform .007; // 

boundaryField
{
    wall1
    {
        type            fixedValue;
        value           uniform .007;
    }

    "(inlet1|outlet1|top)"
    {
        type            zeroGradient;
    }

    "(topEmptyFaces|bottomEmptyFaces)"
    {
        type            empty;
    }
}


// ************************************************************************* //

alphat:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      alphat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    "(inlet1|outlet1|top).*"
    {
        type            calculated;
        value           uniform 0;
    }

    "(bottomEmptyFaces|topEmptyFaces).*"
    {
        type            empty;
    }

    "wall1.*"
    {
        type            compressible::alphatWallFunction;
        Prt             0.85;
        value           uniform 0;
    }
}


// ************************************************************************* //

turbulenceProperties:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType RAS;

RAS
{
    RASModel        kOmegaSST; // test more models
    turbulence      on;
    printCoeffs     on;
}

// ************************************************************************* //

thermophysicalProperties:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            hePsiThermo; // Specifies the type of thermophysical model to be used.
                                 // In this case, it's hePsiThermo. This model is a combination 
                                 // of the h (enthalpy) and psi (specific internal energy) models, 
                                 //and it's often used for compressible flows.
    mixture         pureMixture; // Defines the mixture type. pureMixture indicates that the fluid
                                 // is a pure substance (not a mixture of different species).
    transport       const;       // Specifies the type of transport model to be used. 
                                 // Here, it's set to const, indicating constant transport properties.
    thermo          hConst;      // Indicates the thermo model used for determining thermodynamic 
                                 // properties. hConst specifies that the specific enthalpy (h) will be held constant.
    equationOfState perfectGas;  // Specifies the equation of state to be used. In this case, 
                                 // it's set to perfectGas, indicating that the fluid behaves as an ideal gas.
    specie          specie;      // Indicates the species model. Here, it's set to specie, which
                                 // implies a single-species model (no chemical reactions are considered).
    energy          sensibleInternalEnergy; // Specifies the energy model. sensibleInternalEnergy indicates 
                                            // that the sensible internal energy will be used as the energy variable.
}

// as a tendency: In the source code, typically upper case letters 
// (e.g. Ha) indicate [J/kg] and lower case (e.g. ha) [J/kmol].

mixture

{
    // normalised gas
    specie
    {
        nMoles          1;
        molWeight       28.96; // molar weight of air kg/kmol 
    }
    thermodynamics
    {
        Cp              1004.5; // specific heat J/kg K
        Hf                0; // specific enthalpy
fluid type
    }
    transport
    {
        mu              18e-6; // dynamic viscosity - N s/m^2 or kg/m s, rhoCentralFoam is 0 for inviscid
        Pr              0.71; // prandtl number
        Tref            300;
        Href            250750;
    }
}


// ************************************************************************* //

controlDict:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     rhoCentralFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.5;

deltaT          1e-8;

writeControl    adjustable;

writeInterval   .01;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

graphFormat     raw;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo           0.15;

maxDeltaT       1;


// ************************************************************************* //

fvSolutions:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "(rho|rhoU|rhoE)"
    {
        solver          diagonal;
    }

    "(rho|U|h|k|e|omega)"
    {
        solver          GAMG;
        smoother        GaussSeidel;
        nSweeps         10;
        tolerance       1e-06;
        relTol          0.01;
    }



/*"(U|p|k|nuTilda|omega)"
    {
        solver          GAMG;
        smoother        GaussSeidel;
        nSweeps         10;
        tolerance       1e-6;
        relTol          0.1;
    }
}





PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors 2;
    nCorrectors      2;
    nNonOrthogonalCorrectors 2;
    transonic   yes;
   // rhoMax    10;
}

fvSchemes:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

fluxScheme          Kurganov;

ddtSchemes
{
    default         Euler; 

}

gradSchemes
{
    default         cellLimited Gauss linear 0.333 ;

}

divSchemes
{
    default         Gauss linear;
    div(phi,U) 	    Gauss limitedLinearV 1;
    div(phi,nuTilda)  Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;

    // turbulence      Gauss linearUpwind limited; 
    // div(phi,k)      $turbulence;
    // div(phi,omega)  $turbulence;

}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;

    reconstruct(rho) vanLeer;
    reconstruct(U)  vanLeerV;
    reconstruct(T)  vanLeer;

   }

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method  meshWave;
}

// ************************************************************************* //

meshDict:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      meshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// used in conjunction with the cfMesh, this case, cartesian2DMesh

surfaceFile "Combined.fms"

// minCellSize .0001;
minCellSize .0001;


// maxCellSize .0005;
maxCellSize .006;



// boundaryCellSize .0001;

// ************************************************************************* //

If anyone has any comments, let me know your opinions.

August 29, 2018, 09:06	temperature negative occurred in rhoCentralFoam at high M	#1
longming New Member FENG LM Join Date: Aug 2018 Posts: 1 Rep Power: 0	Hi,guys When I use rhocentralFoam to solve flow over wedge or other simple shape( using RAS model), negative temperature occurred easily when M>5, even dt and maxco are set to very small value. I followed the example in tutorial, it works when computing grid is course. I only refined the grid in boundary layer and the problem just occurred. Is the slover unsteady at high M? Who knows how to slove the problem? Please Help

October 23, 2018, 16:40		#4
TheMadHungarian New Member Andy S. Join Date: Jun 2018 Posts: 17 Rep Power: 8	Usually setting deltaT to a very small number (1e-15 for the 1st time step) will fix this, the solver will adjust the timestep to a larger value automatically. hogsonik and swainwc like this.

January 22, 2024, 04:51		#6
chemak New Member Join Date: Nov 2017 Posts: 6 Rep Power: 9	Hi, I think that the cause of negative temerature in rhoCentralFoam is rather physical than numerical. Check your ground energy in "thermophysicalProperties" (Cp, Tref and Href), setting more physical values will generally prevent this problem. hogsonik likes this.

January 23, 2024, 12:57		#9
chc New Member Join Date: Nov 2023 Posts: 18 Rep Power: 3	I see what you mean. These unfortunately don't change the overall solution. I did simplify the geometry, as I initially had a bump at the outlet location to emulate a the start of a compressor which caused the propagation of a shock upwind. It eventually spills out of my main inlet, at which point the error occurred. When I remove that bump, everything runs smoothly for now, although it takes extremely long. I need to run a type of mesh study to see what resolution is needed. hogsonik likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Negative initial temperature	Lexe	OpenFOAM Running, Solving & CFD	34	November 15, 2024 09:48
[Other] mesh airfoil NACA0012	anand_30	OpenFOAM Meshing & Mesh Conversion	13	March 7, 2022 18:22
channelFoam for a 3D pipe	AlmostSurelyRob	OpenFOAM	3	June 24, 2011 14:06
[blockMesh] BlockMesh FOAM warning	gaottino	OpenFOAM Meshing & Mesh Conversion	7	July 19, 2010 15:11
[blockMesh] Axisymmetrical mesh	Rasmus Gjesing (Gjesing)	OpenFOAM Meshing & Mesh Conversion	10	April 2, 2007 15:00

October 22, 2018, 22:19		#2
hokhay Member Join Date: Nov 2014 Posts: 92 Rep Power: 12	Hi Longming, I am encountering the same problem. Have you found a solution yet? Jason

January 22, 2024, 21:05		#7
chc New Member Join Date: Nov 2023 Posts: 18 Rep Power: 3	My initial set up uses the settings below for that file. I've taken them from relevant articles that reference sea level values for air. I'm solving using dimensional initial values for my velocity, temperature, and pressure, so I'm not sure where to set Tref or Href. If its relevant, I'm setting the outlet to be zeroGradient for each inital value. thermoType type hePsiThermo; mixture pureMixture; transport const; thermo hConst; equationOfState perfectGas; specie specie; energy sensibleInternalEnergy; mixture specie nMoles 1; molWeight 28.96; thermodynamics Cp 1004.5; Hf 0; // 2.544e+06; transport mu 18e-6; Pr .71; I appreciate your help here.

January 23, 2024, 04:30		#8
chemak New Member Join Date: Nov 2017 Posts: 6 Rep Power: 9	You can add Tref and Href in this way: thermodynamics { Cp 1004.5; Hf 0; Tref 250; Href 250750; } Let me know if it solves your problem.

February 25, 2024, 23:51		#10
chc New Member Join Date: Nov 2023 Posts: 18 Rep Power: 3	For anyone still searching for this, I still need to validate the results but what I wound up doing was using the schemes presented in this tutorial: http://www.wolfdynamics.com/training..._tutorials.pdf It's only for laminar flow so you still need to include your turbulent schemes if you don't already have a default to a specific one. I also found that there is a delicate balance between having a small enough maxCellSize and a large enough minCellSize so the results don't blow up or just take forever to complete and it took some time playing around with those to get a simulation to run to completion. I can post my files here later once I figure out how to do the code format thing.