[AndreasP.]

Hello everybody,

I have been trying for quite some time to run a compressible simulation of a NACA 0012 profile, using rhoSimplecFoam and kOmegaSST for turbulence modeling.
My mesh is unstructured, checkMesh indicates no problems.
Running a simulation at Mach 0.7 generates following output:

Starting time loop

Time = 0.1

GAMG: Solving for Ux, Initial residual = 1, Final residual = 0.0615838, No Iterations 1
GAMG: Solving for Uy, Initial residual = 1, Final residual = 0.0598451, No Iterations 1
GAMG: Solving for p, Initial residual = 1, Final residual = 27.4647, No Iterations 1000
time step continuity errors : sum local = 0.793546, global = 0.47889, cumulative = 0.47889
rho max/min : 0.1 0.1
GAMG: Solving for h, Initial residual = 0.627704, Final residual = 0.906171, No Iterations 1000
GAMG: Solving for omega, Initial residual = 0.00118687, Final residual = 5.36443e-06, No Iterations 1
smoothSolver: Solving for k, Initial residual = 1, Final residual = 0.0494026, No Iterations 2
ExecutionTime = 7.73 s ClockTime = 7 s

Time = 0.2

GAMG: Solving for Ux, Initial residual = 0.335287, Final residual = 0.00814846, No Iterations 2
GAMG: Solving for Uy, Initial residual = 0.0892067, Final residual = 0.00223789, No Iterations 2
#0 Foam::error:rintStack(Foam::Ostream&) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2 in "/lib/libc.so.6"
(...
...)
in "/opt/openfoam201/platforms/linux64GccDPOpt/bin/rhoSimplecFoam"
#10 __libc_start_main in "/lib/libc.so.6"
#11
in "/opt/openfoam201/platforms/linux64GccDPOpt/bin/rhoSimplecFoam"
Floating point exception

I presume there's something off when calculating for p and h. I have no clue what the problem is though, an incompressible simulation, using simpleFoam and Spallart-Allmaras works fine.

Switching transonic off in fvSolution seems to stabilize the simulate at least so much to allow a minimal amount of postprocessing; every single variable diverges massively. Changing the relaxation factors has no significant effect.

I am quite new to OpenFoam, and have generated my files according to various threats in this forum as well as to the tutorial cases.
Here are my files fvSolution, fvSchemes, p and U.

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

solvers
{
    
    p0
    {
        solver          GAMG;
        tolerance       1e-10;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }
    p1
    {
        solver          GAMG;
        tolerance       1e-10;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }
    p
    {
        solver          GAMG;
        tolerance       1e-10;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }
    omega0
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          0.1;
    }

    omega1
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         1;
        tolerance       1e-10;
        relTol          0.1;
    }

    omega
    {
        solver          GAMG;
        tolerance       1e-10;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    U0
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          0.1;
    }

    U1
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         1;
        tolerance       1e-10;
        relTol          0.1;
    }

    U
    {
        solver          GAMG;
        tolerance       1e-10;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    h0
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          0.1;
    }

    h1
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         1;
        tolerance       1e-10;
        relTol          0.1;
    }

    h
    {
        solver          GAMG;
        tolerance       1e-10;
        relTol          0.1;
        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    k0
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          0.1;
    }

    k1
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-10;    
        relTol          0.1;
    }

    k
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-10;
        relTol           0.1;
        nSweeps          1;
    }

    

   
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    rhoMin          rhoMin [1 -3 0 0 0] 0.1;
    rhoMax          rhoMax [1 -3 0 0 0] 1.0;
    transonic       yes;                                      
    pRefPoint        (100 0 0);
    pRefValue       11000;
}

relaxationFactors
{
    p               0.5; 
    U               0.7; 
    k               0.5; 
    omega           0.5;            
    rho            1;                       
    h            0.95;
}


relaxationFactors0
{
    p             0.3;
    rho           0.4;                
    U             0.7;
    h             0.7;
    k             0.7;
    omega         0.7;
}

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

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

ddtSchemes
{
    default             steadyState;
}

gradSchemes
{
    default         fourth;
    grad(p)         fourth;
    grad(U)         Gauss linear;
    grad(omega)        Gauss linear;
}

divSchemes
{
    default             none;

    div(phi,U)          Gauss upwind; //limitedLinearV 1; //upwind;
    div((muEff*dev2(T(grad(U)))))      Gauss linear;
    div(phi,h)          Gauss upwind; //limitedLinear 1; //upwind;

    div(phi,k)          Gauss upwind;
    div(phi,omega)    Gauss upwind;
    div(phid,p)         Gauss upwind;
    div(U,p)            Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited 0.5;
}

interpolationSchemes
{
    default                 none;
    interpolate(rho)        linear;
    div(U,p)                upwind phi;
    interpolate((psi*U))    linear;
    interpolate(U)          linear;
    UD                      upwind phid;
    interpolate(p)          linear;
    interpolate(((rho|(A(U)-H(1)))-(rho|A(U)))) linear;
    interpolate((rho*U))    linear;
}

snGradSchemes
{
    default                 corrected;
}

fluxRequired
{
    default                 no;
    p;
    pCorr;
}

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

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

internalField   uniform 11000;

boundaryField
{
    NONE1      //Nearfield1
    {
        type            empty;
        
    }

    NONE2    //Nearfield2
    {
        type            empty;
    }
 
    NONE3     //Pressure Side
    {
        type            zeroGradient;
    }

    NONE4    //Nearfield3
    {
        type            empty;
    }
    
    NONE5    //Nearfield4
    {
        type            empty;
    }

    NONE6    //Suction Side
    {
        type            zeroGradient;
    }

    NONE7    //Wake resolution1
    {
        type            empty;
    }

    NONE8    //Wake resolution2
    {
        type            empty;
    }
    
    NONE9    //Trailing Edge
    {
        type            zeroGradient;
    }

    NONE10    //Farfield1
    {
        type            empty;
    
    }
    NONE11    //Farfield2
    {
        type            empty;
      
    }
    
    NONE12    //Outer rim, lower half
    {
        type            freestreamPressure;    
        //type            fixedValue;
        value           uniform 11000;
    }
    
    NONE13    //Outer rim, upper half
    {
        type            freestreamPressure;
        //type            fixedValue;
        value           uniform 11000;
    }
}

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

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

internalField   uniform (239.266 20.933 0);


boundaryField
{
    NONE1      //Nearfield1
    {
        type            empty;
        
    }

    NONE2    //Nearfield2
    {
        type            empty;
    }
 
    NONE3     //Pressure Side
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    NONE4    //Nearfield3
    {
        type            empty;
    }
    
    NONE5    //Nearfield4
    {
        type            empty;
    }

    NONE6    //Suction Side
    {
       type            fixedValue;
       value           uniform (0 0 0);
    }

    NONE7    //Wake resolution1
    {
        type            empty;
    }

    NONE8    //Wake resolution2
    {
        type            empty;
    }
    
    NONE9    //Trailing Edge
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    NONE10    //Farfield1
    {
        type            empty;
    
    }
    NONE11    //Farfield2
    {
        type            empty;
      
    }
    
    NONE12    //Outer rim, lower half
    {
        type            freestream;
        freestreamValue uniform (239.266 20.933 0);   
    }
    
    NONE13    //Outer rim, upper half
    {
       type            freestream;
       freestreamValue uniform (239.266 20.933 0);
    }
}

Thanks in advance for your help.

[AndreasP.]

*Bump*

Any help would be greatly appreciated

[hannes]

Originally Posted by AndreasP.

Hello everybody,
rho max/min : 0.1 0.1

Looks like you always hit the lower limit of rho (given in fvSolution/SIMPLE-section).
I think you should adapt "rhoMin" to your quite low abient pressure of 11kPa.

[AndreasP.]

Hey Hannes,
thanks for the reply!
Unfortunately, decreasing rhoMin seems to have no effect...

[fredo490]

It's a bit late to answer but I often see thread with problems over rhoMin and rhoMax.

Those two values are some "limiters" to avoid unphysical values of rho. This limiter has to be temporary (often at the beginning of the computation to avoid any divergence).

You need to know first the freestream value of rho. If you are lazy you can use Wolfram: http://www.wolframalpha.com/input/?i...=&equal=Submit

Then you have to choose your rhoMin and rhoMax arround the freestream value. And here you need to be smart and use your experience because the min and max often depend of the geometry. Until Mach 0.5 the density usually change of less than +/- 0.2.