[DoQuocVu]

Dear Foamers,


I implemented a Navier-Stokes solver using a 3 step Runge-Kutta scheme following these steps:
1.png
The solver works fine with orthogonal mesh as I validated it with the lid-driven cavity. However, when it comes to the flow over cylinder problem, where the mesh consists of multi-block and is non-orthogonal, the simulation diverged. The results before it crashed:

2.png
As far as I know with simpleFoam or pisoFoam, nonOrthogonality can be treated by choosing a corrected scheme for the laplacianSchemes and the snGradSchemes, along with nonOrthogonalCorrectors loops.



In my current solver, it clearly doesn't need non-Orthogonal correction loops. I think that when choosing the corrected scheme, the correction term would be automatically evaluated in the discretization process of the laplacian term. So I expected the following setup would work, but it doesn't:

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearlimitedLinearV 1;
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,epsilon) Gauss limitedLinear 1;
    div(phi,R)      Gauss limitedLinear 1;
    div(R)          Gauss linear;
    div(U)            Gauss linear;
    div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear; 
    
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

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

snGradSchemes
{
    default         corrected;
}

Am I wrong about the corrected scheme? I would appreciate any suggestion!

[DoQuocVu]

Here is the main loop of the solver:

while (runTime.run())
{
    #include "readTimeControls.H"
    #include "CourantNo.H"
    #include "setDeltaT.H"
    runTime++;
    Info<< "Time = " << runTime.timeName() << nl << endl;
    const dimensionedScalar dT = runTime.deltaT();

    for(int k=1; k<nStep+1; k++)
    {
        //.1- Explicitly calculate u tildle.
        if (k==1) 
        {
            Uk_1 = U.oldTime();
            pk_1 = p.oldTime();
            surfaceScalarField phik_1
            (
                "phik_1",
                linearInterpolate(Uk_1) & mesh.Sf()
            );
            U_tildle = Uk_1 + dT*(  2.0*alpha[k]*nu*fvc::laplacian(Uk_1)
                                  - 2.0*alpha[k]*fvc::grad(pk_1)
                                  - gamma[k]*fvc::div(phik_1, Uk_1)
                                 );
        }
        else
        {
            surfaceScalarField phik_1
            (
                "phik_1",
                linearInterpolate(Uk_1) & mesh.Sf()
            );
            surfaceScalarField phik_2
            (
                "phik_2",
                linearInterpolate(Uk_2) & mesh.Sf()
            );
            U_tildle = Uk_1 + dT*(   2.0*alpha[k]*nu*fvc::laplacian(Uk_1)
                                   - 2.0*alpha[k]*fvc::grad(pk_1)
                                   - gamma[k]*fvc::div(phik_1, Uk_1)
                                   - zeta[k]*fvc::div(phik_2, Uk_2)
                                 );
        }
        //2. calculating UStar
        volVectorField UStar(U);
        forAll(UStar.internalField(), cellI)
        {
            UStar[cellI] = vector::zero;
        }
        fvVectorMatrix UStarEqn
        (
            fvm::laplacian(UStar)
        );
        for (label i=0; i<mesh.nCells(); i++)
        {
            UStarEqn.diag()[i] -= mesh.V()[i]/(alpha[k]*nu.value()*dT.value());
        }

        solve
        (
            UStarEqn 
            == 
                - (U_tildle/dT + ibForce_ + gradP)/(nu*alpha[k]) 
                + fvc::laplacian(Uk_1) 
        );

        //3. Solve for immediate pressure
        volScalarField pImd(p);
        forAll(pImd.internalField(), cellI)
        pImd[cellI] = 0.0;
                
        fvScalarMatrix pImdEqn
        (
            fvm::laplacian(pImd) == fvc::div(UStar)/(2.0*alpha[k]*dT)
        );

        pImdEqn.solve();

       //4. Correct velocity
        Uk = UStar - 2.0*alpha[k]*dT*fvc::grad(pImd);
            
       //5. Update pressure 
        pk = pk_1 + pImd - alpha[k]*nu*dT*fvc::laplacian(pImd);

     Uk_2 = Uk_1;
     Uk_1 = Uk;
     pk_1 = pk;
     }

U = Uk;
p = pk;
#include "continuityErrs.H"
phi = fvc::flux(U);
U.correctBoundaryConditions();
        
runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
    << "  ClockTime = " << runTime.elapsedClockTime() << " s"
    << nl << endl;
}

[Santiago]

Which is the source/article on which you base the algorithm?

[DoQuocVu]

Hi Santiago,


Its been a while and until now I have time to back to this work.
Actually, I was implementing the immersed boundary method as in this journal:
https://arxiv.org/abs/1809.08170
The immersed boundary technique, however, does not work well with the PISO algorithm as I expected. That's why I wanted to try it with the 3step runge Kutta described in the article.