[dalaron]

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:

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:

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:

#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.

[dalaron]

A semi-progress update. Increasing the Courant number from 0.2 to 0.5 solves this problem and gives me good results.

As far as I understand, the Courant number only affects the stability of the solver. So shouldn't a higher Courant number be more likely to cause problems? What have I done wrong here that a lower Courant number is doing this instead?