[Galactus]

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 ) 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]

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

p

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

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

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

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

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

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

internalField   uniform 0;

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

[Galactus]

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

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

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

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]

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]

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]

Originally Posted by sufjanst

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.

--> 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]

Originally Posted by peterhess

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]

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]

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

}

[Galactus]

Originally Posted by Lookid

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]

Originally Posted by peterhess

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

}

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

[Galactus]

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]

Can you share your files?