[cwang5]

Hi guys,

I'm trying to use the adjustTimeStep option in controDict to run my LES simulation. However, the courant number reported at run time is so high that OpenFOAM was forced to constantly reduce the time step size until it hits ~1e-30 and the whole simulation just crashed. I revised the simulation to run on fixed time step, and was getting the following message in my log:

Courant Number mean: 0.0042170848 max: 136.39391
Time = 75.642

DILUPBiCG:  Solving for Ux, Initial residual = 0.00032185865, Final residual = 3.2946399e-08, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.00029752553, Final residual = 3.1410567e-08, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.00066511108, Final residual = 7.6629127e-08, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.0046290147, Final residual = 7.3302722e-05, No Iterations 2
GAMG:  Solving for p, Initial residual = 7.3305549e-05, Final residual = 9.6448053e-07, No Iterations 13
time step continuity errors : sum local = 7.9685318e-13, global = 8.7932007e-14, cumulative = 1.3817586e-10
ExecutionTime = 53373.4 s  ClockTime = 53430 s

which reported a maximum courant number of ~135 for every time step. However, the value reported by the Co utility showed that the maximum courant number is ~0.04 when I post processed the data:

Time = 75
    Reading phi
    Calculating incompressible Co
Co max : 0.032720504

End


Time = 76
    Reading phi
    Calculating incompressible Co
Co max : 0.032659777

End

Is it a problem with the adjustTimeStep calculation of the courant number that's causing the issue? And what can I do to be able to run the simulation with adjustTimeStep turned on so I can save on simulation time?

Thanks,
John

EDIT: I should add that the simulation was stable for the entire simulation even with maxCo reported at ~135.
EDIT2: I looked into OpenFOAM source code, and the Courant number in Co is calculated using

(0.5*runTime.deltaT())*fvc::surfaceSum(mag(phi)().dimensionedInternalField()/mash.V();

while it's calculated in CourantNo.C with

(0.5*fvc::surfaceSum(mag(phi))().internalField()/rho(phi)().internalField()//mesh.V().field()*mesh.time().deltaTValue();

The biggest difference is the extra rho(phi)().internalField() term, which accounts for compressibility as stated in CourantNo.H

        //- Return 1 if the flux field is volumetric i.e. incompressible flow?, otherwise return rho
        //  from the database
        tmp<volScalarField> rho(const surfaceScalarField& p) const;

everything other difference is between reading from a output database file vs reading from memory at the time of computing as far as I can tell.
On a side note, I ran the same simulation using a finer mesh (2/3 of original grid size), and obtained significantly lower Courant number of ~0.1 in FOAM log, which enabled me to turn on adjustTimeStep.