Development of local time stepping solver for incompressible flow

I'm very much interested in local time stepping method but there's not an LTS version of pimpleFoam (LTSPimpleFoam) in the standard solvers of OpenFOAM.

And so I decided to make my own LTSPimpleFoam.

1. First I copied $FOAM_SOLVERS/incompressible/pimpleFoam directory to some appropriate directory (say $WM_PROJECT_DIR/applications/solvers) and renamed the directory as myLTSPimpleFoam,

2. and deleted pimpleDyMFoam, SRFPimpleFoam, and Make/linux64GccDPOpt and pimpleFoam.dep from myLTSPimpleFoam directory and renamed pimpleFoam.C as myLTSPimpleFoam.C.

3. Next I edited Make/files as follows:

    myLTSPimpleFoam.C
    
    EXE = $(FOAM_USER_APPBIN)/myLTSPimpleFoam

4. From here I started my customizing of the solver. I copied setInitialrDeltaT.H and setrDeltaT.H from $FOAM_SOLVERS/compressible/rhoPimpleFoam/rhoLTSPimpleFoam to myLTSPimpleFoam directory,

5.1. and edited myLTSPimpleFoam.C (differences from $FOAM_SOLVERS/incompressible/pimpleFoam/pimpleFoam.C are emphasized in red),

    #include "fvCFD.H"
    #include "singlePhaseTransportModel.H"
    #include "turbulenceModel.H"
    #include "pimpleControl.H"
    #include "fvIOoptionList.H"
    #include "fvcSmooth.H"
    #include "IOporosityModelList.H"
    #include "IOMRFZoneList.H"
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    int main(int argc, char *argv[])
    {
        #include "setRootCase.H"
        #include "createTime.H"
        #include "createMesh.H"
        #include "createFields.H"
        #include "createFvOptions.H"
        #include "initContinuityErrs.H"
    
        pimpleControl pimple(mesh);
    
        #include "setInitialrDeltaT.H"
    
        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
        Info<< "\nStarting time loop\n" << endl;

        while (runTime.run())
        {
            #include "readTimeControls.H"
            #include "CourantNo.H"
            #include "setDeltaT.H"
    
            runTime++;
    
            Info<< "Time = " << runTime.timeName() << nl << endl;
    
            #include "setrDeltaT.H"
    
            // --- Pressure-velocity PIMPLE corrector loop
            while (pimple.loop())
            {
                #include "UEqn.H"
    
                // --- Pressure corrector loop
                while (pimple.correct())
                {
                    #include "pEqn.H"
                }
    
                if (pimple.turbCorr())
                {
                    turbulence->correct();
                }
            }
    
            runTime.write();
    
            Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                << nl << endl;
        }
    
        Info<< "End\n" << endl;
    
        return 0;
    }


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

5.2. and setrDeltaT.H (differences from $FOAM_SOLVERS/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/setrDeltaT.H are emphasized in red).

    {
        const dictionary& pimpleDict = pimple.dict();
    
        scalar maxCo
        (
            pimpleDict.lookupOrDefault<scalar>("maxCo", 0.8)
        );

        scalar rDeltaTSmoothingCoeff
        (
            pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
        );
    
        scalar rDeltaTDampingCoeff
        (
            pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
        );
    
        scalar maxDeltaT
        (
            pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
        );
    
        volScalarField rDeltaT0("rDeltaT0", rDeltaT);
    
        // Set the reciprocal time-step from the local Courant number
        rDeltaT.dimensionedInternalField() = max
        (
            1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
            fvc::surfaceSum(mag(phi))().dimensionedInternalField()
           /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
        );
    
        if (pimple.transonic())
        {
            surfaceScalarField phid
            (
                "phid",
                fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
            );
    
            rDeltaT.dimensionedInternalField() = max
            (
                rDeltaT.dimensionedInternalField(),
                fvc::surfaceSum(mag(phid))().dimensionedInternalField()
                /((2*maxCo)*mesh.V()*psi.dimensionedInternalField())
            );
        }
    
        // Update tho boundary values of the reciprocal time-step
        rDeltaT.correctBoundaryConditions();
    
        Info<< "Flow time scale min/max = "
            << gMin(1/rDeltaT.internalField())
            << ", " << gMax(1/rDeltaT.internalField()) << endl;
    
        if (rDeltaTSmoothingCoeff < 1.0)
        {
            fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
        }
    
        Info<< "Smoothed flow time scale min/max = "
            << gMin(1/rDeltaT.internalField())
            << ", " << gMax(1/rDeltaT.internalField()) << endl;
    
        // Limit rate of change of time scale
        // - reduce as much as required
        // - only increase at a fraction of old time scale
        if
        (
            rDeltaTDampingCoeff < 1.0
         && runTime.timeIndex() > runTime.startTimeIndex() + 1
        )
        {
            rDeltaT =
                rDeltaT0
               *max(rDeltaT/rDeltaT0, scalar(1) - rDeltaTDampingCoeff);
    
            Info<< "Damped flow time scale min/max = "
                << gMin(1/rDeltaT.internalField())
                << ", " << gMax(1/rDeltaT.internalField()) << endl;
        }
    }

6. Finally I ran wmake on myLTSPimpleFoam directory and $FOAM_USER_APPBIN/myLTSPimpleFoam was made.

I'm now trying this myLTSPimpleFoam and it seems working fine. I will report the results of steady-state simulations of Ahmed body using myLTSPimpleFoam in the next series of articles.

(Postscript added on 19/10/2014)
Source files of myLTSPimpleFoam can be downloaded here.

Comments