CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Steady state simpleFOAM crash (https://www.cfd-online.com/Forums/openfoam-solving/219680-steady-state-simplefoam-crash.html)

Galactus August 5, 2019 07:46

Steady state simpleFOAM crash
 
1 Attachment(s)
Hello fellow foamers. I tried to simulate a flow inside a U shaped pipe whose cross section changes from circular to elliptical suddenly. I've attached a rough drawing of the geometry (done on MS Paint :p) for your reference. The input is 0.001 m^3/s and the diameter is roughly 0.028m. The flow is unsteady and turbulent. I used k-epsilon modelling with ddt scheme set to steady state.

Problem: The steady state simulation was crashing no matter what I tried. The time step continuity error kept increasing leading to a crash or the bounding of k and epsilon kept becoming bigger and bigger (eg: Min = -1e30 and Max = 1e30) or it just crashed with floating point exception error. I did try various things:
  1. I copied the various tutorial cases from the OpenFOAM 4.x simpleFOAM tutorial directory. Incorporated changes only in 0 and constant folder leaving the system folder as is and ran the simulations. But all of the simulations crashed.
  2. Suspecting that it may be due to running in parallel, I tried serial run as well but same errors cropped up.
  3. Tried running without turbulence modelling in laminar. Did not work.
  4. I tried various matrix iterative solvers in fvSolutions along with various smoothers and pre-conditioners. The time at the which the solution crashed changes but the end result was the same - CRASH!
  5. I then played around with relaxation factors. I used factors as high as 0.95 and as low as 0.3 for all the field variables (p, U, k & epsilon). I did not go any lower as I felt such a simulation will not be correct.
  6. I also tried various schemes in fvSchemes. Bounded, unbounded and also cell limited schemes were used but did not make any difference. They all crashed.
  7. I ran the simulation in transient condition (funnily it worked perfectly) until the flow was fully developed. Then I used the results to initialise the flow for my steady state problem (copied the necessary files from the last time step folder of the transient case to the zero folder of the steady state case) but that also crashed. I was able to see the field initialised through paraview but the steady state simulation just crashed.
  8. I then tried to change the solver to pimpleFOAM but it also crashed. This is when I gave up.

The geometry was meshed using snappyHexMesh and checkMesh did not show any errors. I finally had to run the steady state simulation on Ansys CFX and submit my project. Now that I some time I want to investigate why OpenFOAM could not solve the problem when Ansys could. Even more frustrating thing is that transient simulation in OpenFOAM worked without any issue and even giving me similar result to Ansys. Please help in figuring out the cause for the issue. I am currently away for the day and hence not able to upload case files. Moreover, I did not create a new folder for every change I made (such a noob :() and was overwriting the case files with new changes everytime. Please let me know what file you guys will be requiring. Thanks!

Galactus August 7, 2019 03:36

I am playing around with case file to see where I'm making a mistake. My files in the 0 folder:

p
Code:

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

internalField  uniform 97.1477;

boundaryField
{
    wall
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value          uniform 97.1477;
    }
}

U
Code:

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

internalField  uniform (0 0 0);

boundaryField
{
    wall
    {
        type            noSlip;
    }
    inlet
    {
        type                    flowRateInletVelocity;
        volumetricFlowRate      constant 0.001;
        value                  uniform (0 0 0);
    }
    outlet
    {
        type            zeroGradient;
    }
}

epsilon
Code:

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

internalField  uniform 0.0635805;

boundaryField
{
    wall
    {
        type            epsilonWallFunction;
        value          $internalField;
    }
    inlet
    {
        type            fixedValue;
        value          uniform 0.0635805;
    }
    outlet
    {
        type            zeroGradient;
    }
}

k
Code:

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

internalField  uniform 0.00992673;

boundaryField
{
    wall
    {
        type            kqRWallFunction;
        value          uniform 0.00992673;
    }
    inlet
    {
        type            fixedValue;
        value          uniform 0.00992673;
    }
    outlet
    {
        type            zeroGradient;
    }
}

nut
Code:

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

internalField  uniform 0;

boundaryField
{
    wall
    {
        type            nutkWallFunction;
        value          uniform 0;
    }
    inlet
    {
        type            calculated;
        value          uniform 0;
    }
    outlet
    {
        type            calculated;
        value          uniform 0;
    }
}

nuTilda
Code:

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

internalField  uniform 0;

boundaryField
{
    wall
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
}


Galactus August 7, 2019 03:43

Forgot to add my system files. :p There are many things commented out and many changes were made overwriting the files.
controlDict
Code:

application    simpleFoam;

startFrom      startTime;

startTime      0;

stopAt          endTime;

endTime        1;

deltaT          1e-03;

writeControl    adjustableRunTime;

writeInterval  0.1;

purgeWrite      0;

writeFormat    ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision  6;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo          1;

fvSchemes
Code:

ddtSchemes
{
    default        steadyState;
}

gradSchemes
{
    default        Gauss linear;
}

divSchemes
{
    default        none;
    div(phi,U)      bounded Gauss upwind;
    div(phi,p)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div(phi,omega)  bounded Gauss upwind;
    div(U)          Gauss linear;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default        Gauss linear corrected;
}

interpolationSchemes
{
    default        linear;
}

snGradSchemes
{
    default        corrected;
}

fvSolutions
Code:

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

    pFinal
    {
        $p;
        tolerance      1e-06;
        relTol          0;
    }

    "(U|k|epsilon|omega)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance      1e-05;
        relTol          0.01;
    }

    "(U|k|epsilon|omega)Final"
    {
        $U;
        tolerance      1e-06;
        relTol          0;
    }

}
SIMPLE
{
    nNonOrthogonalCorrectors 1;
    pRefCell        0;
    pRefValue      0;

    residualControl
    {
        p              1e-5;
        U              1e-5;
    }
}

/*
PIMPLE
{
    nOuterCorrectors 50;
    nCorrectors    2;
    nNonOrthogonalCorrectors 1;

    residualControl
    {
        "(U|k|epsilon|omega)"
        {
                tolerance      1e-5;
                relTol          0;
        }
        p
        {
                tolerance      1e-5;
                relTol          0;
        }
    }
}
*/
relaxationFactors
{
    fields
    {
        p      0.3;
        pFinal  1;
    }
    equations
    {
        "U|k|epsilon|omega"            0.3;
        "(U|k|epsilon|omega)Final"      1;
    }
}


sufjanst August 8, 2019 07:53

You can try to run potentialFoam first (remember to backup your u in 0, because potentialFoam overwrites it). It will do a initialization of the u-field. Afterwards you can run simpleFoam



Second, you can try to use the localEuler ddt Scheme. It is pseudo-transient and more stable than steadyState.


Third, if nothing helps you can just average you transient results.


Your pressure B/C is really low. Is it supposed to be like that?
Maybe the problem is in you thermophysicalProperties. Can you post the as well?

peterhess August 8, 2019 22:13

Change the velocity at the outlet to:

outlet
{
type inletOutlet;
value uniform (0 0 0);
inletValue uniform (0 0 0);
}

inletOutlet is much stable than zeroGradient!

-------------------------

Try to add referencePressure to the velocity field at the inlet.

referencePressure value;

This advice works well if the flow compressible.

I dont know if it works in simpleFoam...

Galactus August 13, 2019 04:25

Quote:

Originally Posted by sufjanst (Post 741509)
You can try to run potentialFoam first (remember to backup your u in 0, because potentialFoam overwrites it). It will do a initialization of the u-field. Afterwards you can run simpleFoam



Second, you can try to use the localEuler ddt Scheme. It is pseudo-transient and more stable than steadyState.


Third, if nothing helps you can just average you transient results.


Your pressure B/C is really low. Is it supposed to be like that?
Maybe the problem is in you thermophysicalProperties. Can you post the as well?

Thanks sufjanst and sorry for the late reply. I was on a short vacation.
1) I did try potentialFoam but it did not initialise for the whole domain, rather it initialised only for the inlet boundary surface. I tried running it again and again but initialising was not happening.

2) Tried the localEuler but got a weird error.
Code:

--> FOAM FATAL ERROR:

    request for volScalarField rDeltaT from objectRegistry region0 failed
    available objects of type volScalarField are

6
(
nut
pPrevIter
k
nu
p
epsilon
)


    From function const Type& Foam::objectRegistry::lookupObject(const Foam::word&) const [with Type = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]
    in file /home/ufifi/OpenFOAM/OpenFOAM-4.x/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 193.

FOAM aborting

I could not figure out the reason for error so far, but I'm working on it.

3) Do you mean that I should take results from various time step and just average?

For pressure, I wanted to use 1 atm at outlet. According to the forums, OpenFOAM uses p = p/rho. So I divided 1 atm/ rho. 101325 Pa/1043 Kg/m3.
I don't have a thermophysical dict. I guess those are used in heat transfer analysis. In constant all I have are transport properties and turbulenceDict.

Galactus August 13, 2019 04:37

Quote:

Originally Posted by peterhess (Post 741582)
Change the velocity at the outlet to:

outlet
{
type inletOutlet;
value uniform (0 0 0);
inletValue uniform (0 0 0);
}

inletOutlet is much stable than zeroGradient!

-------------------------

Try to add referencePressure to the velocity field at the inlet.

referencePressure value;

This advice works well if the flow compressible.

I dont know if it works in simpleFoam...

Danke Peter, Sorry for the late reply. I did try using inletOutlet though with "value $internalField;" earlier and it did not work. However, your implementation also does not stabilise the solution. Adding referencePressure did not change anything the simulation still crashed.

Lookid August 13, 2019 07:01

Hello,

you can try more stuff :

- start by a laminar simulation, if it runs, introduce turbulence at a later timestep
- your div schemes are very diffusive (that's good to make it run at first), but you can try to change your other schemes (see p.47 of these slides http://www.wolfdynamics.com/wiki/fvm_crash_intro.pdf)
- Increase your mass flow rate gradually to the final value, may help.

good luck, let us know how it goes :)

peterhess August 13, 2019 12:51

outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value $internalField;

}

Galactus August 14, 2019 04:09

Quote:

Originally Posted by Lookid (Post 741952)
Hello,

you can try more stuff :

- start by a laminar simulation, if it runs, introduce turbulence at a later timestep
- your div schemes are very diffusive (that's good to make it run at first), but you can try to change your other schemes (see p.47 of these slides http://www.wolfdynamics.com/wiki/fvm_crash_intro.pdf)
- Increase your mass flow rate gradually to the final value, may help.

good luck, let us know how it goes :)

Thanks Lilian. I already tried laminar simulation but it does not work. I changed the schemes to the most diffusive one for all. I could notice an improvement but it also crashed. As for mass flow rate, I tried as small as 1e-12 with no luck.

I'm now suspecting that it may be due to meshing. CheckMesh did not show any errors but I'm going to tighten the mesh quality parameters and try running the simulation again.

Galactus August 14, 2019 04:10

Quote:

Originally Posted by peterhess (Post 741995)
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value $internalField;

}

Danke Peter. That is exactly what I tried. Unfortunately it did not work :(

Galactus August 14, 2019 04:14

I've found something interesting. When I increase the nNonOrthogonalCorrectors, the simulation crashes at a much later timestep. This leads me to thinks that I may have some problem with my mesh. CheckMesh did not show any errors. However, I'm planning on remeshing the geometry with tighter mesh quality parameters. Can someone suggest what should be the values of parameters to be set in meshQualityDict for snappyHexMesh. I'm currently using default values found in OpenFOAM-4.x/etc/caseDicts.

Lookid August 14, 2019 07:06

Can you share your files?


All times are GMT -4. The time now is 20:43.