[K-al-Eps-o]

Hi,

I have been trying to simulate a very simple nozzle flow with OpenFOAM flow in the past weeks. I have tried many different parameters and hypotheses, changing them only one at a time, but everyone one them has led to the same result: the simulation crashes with the error message "negative temperature inlet T0=-XXXK".

The domain I am considering is a short straight section followed by a simple convergent-divergent with straight walls sections. To take advantage of the geometry, I am solving the flow with axisymmetric BC. My domain is then a simple wedge - of less than 5°.

The domain is meshed with gmsh, creating a structured hexahedral mesh with grading from the center to the walls (growth ratio < 0.95). A grading is also applied from the inlet to the throat; after the throat, the mesh size is constant.

rhoCentralFoam is used as a density-based solver and different turbulent models have been tried - k-eps, RNG k-eps, k-omega-SST. The courant number is set to Co<0.2.

Boundary conditions are set as follow:

• inlet:
  p: total pressure 9.1MPa
  U: InletOutlet
  T: total temperature 305K
  eps: zeroGradient
  omega: zeroGradient
  k: zeroGradient
  alphaT: calculated
• outlet:
  p: wave transmissive 0.101MPa
  U: InletOutlet
  T: zeroGradient
  eps: zeroGradient
  omega: zeroGradient
  k: zeroGradient
  alphaT: calculated
• walls:
  p: zeroGradient
  U: noSlip
  T: zeroGradient
  eps: epsilonWallFunction
  omega: kqRWallFunction
  k: kqRWallFunction
  alphaT: compressible::alphatWallFunction
• wedges(lateral faces of the wedge):
  wedge

I am using CO2 in supercritical conditions as a fluid. I want to use a perfect gas model in the first place. Thus, my thermophysical properties as are follow:

thermoType
{
    type            hePsiThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    // normalised gas
    specie
    {
	nMoles		1;
        molWeight       44.01;
    }
    thermodynamics
    {
        Cp              2327;
        Hf              0;
    }
    transport
    {
	mu		1.444e-5;
	Pr		1.5;
    }
}

For controlDisct, writeFormat is set to binary and writePrecision to 10 because I had issues with the wedge normals being misaligned; maxCo is set to 0.2 and maxDeltaT to 0.0001.

My fvSchemes is set as follow:

fluxScheme          Kurganov;

ddtSchemes
{
    default             Euler;
}

gradSchemes
{
    default             Gauss linear;

    limited             cellLimited Gauss linear 1;
}

divSchemes
{
    default             none;

    div(phi,U)          Gauss linearUpwind limited;
    div(phi,e)          Gauss linearUpwind limited;
    div(phi,K)          Gauss linearUpwind limited;

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

    div(phiv,p)         Gauss linearUpwind limited;
    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;

    div(tauMC)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;

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

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method          meshWave;
}

and my fvSolution as follow:

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-08;
        relTol          0.01;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    "rho.*"
    {
        $p;
        tolerance       1e-05;
        relTol          0;
    }


    "(omega|U|e|h|k|epsilon).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-08;
        relTol          0.01;
    }
}

The problem is solved in parallel using mpi.

Parameters that I tried to modify:

• play with the discretization schemes
• change the initialization value of p and U
• initialize k, eps, omega to high/low values to stabilize the solution
• finner / coarser mesh
• different geometries : change the throat and outlet diameters
• turbulent models: k-eps, RNG k-eps, K-omega, K-omega-SST
• change the turbulent quantities - k, eps, omega - at inlet (instead of zeroGradient BC)
• zeroGradient BC for p at outlet

note: When solving in 2D - 3D mesh, 1 cell wide with empty BC - the simulation converges with k-omegaSST turbulence model. All the other models make the simulation diverge (same error).

I am open to any discussion about my simulation. Please don't hesitate to ask for more information.

Best.

[K-al-Eps-o]

In addition, please find attached pictures of the mesh I am using.

I also want to mention that I tried tu implement limitTemperature with fvOptions. However, it led me to the same "negative inlet temperature" error.

[K-al-Eps-o]

By any chance, does anyone have any idea or insight ?

Best

[K-al-Eps-o]

Is there a way to get people to see this post ? I know that multiple post for the same probleme is not allowed.

Best.

[kinematic_presser]

OpenFOAM is not so good at importing meshes from other softwares. The negative temp at the inlet could be caused by misaligned normals. Check in gmsh that your normals are all pointing the correct direction.

[K-al-Eps-o]

Originally Posted by kinematic_presser

OpenFOAM is not so good at importing meshes from other softwares. The negative temp at the inlet could be caused by misaligned normals. Check in gmsh that your normals are all pointing the correct direction.

You think this issue could come from bad normals ? The fact that the simulation runs for a long time before it crashes doesn't exclude this hypothesis ?

[JBeilke]

Just have a look at this tutorial:
https://github.com/unicfdlab/hybridC...ivergingNozzle

It uses pimpleCentralFoam but might help.

[K-al-Eps-o]

If anyone encounter a similar solution: I solved it by using setFields to initialize de fields in the domain.

The issue was that with the domain being initialized at 1bar and the inlet being at 91bar, it creates a shock in the convergent right at the beginning of the simulation. This shock propagates but unfortunately cannot be damped and creates a divergence when reaching the throat. By simply initializing the domain to highest pressures, and with a jump of pressure AFTER the throat, the simulation converged.

[LuckyTran]

I want to go back to something

Originally Posted by K-al-Eps-o

Parameters that I tried to modify:

• play with the discretization schemes
• change the initialization value of p and U
• initialize k, eps, omega to high/low values to stabilize the solution
• finner / coarser mesh
• different geometries : change the throat and outlet diameters
• turbulent models: k-eps, RNG k-eps, K-omega, K-omega-SST
• change the turbulent quantities - k, eps, omega - at inlet (instead of zeroGradient BC)
• zeroGradient BC for p at outlet

You can imagine it is very hard for anyone to suggest a recommendation when it is already on the list of exclusions. It's like asking for recommendations to lose weight but "I've already tried changing my diet and exercising more and it didn't work." Well gee, I guess I forgot to mention that you need to decrease calorie intake and increase calorie consumption and every recommendation suddenly starts sounding like "just gitgud :4head:"


However, I do applaud you for updating your own post and confirming that initialization errors are still a common occurrence in 2023 and will likely continue to be a pain point for CFDers for the foreseeable future.

[K-al-Eps-o]

Originally Posted by LuckyTran

I want to go back to something



You can imagine it is very hard for anyone to suggest a recommendation when it is already on the list of exclusions. It's like asking for recommendations to lose weight but "I've already tried changing my diet and exercising more and it didn't work." Well gee, I guess I forgot to mention that you need to decrease calorie intake and increase calorie consumption and every recommendation suddenly starts sounding like "just gitgud :4head:"


However, I do applaud you for updating your own post and confirming that initialization errors are still a common occurrence in 2023 and will likely continue to be a pain point for CFDers for the foreseeable future.

I can feel the sarcasm...