[ch_jaehnel]

Hallo Foamers,


I am simulating an vortex power plant with interDyMFoam with AMI and an steady and a rotating mesh. As turbulence model I use SpalartAlmares IDDES After several seconds of smooth simulating the Courant number blows up and the simulation fails. When I look in the latest results I have a few cells (at the water surface, close to the turbine) with velocity magnitudes higher than 1000m/s, which are obviously the problem. The cells and the mesh look quite okay, as I would say. This problem occurs in different meshes at different timesteps. Mostly I can start the simulation at an earlier timestep and it runs through the problematic timestep. Here are some pictures and my boundary conditions.



U-file:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5.0                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      binary;
    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);
    }
    inlet
    {
        type            variableHeightFlowRateInletVelocity;
        flowRate        0.701;
        alpha           alpha.water;
        value           uniform (0 0 0);
    }
    outlet.water
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    outlet.air
    {
        type            zeroGradient;
    }
    SohleOW
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    SohleUW
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    SohleWS
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    WaendeOW
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    WaendeUW
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    WSBecken
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    turbine
    {
        type            movingWallVelocity;
        value           uniform (0 0 0);
    }
    AMI1
    {
        type            cyclicAMI;
        value           uniform (0 0 0);
    }
    AMI2
    {
        type            cyclicAMI;
        value           uniform (0 0 0);
    }
}


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

p_rgh:

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

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

internalField   uniform 0;

boundaryField
{
    top
    {
        type            totalPressure;
        rho             rho;
        psi             none;
        gamma           1;
        p0              uniform 0;
        value           uniform 0;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet.water
    {
        type            zeroGradient;
    }
    outlet.air
    {
        type            totalPressure;
        rho             rho;
        psi             none;
        gamma           1;
        p0              uniform 0;
        value           uniform 0;
    }
    SohleOW
    {
        type            zeroGradient;
    }
    SohleUW
    {
        type            zeroGradient;
    }
    SohleWS
    {
        type            zeroGradient;
    }
    WaendeOW
    {
        type            zeroGradient;
    }
    WaendeUW
    {
        type            zeroGradient;
    }
    WSBecken
    {
        type            zeroGradient;
    }
    turbine
    {
        type            zeroGradient;
    }
    AMI1
    {
        type            cyclicAMI;
        value           uniform 0;
    }
    AMI2
    {
        type            cyclicAMI;
        value           uniform 0;
    }
}


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

dynamicMeshDict:

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

dynamicFvMesh   solidBodyMotionFvMesh;

motionSolverLibs ( "libfvMotionSolvers.so" );

solidBodyMotionFvMeshCoeffs
{
    cellZone        rotating;

    solidBodyMotionFunction  rotatingMotion;
    rotatingMotionCoeffs
    {

        origin        (7.98 2.15 0.9);
        axis          (0 0 1);
        omega         -2.5447; // 5 rad/s
    }
}
// ************************************************************************* //

I would be very grateful if some of you could help me. I'm working on this since weeks and it drives me crazy. Thank you in advance!

Christian

[cojua8]

I don't know what may be the cause of this, but you can try adding

adjustTimeStep true;
maxCo <number less than 1>;

in your case's controlDict

[indy07cz]

Hi, I had similar problem with interFoam. You should use bounding for gradU in fvSchemes, there is topic about it. Definitely use adjustTimestep and Co less then 1. For VOF computations Co less than 0.1 is recommended but it depends on mesh size and cells volume ratio.

[RobertHB]

Furthermore, try to run a checkMesh. If you get to many errors or skewed cells, it's likely that your mesh is the cause for your problems.

[ch_jaehnel]

Thank you for your replies.


here is my checkMesh:

Create time

Create polyMesh for time = constant

Time = constant

Mesh stats
    points:           12571351
    faces:            32443085
    internal faces:   29757718
    cells:            10044785
    faces per cell:   6.19235
    boundary patches: 13
    point zones:      0
    face zones:       1
    cell zones:       1

Overall number of cells of each type:
    hexahedra:     8616751
    prisms:        235466
    wedges:        0
    pyramids:      0
    tet wedges:    671
    tetrahedra:    24
    polyhedra:     1191873
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            4   210236
            5   144169
            6   207571
            7   224
            8   230
            9   422636
           10   19
           11   4
           12   132870
           13   2
           15   73092
           18   820

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
   *Number of regions: 2
    The mesh has multiple regions which are not connected by any face.
  <<Writing region information to "constant/cellToRegion"
  <<Writing region 0 with 3814383 cells to cellSet region0
  <<Writing region 1 with 6230402 cells to cellSet region1

Checking patch topology for multiply connected surfaces...
                   Patch    Faces   Points                  Surface topology
                     top    33807    35382  ok (non-closed singly connected)
                   inlet     1212     1361  ok (non-closed singly connected)
            outlet.water      482      586  ok (non-closed singly connected)
              outlet.air      558      639  ok (non-closed singly connected)
                 SohleOW    16892    17262  ok (non-closed singly connected)
                 SohleUW    13280    13603  ok (non-closed singly connected)
                 SohleWS   459719   466049  ok (non-closed singly connected)
                WaendeOW    26801    27627  ok (non-closed singly connected)
                WaendeUW    16280    16974  ok (non-closed singly connected)
                WSBecken   133829   135355  ok (non-closed singly connected)
                 turbine  1267887  1480133      ok (closed singly connected)
                    AMI1   357310   358590      ok (closed singly connected)
                    AMI2   357310   358590      ok (closed singly connected)

Checking geometry...
    Overall domain bounding box (0 0.27509 0.100001) (15 4.41 2.3)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-1.70651e-16 1.27318e-16 -2.97151e-13) OK.
    Max cell openness = 4.24293e-16 OK.
    Max aspect ratio = 12.0035 OK.
    Minimum face area = 1.25745e-08. Maximum face area = 0.00741749.  Face area magnitudes OK.
    Min volume = 6.17915e-10. Max volume = 0.000250245.  Total volume = 41.298.  Cell volumes OK.
    Mesh non-orthogonality Max: 78.2713 average: 10.0852
   *Number of severely non-orthogonal (> 70 degrees) faces: 4.
    Non-orthogonality check OK.
  <<Writing 4 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
 ***Max skewness = 5.36521, 117 highly skew faces detected which may impair the quality of the results
  <<Writing 117 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

Failed 1 mesh checks.

Indeed there are a few skewed faces. Due to a bit difficult geometry, it was the best result I could get. I am just wondering why the problems occurs after 9seconds of simulation. So the turbine already made 4-5 rotations.



The other tip with fv-Schemes sounds interesting. Maybe you can find a problem in my fv-Schemes:

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default             none;

    div(rhoPhi,U)      Gauss linearUpwind grad(U);
    div(phi,alpha)      Gauss vanLeer;
    div(phirb,alpha)    Gauss linear;
    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((muEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

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


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

I will try to find out more about "bounding for gradU in fvSchemes"



The simulation already runs with adjustableTimeStep and a maxCo of 0.95. Though I recongnized that mostly is possible to compute over these problematic time steps without an adjustableTimeStep and quite small timesteps of maybe 5e-5 or less. With adjustableTimeStep thats mostly not possible as it becomes less than 1e-100 with velocity magnitudes of 1e5 or even more.

[ch_jaehnel]

Originally Posted by indy07cz

Hi, I had similar problem with interFoam. You should use bounding for gradU in fvSchemes, there is topic about it. Definitely use adjustTimestep and Co less then 1. For VOF computations Co less than 0.1 is recommended but it depends on mesh size and cells volume ratio.

Why would you say Co has to be less than 0.1 in VOF? Never heard of this criteria... Is'nt it just wasted computation time?

[indy07cz]

Hi, here is part of fvSchemes you can use for your case

gradSchemes
{
    default         Gauss linear;
    grad(U)        cellLimited Gauss linear 1; 

}

About the Co criteria, I've made several simulations with interFoam with large geometries (spillway, dam chute and so on) and mostly I had to use Co about 0.1-0.05 because of high velocities and high cell volume ratio. If I used high Co I got alpha.water from 0 to 2 which is wrong. But as I said, it depends on your mesh quality, ratio between smallest and largest cell, velocity etc. Just watch your p_rgh residuals (decreasing Co=decreasing residuals) and simply compute initial Co (or time step) for your mesh from the equation for 3 dimensions https://en.wikipedia.org/wiki/Couran...Lewy_condition

I recommend to you:
1. Use bounding for grad(U) and check mesh again (you can watch skewed cells in Paraview and improve mesh there)
2. Check BC again
3. Use lower Co

[akidess]

Originally Posted by ch_jaehnel

Why would you say Co has to be less than 0.1 in VOF? Never heard of this criteria... Is'nt it just wasted computation time?

It is well known that the stability criterion for the standard VoF algorithm is much more severe for fine meshes than the Courant criterion for pressure-velocity coupling. Check the original paper of Brackbill, and for even more details the paper by Deshpande (everyone using interfoam should know these two papers).