[Igor_2011]

Hi everyone!

I'm rookie on the forum and I'd like to start by greeting all the members for the very interesting topics that are being discussed here.

I recently tried to simulate the unsteady around a basic 2-D cylinder using the solvers available in OpenFoam, in particular pimpleDyMFoam. Basically the code I use is the looks like:

        // Do any mesh changes
        mesh.update();

        // Make the fluxes relative to the mesh motion
        fvc::makeRelative(phi, U);

        if (mesh.changing() )
        {
            #include "meshCourantNo.H"
        }

       for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
        {
        // Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H"
        
            #include "pEqn.H"

        }
        
        turbulence->correct();
        }

where the velocity and pressure are solved in the following parts:

    // Solve the Momentum equation

    fvVectorMatrix UEqn
    (
    fvm::ddt(U)
       +fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );

    UEqn.relax();

    eqnResidual = solve
    (
        UEqn == -fvc::grad(p)
    ).initialResidual();

    maxResidual = max(eqnResidual, maxResidual);

    p.boundaryField().updateCoeffs();

    volScalarField AU = UEqn.A();
    U = UEqn.H()/AU;
    //UEqn.clear();
    phi = fvc::interpolate(U) & mesh.Sf();
    adjustPhi(phi, U, p);

    // Non-orthogonal pressure corrector loop
    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(1.0/AU, p) == fvc::div(phi)
        );

        pEqn.setReference(pRefCell, pRefValue);

        if (nonOrth == 0)
        {
            eqnResidual = pEqn.solve().initialResidual();
            maxResidual = max(eqnResidual, maxResidual);
        }
        else if( ocorr == nOuterCorr-1 && nonOrth == nNonOrthCorr)
        {
                        pEqn.solve(mesh.solver(p.name() + "Final"));
        }
        else
        {
            pEqn.solve();
        }

        if (nonOrth == nNonOrthCorr)
        {
            phi -= pEqn.flux();
        }
    }

#   include "continuityErrs.H"

    // Explicitly relax pressure for momentum corrector

    p.relax();

    // Momentum corrector
    U -= fvc::grad(p)/AU;
    U.correctBoundaryConditions();

As you can see, two inner loops are performed during each iterations. From top to bottom:

1. loop for non-orthogonality correction on the pressure equation (nCorr
2. velocity- pressure correcion

My question here concerns the initial residual at the beginning of each inner iterations. Typically for nNonOrthCorr=3 and nOuterCorr=2 the following residual history is obtained:

Mesh Courant Number mean: 1.03403 max: 1645.15
smoothSolver:  Solving for Ux, Initial residual = 0.000663609, Final residual = 2.88668e-06, No Iterations 8
smoothSolver:  Solving for Uz, Initial residual = 0.00548648, Final residual = 2.28391e-05, No Iterations 8
DICPCG:  Solving for p, Initial residual = 0.828448, Final residual = 9.83331e-07, No Iterations 506
DICPCG:  Solving for p, Initial residual = 0.211827, Final residual = 9.78121e-07, No Iterations 478
DICPCG:  Solving for p, Initial residual = 0.0673258, Final residual = 9.56774e-07, No Iterations 468
time step continuity errors : sum local = 5.6259e-12, global = 4.41739e-20, cumulative = 1.4129e-19
smoothSolver:  Solving for omega, Initial residual = 7.3622e-06, Final residual = 7.40485e-09, No Iterations 4
smoothSolver:  Solving for k, Initial residual = 0.0114753, Final residual = 3.82431e-06, No Iterations 3

smoothSolver:  Solving for Ux, Initial residual = 0.000132713, Final residual = 1.01532e-06, No Iterations 8
smoothSolver:  Solving for Uz, Initial residual = 0.000974607, Final residual = 6.52005e-06, No Iterations 8
DICPCG:  Solving for p, Initial residual = 0.804658, Final residual = 9.75985e-07, No Iterations 513
DICPCG:  Solving for p, Initial residual = 0.241636, Final residual = 9.50152e-07, No Iterations 481
DICPCG:  Solving for p, Initial residual = 0.0714927, Final residual = 9.76869e-07, No Iterations 458
time step continuity errors : sum local = 5.33943e-12, global = 1.72288e-20, cumulative = 1.58519e-19
smoothSolver:  Solving for omega, Initial residual = 1.39838e-06, Final residual = 8.80479e-09, No Iterations 3
smoothSolver:  Solving for k, Initial residual = 0.00313893, Final residual = 1.69267e-06, No Iterations 3

and I don't understand why my initial residual for the pressure is not equal to the final residual at the previous sub-iteration. In particular, should I consider the case to be converged when the first initial residual is below a pre-defined threshold?

I'm a beginner in CFD so do not hesitate to correct or highlight any misinterpretations or vocabulary misuses.

Great thanks!

Igor

[L1011]

I'm also interested by this post, nobody has a clue?

[Yahoo]

I am trying to solve equation

fvScalarMatrix CEqn
(
fvm::ddt(C)
+ fvm::div(phi, C)

);
CeqnResidual = CEqn.solve().initialResidual();

In case where phi is zero (i.e. no flow), my initial residuals are always high (~ 0.3). When I am forcing the solver to exit without satisfying the convergence criterion (by setting a maximum for the iterations), I see that phi is zero everywhere and C is constant, as I expect. But, I don't know why CeqnResidual is high!