[mgdenno]

Hello All,

I am hoping someone here can provide some guidance for a simulation that I am trying to run. I am trying to analyze a Labyrinth Spillway using the interFoam solver. Attached are a few screen captures to give an idea of the overall setup.

I have been able to run approximately 21.7 seconds of simulation before it blows up (I don’t know if it matters but I stopped it after 19 seconds to reconstructPar and view it, and then restarted it from 19 seconds). My fvSolution and fvSchemes files are the same as those used for the dam break model except for nNonOrthogonalCorrectors = 1 (I have also tried 2, 3 and 4 with no effect). I have attached the log file starting from 21 seconds as well as my 0 directory with my boundary conditions in it. The deltaT gets very small and eventually the simulation stops running. It looks to me that the time step effectively goes to 0 and causes a floating point exception, but I am not sure how to pin point the cause (bad mesh, boundary condition, etc.). I suspect that it is related to the k and/or epsilon but I am not sure how to correct it. I have run checkMesh and it returns OK.

I would be very appreciative if someone could point me in the right direction to get this model running.

Regards,

MD

[pablodecastillo]

did you try with calpha=0.5 or less?

[mgdenno]

No, I haven't tried that, but I will now.

Thanks,

MD

[mgdenno]

Pablo,

Thanks again for the idea, unfortunately it didn't fix the problem.

I included a few plots of the output; maybe they will mean more to you than they do to me. It seems to run fine for a few tenths of a second before going unstable then eventually quits.

Any other suggestions?

Thanks,

MD

[kumar]

I would suggest you try to change your pressure boundary condition. If yuo give me some details about your pressure and velocity boundary conditions, may be I can help.

Do you have an Inlet and Outlet for your domain. It is likely that if your definition of pressure at the outlet is zerogradient then it is possible to have this kind of problems.


regards
K.Suresh kumar

[pablodecastillo]

Like Suresh is pointing, usually i got better results with totalpressure at the outlet and not zerogradient.

Pablo

[mgdenno]

Thank you both of you for your responses.

I am using inletOutlet for the outlet, so yes effectively zeroGradient for flows out of the domain. My p_rgh and U files are below.

I tried to look at using pressure inlets and outlets previously, but couldn't make them work (I am pretty new to CFD) - I would be very interested to hear your suggestions.

p_rgh:

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

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

internalField   uniform 0;

boundaryField
{
    front
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    back
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    top
    {
        type            totalPressure;
        p0              uniform 0;
        U               U;
        phi             phi;
        rho             rho;
        psi             none;
        gamma           1;
        value           uniform 0;
    }
    bottom
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    outlet
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0;
    }
    labyrinth_labyrinth
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    inlet
    {
        type            buoyantPressure;
        value           uniform 0;
    }
    aboveInlet
    {
        type            buoyantPressure;
        value           uniform 0;
    }
}

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

U:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform (0 0 0);

boundaryField
{
    top
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }
    bottom
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    inlet
    {
        type            fixedValue;
        value           uniform (2.407 0 0);
    }
    aboveInlet
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    outlet
    {
        type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           uniform (0 0 0);
    }
    back
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    front
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    labyrinth_labyrinth
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}


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

[mgdenno]

Also, I don't know if you need this information, but the floor (bottom) of the domain is at elevation 1 m, the spillway crest is at elevation 4.66 m (so it is 3.66 m tall), the water level at the inlet is at about 6.49 m (5.49 m deep) and the water level at the outlet is not known.

[Eren10]

maybe you should not use inletoutlet things for your case, you have inlet and outlet patch, define it separately

[kumar]

Hello ,
The elevation may be required to define the initial conditions, like using setFields to define the volume fraction =1 on certain cells.

And regarding the boundary conditions for inlet and outlet, I would say since you know the velocity at your inlet you should use total pressure boundary condition at the outlet as well.

Have a look at the tutorial of les in the interFoam tutotials directory. It is a nozzle with velocity inlet and total pressure at other patches.

Try to understand this tutorial case and define similar inlet and outlet conditions.

goodluck
regards
K.Suresh kumar

[mgdenno]

Hi Kumar,

Thanks for the tips. I will try changing the outlet boundary condition and report back how it works out.

MD

[mgdenno]

Hello,

As suggested, I changed the outlet pressure boundary condition to totalPressure to match the atmosphere (same parameters as the atmosphere patch. Is that what you had in mind kumar?, Pablo?), but unfortunately it is still blowing up. It seems to be running more smoothly for the first 16 seconds or so, but ultimately deltaT gets very small and it fails.

Is there a way to tell where in the domain the problem is occurring?

Could this behavior be caused by a mesh that is too course? My mesh is a little course right now as I was hoping to run it this way to a steady state ( to save computation time) and then refine the mesh, map the fields and run it some more. Maybe this is a flawed approach?

As before if there are any particular files (besides those above) that would assisting in troubleshooting, I would be happy to provide them.

Any other suggestions?

Thanks for all you help so far.

MD

[pablodecastillo]

totalPressure for the outlet was the idea.

If you have a good mesh quality, interFoam must work okay, so check mesh. And let us to known your schemes for laplacian etc..., maxCo .........

If you want the steadystate, you can force to PISO with only one step and using relaxation factors.

Pablo

[mgdenno]

Hi Pablo,

Thanks for your help.

I believe that the mesh is OK. I have been using checkMesh and it says that it is OK...I am not sure what the most important result is other than it says it is OK, but I included the checkMesh log if you want to take a look at it.

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.0.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.0.0-a317a4e7cd55
Exec   : checkMesh
Date   : Jul 10 2011
Time   : 09:16:38
Host   : vm-xtide64
PID    : 25141
Case   : /media/Storage/OpenFOAM/mdenno-2.0.0/rasLabyrinthCourse
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           492088
    faces:            1342113
    internal faces:   1259517
    cells:            425572
    boundary patches: 8
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     400296
    prisms:        6465
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     18811

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Surface topology                  
    top                 15840    16193    ok (non-closed singly connected)  
    bottom              9910     11318    ok (non-closed singly connected)  
    inlet               432      481      ok (non-closed singly connected)  
    aboveInlet          144      185      ok (non-closed singly connected)  
    outlet              576      629      ok (non-closed singly connected)  
    back                6345     6710     ok (non-closed singly connected)  
    front               6399     6765     ok (non-closed singly connected)  
    labyrinth_labyrinth 42950    44045    ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-50 1 1) (50 9 19.28)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (2.23571e-16 -2.4443e-14 2.04545e-16) OK.
    Max cell openness = 4.36123e-16 OK.
    Max aspect ratio = 7.2213 OK.
    Minumum face area = 0.00252364. Maximum face area = 0.279707.  Face area magnitudes OK.
    Min volume = 0.000315588. Max volume = 0.139653.  Total volume = 14462.3.  Cell volumes OK.
    Mesh non-orthogonality Max: 61.9855 average: 7.03187
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.18436 OK.

Mesh OK.

End

Here are my fvSchemes and my fvSolution files. I am still figuring it all out...at this point they are mostly the same as the damBreak tutorial. Since the damBreak tutorial is not steady state it seems there maybe some efficiency to be gained here...not sure how though.

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    div(rho*phi,U)  Gauss linear;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,R)      Gauss upwind;
    div(R)          Gauss linear;
    div(phi,nuTilda) Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha;
}


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

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

solvers
{
    pcorr
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-10;
        relTol          0;
    }

    p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-07;
        relTol          0.05;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-07;
        relTol          0;
    }

    "(U|k|epsilon)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }

    "(U|k|epsilon)Final"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor no;
    nCorrectors     3;
    nNonOrthogonalCorrectors 1;
    nAlphaCorr      1;
    nAlphaSubCycles 4;
    
    // level of interface compression. 0 = none, 1 = conservative, >1 = enhanced
    // cAlpha = 1 is recommended
    cAlpha          1;
}


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

I am really interested in the steady state for this simulation. How would I go about forcing PISO with only one step and what settings would you recommend for relaxation factors (which variables and what values)?

Would changing the PISO and relaxation factors help it run better and not crash or just faster?

Thanks again for your help.

MD

[pablodecastillo]

You can try gradSchemes
default cellMDLimited Gauss linear 1; and div(rho*phi,U) Gauss linearUpwind cellLimited Gauss linear 1;, i hope that it is going to improve.

[kumar]

Hi,
If I remember correctly, I used to have similar problems because of three reasons:
1) The boundary conditions.

2) or the mesh .

3) or if the velocity increases suddenly very high (which has something to do with mesh as well)
I think yu should try to make a simple mesh with almost uniform cells and see if you are having simillar problem. Atleast that will let you know if your boundary conditions are completely correct or not.

I am telling this because sometimes the definition of initial time step (max time step) in the controldict file should also be appropriate.

regard
K.Suresh kumar

[kumar]

Just one question I thought thatyou are using the same settings as damBreak case. But in your fvSolution file you are also using some k epsilon settings.

I am a bit confused are you using some turbulence model .

regards
K.Suresh kumar

[mgdenno]

Hello Pablo and kumar,

Thanks so much for your help so far. I really apprecieate it and I think we/I am making progress.

Pablo - I changed my fvSchemes as you suggested (I think). It is definitly running more smoothly/more stable now. However, I think that I may not have the fvSchemes defined incorrectly as I am getting a warning when I run it. Also, I am not so sure about the results...it is showing no water coming over the weir where the weir meets the wall (see attached picture). I don't recall seeing that in the pysical model pictures, but I will have to see if there are any pictures of that location available.

My fvSchemes are now:

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    //default         Gauss linear;
    default         cellMDLimited Gauss linear 1;
}

divSchemes
{
    //div(rho*phi,U)  Gauss linear;
    div(rho*phi,U) Gauss linearUpwind cellLimited Gauss linear 1;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,R)      Gauss upwind;
    div(R)          Gauss linear;
    div(phi,nuTilda) Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha;
}


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

The warning is:

From function linearUpwind(const fvMesh&, Istream&)in file interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.H at line 146 Reading "/media/Storage/OpenFOAM/mdenno-2.0.0/rasLabyrinthCourse/processor0/../system/fvSchemes::divSchemes::div(rho*phi,U)" at line 32 unexpected additional entries in stream.
Only the name of the gradient scheme in the 'gradSchemes' dictionary should be specified.

kumar - I am indeed using the k-epsilon turbulance model. my fvSolution and fvSchemes files are from the ../tutorials/multiphase/interFoam/ras/damBreak case. I should have been more specific in my first post as there are multiple damBreak tutorial cases. I previously suspected that maybe the turbulance was my problem so I tried to also run it as laminar but it still blew up. Does turbulance change your thoughts?

Regarding the simple mesh, are you suggesting a simpler geometry (i.e. a straight wall) or just trying to simplify the mesh for the current geometry? I ran a simple 2D unit width spillway case before trying this 3D case with these same boundary conditions and it worked quite well.

Thanks again,

MD

[mgdenno]

Hello,

Just to provide a little extra information, while I have some hesitation regarding the small "dry" areas at the walls, at time = 50 seconds the flow vs depth upstream of the weir is very close to the other published CFD results and the physical model. I think that it needs further refinement but it is getting there.

MD

[JonW]

Originally Posted by mgdenno

Hello,

Just to provide a little extra information, while I have some hesitation regarding the small "dry" areas at the walls, at time = 50 seconds the flow vs depth upstream of the weir is very close to the other published CFD results and the physical model. I think that it needs further refinement but it is getting there.

MD

Hi Matthew

I would like to suggest one test.
Try cAlpha = 0 with our original case (the cfd result will be useless, but that is not the point).

If it runs, then I think I know what the problem is (at least where it is pointing to).

Cheers
Jon