[David*]

Hey Foamers,

currently I am using groovyBC to simulate Airy waves in a channel with interFoam. For a small wave height (H=0.1m, T=2s --> L=4,74m) everything works fine. But when I try to change to bigger waves (H=0.22m, T=3s --> L=7,68m) my case crashes after 3.48 s.

The stunning part is, that the residuals look quite good until this point and it crashes all out of a sudden.

I guess that my schemes and solution settings are responsible for this behavior, and ask you kindly, to provide me with some tips for improving those settings.

The residual plots are:




And the settings are:
controldict

adjustTimeStep  yes;

maxCo           0.6;
maxAlphaCo    0.6;

maxDeltaT       0.01;

fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         Gauss linear;
}


//V=Vectorfield
divSchemes
{
    div(rho*phi,U)  Gauss limitedLinearV 1.0;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;
    div(phi,k)      Gauss linear;
    div(phi,omega)  Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}
/* corrected: unbounded, 2. order, conservative
   uncorrected: bounded, 1. order, non-conservative
   limited 0.5: blend corrected, uncorrected
*/

interpolationSchemes
{
    default         linear;
    //interpolate(U)  linear;
}

snGradSchemes
{
    default         corrected;
}

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

fvSolution

solvers
{
    pcorr
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-05;
            relTol          0;
            smoother        DICGaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nBottomSweeps   2;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 20;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }
        tolerance       1e-06;
        relTol          0;
        maxIter         100;
    }

    p_rgh
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-05;
            relTol          0;
            smoother        DICGaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nBottomSweeps   2;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 20;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }
        tolerance       1e-06;
        relTol          0;
        maxIter         100;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-06;
            relTol          0;
            nVcycles        2;
            smoother        DICGaussSeidel;
            nPreSweeps      2;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 20;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       1e-06;
        relTol          0;
        maxIter         100;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-06;
        relTol          0;
        nSweeps         1;
    }

    k
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }

    omega
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }
}

PISO
{
    momentumPredictor no;
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;        // non-orthogonality Max < 70: 0, otherwise 1-2, max. 3
    nAlphaCorr      1;
    nAlphaSubCycles 3;                // # sub-cycles, 2: 2x half lengh time step within each actual time step
    cAlpha          1.0;            // 0: no compression, 1: conservative compression, >1: enhanced compression
}

Thank you very much!

David

[David*]

For the record: The crash was due to a poorly set inlet boundary condition.

But nevertheless, it would be interesting to know what's the difference between pcorr, p_rgh and p_rghFinal. Does somebody know that?

Thanks in advance,
David

[tfuwa]

Hi David,

I am also using groovyBC to generate waves. The simulation crashed at a certain point when high wave height or 3D wave tank were tried. The simulation for low wave height and 2D can run smoothly. I suspect I am facing the same problem as you did. Can you please post here the solution to your problem? Many thanks.

Kind regards,
Albert

[akidess]

Originally Posted by David*

But nevertheless, it would be interesting to know what's the difference between pcorr, p_rgh and p_rghFinal. Does somebody know that?

pCorr is a single correction done before the start of the time loop to improve the initial value of the pressure field. p_rgh refers to the pressure solution excluding the final correction, and p_rghFinal to the previously excluded final correction. It is beneficial for computational performance to solve p_rgh with a lower tolerance than p_rghFinal.

[David*]

Originally Posted by tfuwa

Hi David,

I am also using groovyBC to generate waves. The simulation crashed at a certain point when high wave height or 3D wave tank were tried. The simulation for low wave height and 2D can run smoothly. I suspect I am facing the same problem as you did. Can you please post here the solution to your problem? Many thanks.

Kind regards,
Albert

Hi Albert,
as far as I remember, it was a typo in the inlet BC. Re-check it with the sources, e.g. USACE CEM. Another tip is to translate the function by pi/2. By doing this, you start with a neutral water level and not with a crest or through at the inlet.
Regards,
David

edit: thanks @ anton

[eugene]

It doesn't look like it was the problem in this case, but in general a lot of the crashes in interFoam are caused by collapsing air pockets which, due to the large density differences, can drive the velocity to increase very rapidly. To address this you can instead use the "compressible" interFoam solver or put an explicit limiter on the air side velocity magnitude.

[tfuwa]

Hi Eugene,

I saw a treatment suggested by Eric Paterson to ignore the convective term on air-side by changing Ueqn to

fvVectorMatrix UEqn
 (
 fvm::ddt(rho, U)
 + gamma*fvm::div(rhoPhi, U)  // change occurs here 
 - fvm::laplacian(muf, U, "laplacian(mut,U)")
 - (fvc::grad(U) & fvc::grad(muf))
 );

But, I failed to get any improvements. The velocity of air at free-surface is still very large pushing the time step to be small and always causing crashes.

Can you please explain how to "use the compressible interFoam" and especially how to "put an explicit limiter on the air side" to address this problem? Or please give some references.

Kind regards,
Albert

Edit: Thanks, David.