[djbungee]

Dear all,

as already discussed in the forum "Questions & Answers related to OpenFOAM®", in the article thread "pimpleFoam, transientSimpleFoam, simpleFoam, turbFoam, oodles, transientSimpleOodles", in the XING OpenFOAM-group, the CFL-number in most of the transient solvers provided with the OpenFOAM library is not required, and to limit SIMPLE to steady simulations is also wrong. Moreover, in these solvers, you do actually not get convergence on an (arbitrarily) large time step and, therefore, experience convergence problems for larger time steps.

What I consider important in order to achieve an implicit procedure with arbitrary time-step sizes is that all discretization is done in an implicit manner (fvm::....) and that the sequentially solved equations are solved for _several_ times in a single time step to reach a converged _SYSTEM_ of equations on a time step. Below and attached you will find an early development version of our code with fixed convergence criteria and some features that still need some improvements. However, the basic idea should become apparent. With this code you can achieve large time steps, and if the initial residuals in the inner loop below the time-step loop do not decrease, the time-step is probably physically too 'coarse' or too large and should be decreased.

Suggestions and discussions are welcome. We want to discuss this on the Stammtisch in Berlin tomorrow, for anyone interested in the bof-session 'Hybrid RANS/LES and DES'.

Best regards, Ulf.

 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
    transientSimpleOodles

Description
    Incompressible LES solver using transient SIMPLE.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "incompressible/transportModel/transportModel.H"
#include "LESModel.H"
#include "IFstream.H"
#include "OFstream.H"
#include "Random.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"//#include "createMeshNoClear.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"
    #include "readTransportProperties.H"

    Info<< "\nStarting time loop\n" << endl;

    for (runTime++; !runTime.end(); runTime++)
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "readPISOControls.H"
        #include "CourantNo.H"
    scalar eqnResidual = 1;
    scalar convergenceCriterion = 1.E-20;
    scalar maxResidual = convergenceCriterion + 1;
    scalar MaxResidualChange = 0;
    scalar MaxResidualPrev = 0;
    int IIcount=0;
        // Pressure-velocity SIMPLE corrector
        //for (int corr=0; corr<nCorr; corr++)
    while(maxResidual > convergenceCriterion)
        {
         // Momentum predictor
         eqnResidual = 1; 
         MaxResidualPrev=maxResidual;
         maxResidual = 0;
         IIcount++;

         tmp<fvVectorMatrix> UEqn
           (
                  fvm::ddt(U)
        + fvm::div(phi, U)
        + sgsModel->divDevBeff(U)
        );

         UEqn().relax();

        eqnResidual = solve
          (
           UEqn() == -fvc::grad(p)
           ).initialResidual();
        
        //  Retain maximum residual for global convergence check
        maxResidual = max(eqnResidual, maxResidual);

            p.boundaryField().updateCoeffs();
        
        volScalarField rUA = 1.0/UEqn().A();
        U = rUA*UEqn().H();
        UEqn.clear();
        
        phi = (fvc::interpolate(U) & mesh.Sf())
          + fvc::ddtPhiCorr(rUA, U, phi);
        
        adjustPhi(phi, U, p);
        
        // Non-orthogonal pressure corrector loop
        if (nCorr != 1) //Ulf aus pimple
          {
        p.storePrevIter();
          }
        //////--------------------------------------------------------------------------------------
        scalar StoppPressCorr =1;
        int CountPressCorr =0;
        scalar PrevPressEqnResidual=1.E20;
        //for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
        while ( (StoppPressCorr > convergenceCriterion) && (CountPressCorr < 10) )
          {
        CountPressCorr++;
        fvScalarMatrix pEqn
          (
           fvm::laplacian(rUA, p) == fvc::div(phi)
           );
        
        pEqn.setReference(pRefCell, pRefValue);

        eqnResidual = pEqn.solve().initialResidual();
        StoppPressCorr = eqnResidual;

        // Stop loop if residual increases, i.e. set loop counter to 10
        if (eqnResidual > PrevPressEqnResidual) CountPressCorr=10;
        PrevPressEqnResidual = eqnResidual;

        if (StoppPressCorr <= convergenceCriterion)
          {
            phi -= pEqn.flux();
            maxResidual = max(eqnResidual, maxResidual);
          }
        if (CountPressCorr == 10)
          {
            // Do not save maxResidual so that the outer, U-based maxResidual is not influenced if convergence is not reached 
            phi -= pEqn.flux();
          }
          }
        ////----------------------------------------------------------------------------------------
        
#           include "continuityErrs.H"

        // Explicitly relax pressure for momentum corrector
        p.relax();

        // Ulf, nach unten geschoben: DAS SOLLTE HINTER MOMENTUM CORRECTOR:  sgsModel->correct();
        
        // Momentum corrector
        U -= rUA*fvc::grad(p);
        U.correctBoundaryConditions();
        
        sgsModel->correct();
        // Info << "   ->  maxResidual(U,p) " << maxResidual << endl;


        // CONVERGENCE CHECK 

        MaxResidualChange=fabs(MaxResidualPrev - maxResidual);
        if (IIcount > 100 ) MaxResidualChange = MaxResidualChange * 100./IIcount;
        // If residual does not change anymore, stop loop 
        if ( (MaxResidualChange/maxResidual < 0.001) && (IIcount!=1) )
          {
        Info << "   ->  maxResidual(U,p) " << maxResidual << endl;
        Info << "   ->  MaxResidualChange (U,p)   " <<MaxResidualChange <<endl;
        Info << "   ->  MaxResidualChange/maxResidual   " <<MaxResidualChange/maxResidual <<endl;
        Info << "   ->  NOT reached convergence criterion " << convergenceCriterion << " on time step " << runTime.timeName() << " after " << IIcount <<" cycles " <<endl;
                Info << "   ->  maxResidual (U,p) set to 0 to end loop due to low MaxResidualChange" <<endl;
        maxResidual=0;
          }
        
        // FINALLY REACHED CONVERGENCE
        if (maxResidual < convergenceCriterion)
          {
        Info << "   ->  reached convergence criterion " << convergenceCriterion << " on time step " << runTime.timeName() << " after " << IIcount <<" cycles " <<endl;
          }

        }
    
        runTime.write();

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

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

    return(0);
}


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