CFD Online Discussion Forums - Modifying sonicFoam to run with variable timestep

I have modified the sonicFoam solver (version 2.1.1) to run a variable timestep and with maximum Courant number. The reason I want to have a variable timestep sonicFoam solver is to perform a grid dependency study, and doing so with a static timestep was very frustrating.

The modified solver works for some grid resolutions and not for others. I don't know why. The test case, and some (incorrect) pressure plots are attached.

At the start of the simulation I set MaxCo = 0.2 in my controlDict file, run own_sonicFoam, and everything appears to be proceeding okay. A typical output near the start of the simulation is:

Code:

Courant Number mean: 0.0629867 max: 0.200099

deltaT = 2.48178e-05

TIME = 0.127543

The max courant number is slightly higher than the specified MaxCo but I think that's okay for now. The timestep updates nicely, increasing or decreasing to keep the max Courant number ~< 0.2.

Then at a certain point the timestep drops rapidly and there is a huge discrepancy between the Courant number mean and max:

Code:

Courant Number mean: 8.67688e-05 max: 0.2

deltaT = 1.53002e-09

TIME = 0.199746794

Looking at the solution for pressure develop over time in paraview, there is a transition from being symmetric about the central line (everything updating nicely, deltaT nice size) to unsymmetrical flow (timestep rapidly drops). What would cause this transition to unsymmetrical? Is there anything I can do to stop it?

One of the more confusing things about this was that for one test case (refer blockMeshDict file) I found that an odd number of cells in the vertical direction caused problems, but an even number gave a good solution.

The changes I made to the sonicFoam solver are:

Code:



#include "fvCFD.H"

#include "basicPsiThermo.H"

#include "turbulenceModel.H"



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



int main(int argc, char *argv[])

{

    #include "setRootCase.H"        

    #include "createTime.H"             

    #include "createMesh.H"            

    #include "createFields.H"

    #include "initContinuityErrs.H"



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



    Info<< "\nSTARTING TIME LOOP\n" << endl;



    while (runTime.loop()) //runTime.run() ?

    {



        #include "readPISOControls.H"

        #include "compressibleCourantNo.H"

        #include "readTimeControls.H" 

        #include "setDeltaT.H" 

        runTime++;    

        

        Info<< "TIME = " << runTime.timeName() << nl << endl;



        //#include "readPISOControls.H"

        //#include "compressibleCourantNo.H"



        #include "rhoEqn.H"



        #include "UEqn.H"



        #include "eEqn.H"

        // --- PISO loop



        for (int corr=0; corr<nCorr; corr++)

        {

            #include "pEqn.H"

        }



        turbulence->correct();



        rho = thermo.rho(); // thermodynamics call.



        runTime.write();



        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

            << "  ClockTime = " << runTime.elapsedClockTime() << " s"

            << nl << endl;

    }



    Info<< "End\n" << endl;



    return 0;

}

My programming experience with C++ and OpenFOAM is limited, so any insights or suggestions would be appreciated.