rhoSimpleFoam vs rhoPimpleFoam, steady state vs transient solver

Oliver Meng · June 16, 2022, 07:23

Hi lovely CFD people,

I have a question about steady state solver and transient solver (rhoSimpleFoam vs. rhoPimpleFaom). My case is a train running in a tunnel, with a velocity of 40 m/s. I run the case both with rhoSimpleFoam and rhoPimpleFoam. Of course, I would prefer rhoSimpleFoam if it converges (residuals under 1e-7 or -8).

My questions are:
a.What could be the reasons, why the rhoSimpleFoam does not converge?
b.Can I make the rhoSimpleFoam converge by manipulating fvSchemes or FvSolutions?

Some more information:
My checkMesh shows that the Mesh is 100% ok.
The tunnel and the ground are slip. And the train is no slip. Velocity bc at inlet and pressure bc at outlet.

U:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2106                                  |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    arch        "LSB;label=32;scalar=64";
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform (40 0 0);

boundaryField
{
    "tub.*"
    {
        type            slip;
        //value           uniform (0 0 0);
    }
    
    pod
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    
    noSlipGroundPod
    {
        type            slip;
        //value           uniform (0 0 0);
    }

    up
    {
        type            zeroGradient;
    }

    sides
    {
        type            zeroGradient;
    }


    front
    {
        type            fixedValue;
        value           $internalField;
    }


    back
    {
        type            zeroGradient;
    }
    
    "proc.*"
    {
        type    processor;
    }
}


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

p:

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 1e5;

boundaryField
{
 #includeEtc "caseDicts/setConstraintTypes"
    
    "tub.*"
    {
         type           zeroGradient;
    }
    
    pod
    {
         type           zeroGradient;
    }

    noSlipGroundPod
    {
         type           zeroGradient;
    }

    up
    {
        type            zeroGradient;
    }

    sides
    {
        type            zeroGradient;
    }

    front
    {
        type            zeroGradient;
    }

    back
    {
        type            fixedValue;
        value           $internalField;
    }

    "proc.*"
    {
        type    processor;
    }
}


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

k:

Code:

/*--------------------------------*- C++ -*----------------------------------*\

\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// k = 3/2*(I*U)^2
// Re = 4.75e5, U = 70, I = 0.01

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

internalField   uniform 0.24;

boundaryField
{
    "tub.*"
    {
        type            fixedValue;
        value           $internalField;
    }
    
    pod
    {
        type            kqRWallFunction;
        value           uniform 0;
    }

    noSlipGroundPod
    {
        type            kqRWallFunction;
        value           uniform 0;
    }

    up
    {
        type            zeroGradient;
    }
    
    sides
    {
        type            zeroGradient;
    }

    front
    {
        type            fixedValue;
        value           $internalField;
    }

    back
    {
        type            zeroGradient;
    }

    "proc.*"
    {
        type    processor;
    }
    #includeEtc "caseDicts/setConstraintTypes"
}


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

omega:

Code:

/*--------------------------------*- C++ -*----------------------------------*\

\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      omega;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


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

internalField   uniform 2.4e03;

boundaryField
{
    "tub.*"
    {
        type            omegaWallFunction;
        value           $internalField;
    }
    
    pod
    {
        type            omegaWallFunction;
        value           $internalField;
    }

    noSlipGroundPod
    {
        type            omegaWallFunction;
        value           $internalField;
    }

    up
    {
        type            zeroGradient;
    }
    sides
    {
        type            zeroGradient;
    }

    front
    {
        type            fixedValue;
        value           $internalField;
    }

    back
    {
        type            zeroGradient;
    }

    "proc.*"
    {
        type    processor;
    }

    #includeEtc "caseDicts/setConstraintTypes"
}


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

T:

Code:

/*--------------------------------*- C++ -*----------------------------------*\

\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 298;

boundaryField
{
    "tub.*"
    {
        type            fixedValue;
        value           $internalField;
    }
    
    pod
    {
        type            fixedValue;
        value           $internalField;
    }

    noSlipGroundPod
    {
        type            fixedValue;
        value           $internalField;
    }

    up
    {
        type            zeroGradient;
    }

    sides
    {
        type            zeroGradient;
    }

    front
    {
        type            fixedValue;
        value           $internalField;
    }

    back
    {
        type            zeroGradient;
    }

    "proc.*"
    {
        type    processor;
    }

    #includeEtc "caseDicts/setConstraintTypes"
}


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

alphat:

Code:

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

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

internalField   uniform 0;

boundaryField
{

    pod
    {
        type            compressible::alphatWallFunction;// calculated; 
        value           uniform 0;
    }
    
    "tub.*"
    {
        type            compressible::alphatWallFunction;// calculated; 
        value           uniform 0;
    }

    noSlipGroundPod
    {
        type            compressible::alphatWallFunction;// calculated; 
        value           uniform 0;
    }

    up
    {
        type            calculated;
        value           uniform 0;
    }
    sides
    {
        type            calculated;
        value           uniform 0;
    }

    front
    {
        type            calculated;
        value           uniform 0;
    }

    back
    {
        type            calculated;
        value           uniform 0;
    }

    "proc.*"
    {
        type    processor;
    }

    #includeEtc "caseDicts/setConstraintTypes"
}


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

nut:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2106                                 |
|   \\  /    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
{

    "tub.*"
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }
    
    pod
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }

    noSlipGroundPod
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }

    up
    {
        type            calculated;
        value           uniform 0;
    }

    sides
    {
        type            calculated;
        value           uniform 0;
    }

    front
    {
        type            calculated;
        value           uniform 0;
    }

    back
    {
        type            calculated;
        value           uniform 0;
    }

    "proc.*"
    {
        type    processor;
    }

    #includeEtc "caseDicts/setConstraintTypes"
}


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

SIMPLE-fvSchemes:

Code:

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

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;

    limited         cellLimited Gauss linear 1;
    grad(U)         $limited;
    grad(k)         $limited;
    grad(omega)     $limited;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss linearUpwind limited;

    energy          Gauss linearUpwind limited;
    div(phi,e)      $energy;
    div(phi,K)      $energy;
    div(phi,Ekp)    $energy;

    turbulence      Gauss upwind;
    div(phi,k)      $turbulence;
    div(phi,omega)  $turbulence;

    div(phiv,p)     Gauss upwind;
    div((phi|interpolate(rho)),p) Gauss upwind;

    div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method          meshWave;
}


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

SIMPLE-fvSolution:

Code:

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

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.01;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    "(rho|U|k|omega|e)"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-7;
        relTol          1e-3;
    }

    "(rho|U|k|omega|e)Final"
    {
        $U;
        relTol          0;
    }
    
    cellDisplacement
    {
    solver          GAMG;
    smoother        GaussSeidel;
    tolerance       1e-7;
    relTol          0.01;
    }
    
}

SIMPLE
{
    residualControl
    {
        p               1e-8;
        U               1e-8;
        "(k|omega|e|nut)"   1e-8;
    }

    nNonOrthogonalCorrectors 1;
    pMinFactor      0.05;
    pMaxFactor      5;
}

PIMPLE
{
    nCorrectors              2;
    nNonOrthogonalCorrectors 1;
    nOuterCorrectors         1;
    pMinFactor      0.1;
    pMaxFactor      2;
}

relaxationFactors
{
    fields
    {
        p               0.25;
        rho             0.01; //0.01
    }
    equations
    {
        U               0.7;
        "(k|omega)"     0.7;
        e               0.01;
    }
}


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

PIMPLE-fvSchemes:

Code:

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

ddtSchemes
{
    default         CrankNicolson 0.9;
}

gradSchemes
{
    default         Gauss linear;

    limited         cellLimited Gauss linear 1;
    grad(U)         $limited;
    grad(k)         $limited;
    grad(omega)     $limited;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss linearUpwind limited;

    energy          Gauss linearUpwind limited;
    div(phi,e)      $energy;
    div(phi,K)      $energy;
    div(phi,Ekp)    $energy;

    turbulence      Gauss upwind;
    div(phi,k)      $turbulence;
    div(phi,omega)  $turbulence;

    div(phiv,p)     Gauss upwind;
    div((phi|interpolate(rho)),p) Gauss upwind;

    div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method          meshWave;
}


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

PIMPLE-fvSolution:

Code:

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

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }

    pFinal
    {
        $p;
        relTol          0.1;
    }

    "(rho|U|k|omega|e)"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-8;
        relTol          1e-3;
    }

    "(rho|U|k|omega|e)Final"
    {
        $U;
        relTol          1e-3;
    }
}

SIMPLE
{
    residualControl
    {
        p               1e-8;
        U               1e-8;
        "(k|omega|e)"   1e-8;
    }

    nNonOrthogonalCorrectors 0;
    pMinFactor      0.01;
    pMaxFactor      5;
}

PIMPLE
{
    nCorrectors              2;
    nNonOrthogonalCorrectors 1;
    nOuterCorrectors         10;
    pMinFactor      0.01;
    pMaxFactor      5;
}

relaxationFactors
{
    fields
    {
        p               0.3;
        rho             0.01; //0.01
    }
    equations
    {
        U               0.7;
        "(k|omega)"     0.7;
        e               0.7;
    }
}


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

Tobermory · June 17, 2022, 14:26

One of the difficulties I have observed is that the OpenFOAM solver schemes are very "clean" compared to commercial codes like StarCCM+ of Fluent, and this means that if you have a fine mesh, then there is very little numerical diffusion to stabilise the simulations. This means that flows over bluff bodies will start to show flow unsteadiness & shed turbulent structures into the wake.

I am not sure if this is relevant for your scenario, since we don't have a picture of your geometry, but have a think - is the solver just picking up inherent flow instability? If so, then you won't get full convergence with a steady solver ... unless you coarsen the mesh, which is not ideal! In these circumstances, you may have to run as a transient and field average over a number of shedding cycles.

Also, looking at your relaxation factors, I see you have tightened rho and e down to almost zero ... I know one of the tutorials has 0.01 for rho, but I find it difficult to believe that this is a good way to run the solver, i.e. basically to not update the rho or e equations ...

Good luck!

joshwilliams · June 18, 2022, 06:43

Quote:

Originally Posted by Tobermory

Also, looking at your relaxation factors, I see you have tightened rho and e down to almost zero ... I know one of the tutorials has 0.01 for rho, but I find it difficult to believe that this is a good way to run the solver, i.e. basically to not update the rho or e equations ...

I agree with Tobermory. It seems like a strange concept to use a solver for rho but not really update rho (otherwise, why not just use simpleFoam?).

dlahaye · June 18, 2022, 07:35

@tobermory: what criterium have you employed to access how "clean" OpenFoam is? Are you aware of any documentation that shows this "cleanness"?

@joshwilliams: the confusion might be that rho is updated (using equation of state) despite zero iterations for the rho-solve (recorded in the logfile).

Thx.

Tobermory · June 18, 2022, 08:00

Domenico - the numerical cleanliness of OF depends on the schemes (divScheme, etc.) that you choose, of course, but in OF you can at least choose to run with very little artificial stability by using a good quality orthogonal mesh with little or no stretching, second order schemes etc. In the commercial codes even the second order convection schemes are typically blended 2nd/upwind schemes, to ensure stability ... which ofc is often more important in an industrial CFD situation, where a slightly less precise solution is better than no solution.

Indeed, after moving from StarCCM+ to OF, I was shocked at how often my simulations would blow up - it was a real trial at the start!

arjun · June 18, 2022, 18:29

Quote:

Originally Posted by Tobermory

Domenico - the numerical cleanliness of OF depends on the schemes (divScheme, etc.) that you choose, of course, but in OF you can at least choose to run with very little artificial stability by using a good quality orthogonal mesh with little or no stretching, second order schemes etc. In the commercial codes even the second order convection schemes are typically blended 2nd/upwind schemes, to ensure stability ... which ofc is often more important in an industrial CFD situation, where a slightly less precise solution is better than no solution.

Indeed, after moving from StarCCM+ to OF, I was shocked at how often my simulations would blow up - it was a real trial at the start!

the idea that if the solver is unstable than it must be more accurate compared to stable solver is completely not true.

It is more so not true when you compare OF against StarCCM. CCM does not add sort of artifical dissipation to make it stable and not only that its descretization is actually more accurate.

I personally have never found a case where starccm and of are compared and starccm came out to be less accurate.

Mostly if you benchmark you will find that starccm does better than OF.

I even have a study where most scheme of OF for convection terms were compared in a validation case and starccm was also one of the solver. Guess what there was hardly any scheme that was as good as starccm and around 40 scheme were compared.

Give me email i will send you the results.

dlahaye · June 19, 2022, 16:10

@tobermory: thank you for sharing your thought on the matter. Much appreciated.

@arjun: thank you as well. Is the concern you share related to the fact the starccm is able to take advantage of polyhedral meshes in a way that OF is not?

arjun · June 20, 2022, 03:05

Quote:

Originally Posted by dlahaye

@tobermory: thank you for sharing your thought on the matter. Much appreciated.

@arjun: thank you as well. Is the concern you share related to the fact the starccm is able to take advantage of polyhedral meshes in a way that OF is not?

I used to be developer for starccm and I am probably only one to benchmark its flow model against manufactured solution and study it.

Later when i am writing Wildkatze solver, I am always benchmarking and validating solvers (last 7 years).

Having seen the code myself i know for sure that there is no artifical dissipation added the stability is due to much more accurate schemes used.

We have done lots of benchmarks sometimes we compare OF too just to know how it does.

For example for this case https://fvus.github.io/wildkatze/ver...itverification Tobi (who posts frequently on openfoam forums here too) was kind enough to test many convection schemes. The results do show that openfoam is more dissipative actually (this case actually is the test of it).

Even last two weeks we are testing a test problem and it turned out OF is most dissipative and starccm is least to the point that it ends up generating lots of noise in solution due to it. Starccm tries to maintain the signal while OF pretty much diffuses it.

dlahaye · June 20, 2022, 04:53

Dear Arjun,

Thank you so much for your very valuable input.

I am very happy to see more details of the arguments raised so as gain a better understanding.

Are you able to share more details on the benchmarks and the results obtained with both solvers?

Thank you. Domenico.

Tobermory · June 20, 2022, 05:30

Dear Arjun - thanks for your input, and make no mistake - I am a BIG fan of starccm+ ... and have used starccm+ and StarCD before that since the early 90s. My earlier comments about stability were not intended as a discussion of the accuracy of starccm+ over any other code, but I stand by the observation that it (and other leading commercial CFD codes) is much more stable than OF - that must surely come from more robust & forgiving numerical schemes, and perhaps also more forgiving boundary condition coding? Perhaps you can enlighten us, with your experience at CD-Adapco.

Thanks also for the link to Tobi's test case. Regarding validation cases, my view is that you need to be VERY careful about mesh effects - the mesh in the test case, 2 cells across a duct, is clearly not sufficient to resolve anything physical, and so my suggestion is that it is difficult to learn anything real from this test case. I remember a 3D validation case I did against wind tunnel data where brilliant agreement was found from a coarse mesh solution, but the agreement got worse the more I refined the mesh .. until I got to the "real" mesh-independent solution ... that's the one I should have been using at the start to judge the numerical schemes.

Anyway - thank you for the discussion.

arjun · June 20, 2022, 05:44

Quote:

Originally Posted by dlahaye

Dear Arjun,

Thank you so much for your very valuable input.

I am very happy to see more details of the arguments raised so as gain a better understanding.

Are you able to share more details on the benchmarks and the results obtained with both solvers?

Thank you. Domenico.

It was long time ago so give me a day or two to find out his files. It can be useful for people who use openfoam.

arjun · June 20, 2022, 05:49

Quote:

Originally Posted by Tobermory

Thanks also for the link to Tobi's test case. Regarding validation cases, my view is that you need to be VERY careful about mesh effects - the mesh in the test case, 2 cells across a duct, is clearly not sufficient to resolve anything physical, .

1. The benchmark solution was generated from extremely fine mesh. So it was very physical. This is actually one of the validation cases of Ansys Fluent.

2. In case you did not notice third order solver with those 2 cells width able to reproduce solution generated by very very fine mesh.

3. All solvers used the same mesh. Starccm and Wildkatze did very good but you see if they were as diffusive as you believe they would have done worse.

4. As I said we have done multiple tests. I share you image from last test. OpenFOAM is most diffusive here and Starccm the least. Wildkatze has balanced solution with least amount of noise thanks to new specialised solver in it.

Tobermory · June 20, 2022, 05:59

Quote:

Originally Posted by arjun

2. In case you did not notice third order solver with those 2 cells width able to reproduce solution generated by very very fine mesh.

This is the bit that puzzles me - why would anyone run a CFD simulation with just two cells across the domain? That fails the basics of Best Practice. I am not really interested in what that solution gives, since as I said before it contains no real physics. Much more interesting is the grid independent solution ... you should only be comparing grid independent solutions when examining the effects of solver parameters.

Oliver Meng · June 20, 2022, 17:55

Thanks for your reply.

A little more information about my case. The geometry is actually like a nozzle, wider at the inlet and outlet, very tight in the middle. And the maximum velocity will lead to the chocking effect.

i have tried different cases with different velocities at the inlet. I successfully made all simulations converge (residuals less then 1e-6)with rhoSimplefoam, as long as the max. velocity didn't reach 1 Ma. But for the cases, which should have a higher max. local velocity than 1 Ma (choked), it still did not work.

So my new question is: maybe a steady state solver like rhoSimpleFoam is not capable of a transonic or supersonic compressible case? Which means i really have to switch to rhoPimpleFoam for example?

Quote:

Originally Posted by Tobermory

One of the difficulties I have observed is that the OpenFOAM solver schemes are very "clean" compared to commercial codes like StarCCM+ of Fluent, and this means that if you have a fine mesh, then there is very little numerical diffusion to stabilise the simulations. This means that flows over bluff bodies will start to show flow unsteadiness & shed turbulent structures into the wake.

I am not sure if this is relevant for your scenario, since we don't have a picture of your geometry, but have a think - is the solver just picking up inherent flow instability? If so, then you won't get full convergence with a steady solver ... unless you coarsen the mesh, which is not ideal! In these circumstances, you may have to run as a transient and field average over a number of shedding cycles.

Also, looking at your relaxation factors, I see you have tightened rho and e down to almost zero ... I know one of the tutorials has 0.01 for rho, but I find it difficult to believe that this is a good way to run the solver, i.e. basically to not update the rho or e equations ...

Good luck!

Tobermory · June 20, 2022, 18:06

Oliver - it is possible to use rhoSimpleFoam for transonic/sonic cases - see the attached results from a nozzle flow at a pressure ratio of 6.8, run with rhoSimpleFoam. I always prefer to use a steady solver where I can, but as I mentioned in an earlier post, it can sometimes be difficult to get a stable solution.

Incidentally, have you turned on the transonic switch in fvSolution/SIMPLE? I also run with the consistent flag set to yes.

Oliver Meng · June 20, 2022, 18:24

My pressure ratio should also be round 6.

This is my new fvSolution file with transonic and consistent turning on.

But the simulation crashed at the second iteration with residuals of p ending up at 1e47

Code:

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

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-10;
        relTol          0.1;
    }


    "(rho|U|k|omega|e)"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-10;
        relTol          1e-3;
    }
    
}

SIMPLE
{
    residualControl
    {
        p               1e-10;
        U               1e-10;
        "(k|omega|e|nut)"   1e-10;
    }
    transonic true;
    consistent true;
    nNonOrthogonalCorrectors 1;
    pMinFactor      0.05;
    pMaxFactor      6;
}

PIMPLE
{
    nCorrectors              2;
    nNonOrthogonalCorrectors 1;
    nOuterCorrectors         1;
    pMinFactor      0.05;
    pMaxFactor      6;
}

relaxationFactors
{
    fields
    {
        p               0.2;
        rho             0.1;
        pEqn            0.3;
    }
    equations
    {
        U               0.7;
        "(k|omega)"     0.7;
        e               0.7;
    }
}


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

Code:

DILUPBiCGStab:  Solving for Ux, Initial residual = 0.02137706153, Final residual = 2.082396725e-05, No Iterations 2
DILUPBiCGStab:  Solving for Uy, Initial residual = 0.1155308207, Final residual = 0.0001081439417, No Iterations 2
DILUPBiCGStab:  Solving for Uz, Initial residual = 0.1310687182, Final residual = 0.0001163895147, No Iterations 2
DILUPBiCGStab:  Solving for e, Initial residual = 0.2772536778, Final residual = 0.0001556240328, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.2028451406, Final residual = 3.920168245e+47, No Iterations 1000
GAMG:  Solving for p, Initial residual = 0.0337139218, Final residual = 1.233059617e+48, No Iterations 1000
time step continuity errors : sum local = 1.154070612e+101, global = -6.316817843e+100, cumulative = -6.316817843e+100
pressureControl: p max 1.383434598e+104
DILUPBiCGStab:  Solving for omega, Initial residual = 0.004262390935, Final residual = 1.678025594e-06, No Iterations 3
DILUPBiCGStab:  Solving for k, Initial residual = 0.03920934066, Final residual = 2.115743412e-05, No Iterations 3
ExecutionTime = 29.61 s  ClockTime = 30 s

Quote:

Originally Posted by Tobermory

Oliver - it is possible to use rhoSimpleFoam for transonic/sonic cases - see the attached results from a nozzle flow at a pressure ratio of 6.8, run with rhoSimpleFoam. I always prefer to use a steady solver where I can, but as I mentioned in an earlier post, it can sometimes be difficult to get a stable solution.

Incidentally, have you turned on the transonic switch in fvSolution/SIMPLE? I also run with the consistent flag set to yes.

arjun · June 21, 2022, 01:48

Quote:

Originally Posted by Tobermory

This is the bit that puzzles me - why would anyone run a CFD simulation with just two cells across the domain? .

To demonstrate that high order solver is more accurate. That was point of that validation no.

Quote:

Originally Posted by Tobermory

I am not really interested in what that solution gives, since as I said before it contains no real physics.

May be this is why you think the way you think since you completely missed that the solution produced was very physical. You just ignore the facts in front of you for the myths you believe in.

Quote:

Originally Posted by Tobermory

Much more interesting is the grid independent solution ... you should only be comparing grid independent solutions when examining the effects of solver parameters.

Not if my higher order solver could get me more accurate results on coarser grids. I am not stupid to keep running simulations for days when higher order would do the same 10 times faster.

BTW the last picture was not some very coarse simulation. That mesh has enough cells and you can see that OpenFOAM was far more diffusive than other two solvers and Starccm was least diffusive. Goes completely against the idea that OpenFOAM is not using any artifical diffusivity while Starccm is using it. Another myth that people keep repeating here.

Tobermory · June 21, 2022, 05:40

Quote:

Originally Posted by Oliver Meng

My pressure ratio should also be round 6.

This is my new fvSolution file with transonic and consistent turning on.

But the simulation crashed at the second iteration with residuals of p ending up at 1e47

Did you apply the full pressure from the start? I needed to ramp up my pressure over several thousand iterations, to give the pressure solver a chance. I used a uniformTotalPressure condition:

Code:

    upStream {
        type            uniformTotalPressure;
        p0              table (
                            (     0   $p0)
                            (  1000   150000)
                            (  2000   300000)
                            (  3000   $pInlet)
                            (999999   $pInlet)
                        );
    }

and that seemed to do the trick for me, with a totalPressure condition on the downstream outley boundary. For relaxation I had:

Code:

    fields
    {
        p               0.3;
        rho             0.2;
    }
    equations
    {
        velocity        0.5;
        pressure        0.9;
        energy          0.5;
        turbulence      0.7;
        U               $velocity;
        p               $pressure;
        "(k|epsilon|omega)" $turbulence;
        "(h|e)"         $energy;
    }

and I also used fvOptions to constrain T in some of my later simulations:

Code:

    limitT
    {
        type            limitTemperature;
        active          yes;
        selectionMode   all;
        min             50;
        max             750;
    }

to help with some of the early initial transients, where the temperature/h would swing wildly. I hope this helps.

Tobermory · June 21, 2022, 05:46

Quote:

Originally Posted by arjun

To demonstrate that high order solver is more accurate. That was point of that validation no.

Thank you for giving me a giggle so early on in the morning. I love the idea that you think that two to three cells are enough to resolve fully the complex flow features, including free shear layers and turbulent entrainment, that are required to reproduce the physics involved in your test case.

We are not going to agree on this, so perhaps let us just focus on Oliver's original question, my friend. All the best.

Oliver Meng · June 21, 2022, 07:49

Thanks you so much! I am trying these methods from you.

Where is exactly your upStream patch in the geometry?

Quote:

Originally Posted by Tobermory

Did you apply the full pressure from the start? I needed to ramp up my pressure over several thousand iterations, to give the pressure solver a chance. I used a uniformTotalPressure condition:

Code:

    upStream {
        type            uniformTotalPressure;
        p0              table (
                            (     0   $p0)
                            (  1000   150000)
                            (  2000   300000)
                            (  3000   $pInlet)
                            (999999   $pInlet)
                        );
    }

and that seemed to do the trick for me, with a totalPressure condition on the downstream outley boundary. For relaxation I had:

Code:

    fields
    {
        p               0.3;
        rho             0.2;
    }
    equations
    {
        velocity        0.5;
        pressure        0.9;
        energy          0.5;
        turbulence      0.7;
        U               $velocity;
        p               $pressure;
        "(k|epsilon|omega)" $turbulence;
        "(h|e)"         $energy;
    }

and I also used fvOptions to constrain T in some of my later simulations:

Code:

    limitT
    {
        type            limitTemperature;
        active          yes;
        selectionMode   all;
        min             50;
        max             750;
    }

to help with some of the early initial transients, where the temperature/h would swing wildly. I hope this helps.

June 17, 2022, 14:26		#2
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 736 Rep Power: 14	One of the difficulties I have observed is that the OpenFOAM solver schemes are very "clean" compared to commercial codes like StarCCM+ of Fluent, and this means that if you have a fine mesh, then there is very little numerical diffusion to stabilise the simulations. This means that flows over bluff bodies will start to show flow unsteadiness & shed turbulent structures into the wake. I am not sure if this is relevant for your scenario, since we don't have a picture of your geometry, but have a think - is the solver just picking up inherent flow instability? If so, then you won't get full convergence with a steady solver ... unless you coarsen the mesh, which is not ideal! In these circumstances, you may have to run as a transient and field average over a number of shedding cycles. Also, looking at your relaxation factors, I see you have tightened rho and e down to almost zero ... I know one of the tutorials has 0.01 for rho, but I find it difficult to believe that this is a good way to run the solver, i.e. basically to not update the rho or e equations ... Good luck! dlahaye, hogsonik and joshwilliams like this.

June 18, 2022, 07:35		#4
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 798 Blog Entries: 1 Rep Power: 17	@tobermory: what criterium have you employed to access how "clean" OpenFoam is? Are you aware of any documentation that shows this "cleanness"? @joshwilliams: the confusion might be that rho is updated (using equation of state) despite zero iterations for the rho-solve (recorded in the logfile). Thx. hogsonik and Tobermory like this.

June 18, 2022, 08:00		#5
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 736 Rep Power: 14	Domenico - the numerical cleanliness of OF depends on the schemes (divScheme, etc.) that you choose, of course, but in OF you can at least choose to run with very little artificial stability by using a good quality orthogonal mesh with little or no stretching, second order schemes etc. In the commercial codes even the second order convection schemes are typically blended 2nd/upwind schemes, to ensure stability ... which ofc is often more important in an industrial CFD situation, where a slightly less precise solution is better than no solution. Indeed, after moving from StarCCM+ to OF, I was shocked at how often my simulations would blow up - it was a real trial at the start! dlahaye likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
From steady state to transient simulation	Mechand	FLUENT	0	December 24, 2020 13:17
Will the results of steady state solver and transient solver be same?	carye	OpenFOAM Running, Solving & CFD	9	December 28, 2019 06:21
Effect of initial condition for steady state vs Transient	prasa	ANSYS	0	August 22, 2018 05:45
accelarate simulation using steady state solver	tjliang	OpenFOAM Running, Solving & CFD	0	October 7, 2016 15:42
Constant velocity of the material	Sas	CFX	15	July 13, 2010 09:56

June 19, 2022, 16:10		#7
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 798 Blog Entries: 1 Rep Power: 17	@tobermory: thank you for sharing your thought on the matter. Much appreciated. @arjun: thank you as well. Is the concern you share related to the fact the starccm is able to take advantage of polyhedral meshes in a way that OF is not?

June 20, 2022, 04:53		#9
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 798 Blog Entries: 1 Rep Power: 17	Dear Arjun, Thank you so much for your very valuable input. I am very happy to see more details of the arguments raised so as gain a better understanding. Are you able to share more details on the benchmarks and the results obtained with both solvers? Thank you. Domenico.

June 20, 2022, 05:30		#10
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 736 Rep Power: 14	Dear Arjun - thanks for your input, and make no mistake - I am a BIG fan of starccm+ ... and have used starccm+ and StarCD before that since the early 90s. My earlier comments about stability were not intended as a discussion of the accuracy of starccm+ over any other code, but I stand by the observation that it (and other leading commercial CFD codes) is much more stable than OF - that must surely come from more robust & forgiving numerical schemes, and perhaps also more forgiving boundary condition coding? Perhaps you can enlighten us, with your experience at CD-Adapco. Thanks also for the link to Tobi's test case. Regarding validation cases, my view is that you need to be VERY careful about mesh effects - the mesh in the test case, 2 cells across a duct, is clearly not sufficient to resolve anything physical, and so my suggestion is that it is difficult to learn anything real from this test case. I remember a 3D validation case I did against wind tunnel data where brilliant agreement was found from a coarse mesh solution, but the agreement got worse the more I refined the mesh .. until I got to the "real" mesh-independent solution ... that's the one I should have been using at the start to judge the numerical schemes. Anyway - thank you for the discussion.