Projection Methods

This post is for projection methods in CFD.

PISO Algorithm

// Pressure-velocity PISO corrector
{
    // Solve the Momentum equation

    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
        (
        fvm::ddt(U) + fvm::div(phi, U)
        + MRF.DDt(U)
        + turbulence->divDevSigma(U)
        ==
        fvOptions(U)
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (piso.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
    }

    // --- PISO loop
    while (piso.correct())
    {
        volScalarField rAU(1.0/UEqn.A());
        volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
        surfaceScalarField phiHbyA
            (
            "phiHbyA",
            fvc::flux(HbyA)
            + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi))
        );

        MRF.makeRelative(phiHbyA);

        adjustPhi(phiHbyA, U, p);

        // Update the pressure BCs to ensure flux consistency
        constrainPressure(p, U, phiHbyA, rAU, MRF);

        // Non-orthogonal pressure corrector loop
        while (piso.correctNonOrthogonal())
        {
            // Pressure corrector

            fvScalarMatrix pEqn
            (
                fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
            );

            pEqn.setReference(pRefCell, pRefValue);

            pEqn.solve();

            if (piso.finalNonOrthogonalIter())
            {
                phi = phiHbyA - pEqn.flux();
            }
        }

        #include "continuityErrs.H"

        U = HbyA - rAU*fvc::grad(p);
        U.correctBoundaryConditions();
        fvOptions.correct(U);

    }
}

PIMPLE Algorithm

while (pimple.run(runTime))
{
    #include "readDyMControls.H"

    if (LTS)
    {
        #include "setRDeltaT.H"
    }
    else
    {
        #include "CourantNo.H"
        #include "setDeltaT.H"
    }

    runTime++;

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

    // --- Pressure-velocity PIMPLE corrector loop
    while (pimple.loop())
    {
        if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
        {
            mesh.update();

            if (mesh.changing())
            {
                MRF.update();

                if (correctPhi)
                {
                    // Calculate absolute flux
                    // from the mapped surface velocity
                    phi = mesh.Sf() & Uf();

                    CorrectPhi
                    (
                        U,
                        phi,
                        p,
                        dimensionedScalar("rAUf", dimTime, 1),
                        geometricZeroField(),
                        pimple,
                        true
                    );

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

                if (checkMeshCourantNo)
                {
                    #include "meshCourantNo.H"
                }
            }
        }

        MRF.correctBoundaryVelocity(U);

        tmp<fvVectorMatrix> tUEqn
        (
            fvm::ddt(U) + fvm::div(phi, U)
          + MRF.DDt(U)
          + turbulence->divDevSigma(U)
         ==
            fvOptions(U)
        );
        fvVectorMatrix& UEqn = tUEqn.ref();

        UEqn.relax();

        fvOptions.constrain(UEqn);

        if (pimple.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));

            fvOptions.correct(U);
        }

        // --- Pressure corrector loop
        while (pimple.correct())
        {
            volScalarField rAU(1.0/UEqn.A());
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA)
              + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf))
            );

            MRF.makeRelative(phiHbyA);

            if (p.needReference())
            {
                fvc::makeRelative(phiHbyA, U);
                adjustPhi(phiHbyA, U, p);
                fvc::makeAbsolute(phiHbyA, U);
            }

            tmp<volScalarField> rAtU(rAU);

            if (pimple.consistent())
            {
                rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
                phiHbyA +=
                    fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
                HbyA -= (rAU - rAtU())*fvc::grad(p);
            }

            if (pimple.nCorrPiso() <= 1)
            {
                tUEqn.clear();
            }

            // Update the pressure BCs to ensure flux consistency
            constrainPressure(p, U, phiHbyA, rAtU(), MRF);

            // Non-orthogonal pressure corrector loop
            while (pimple.correctNonOrthogonal())
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);

                pEqn.solve();

                if (pimple.finalNonOrthogonalIter())
                {
                    phi = phiHbyA - pEqn.flux();
                }
            }

            #include "continuityErrs.H"

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

            U = HbyA - rAtU*fvc::grad(p);
            U.correctBoundaryConditions();
            fvOptions.correct(U);

            // Correct Uf if the mesh is moving
            fvc::correctUf(Uf, U, phi);

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

        if (pimple.turbCorr())
        {
            laminarTransport.correct();
            turbulence->correct();
        }
    }
}

Runge kutta Algorithm

U 4th order p 2nd order

for (int cycle =0; cycle < rkCoeff.size(); cycle++)
{
    phi = (fvc::interpolate(U) & mesh.Sf());  // evaluate flux
    dU  = runTime.deltaT()*(fvc::laplacian(turbulence->nuEff(), U) - fvc::div(phi,U));
    Uc  = Uc + rkCoeff[cycle]*dU;
    U   = Uold + rkCoeff1[cycle]*dU;
    if(cycle==3) U=Uc;
    U.correctBoundaryConditions();
    solve( pEqn == fvc::div(U)/runTime.deltaT());
    U = U -fvc::grad(p)*runTime.deltaT();
    U.correctBoundaryConditions();
    if(cycle==1) p2 = p;
    if(cycle==3) p = 4*p - 3*p2;
}

phi = (fvc::interpolate(U) & mesh.Sf());
turbulence->correct();
#include "continuityErrs.H"

Rhie-Chow Algorithm

U p 2nd order

https://github.com/asimonder/projectionFoam

while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
	
        #include "readPISOControls.H"
        #include "CourantNo.H"

        // Non-iterative projection scheme
        {
	    //linearization of convective flux using second order extrapolation  	  
	    surfaceScalarField phi_o(0.5*(3.0*phi.oldTime()-phi.oldTime().oldTime()));
	  
	    // Momentum step

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

	    solve(UEqn == -fvc::grad(p));
 
            // Projection step

	    dimensionedScalar dt=runTime.deltaT();
	    scalar rDeltaT=1.0/dt.value();

	    U += dt*fvc::grad(p);
	    U.correctBoundaryConditions();
      
	    phi= (fvc::interpolate(U) & mesh.Sf());
      
	    adjustPhi(phi, U, p);
	  
	    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
	      { 
		fvScalarMatrix pEqn
		  (
		   fvm::laplacian(dt, p) == fvc::div(phi)
		   );
	  
		pEqn.setReference(pRefCell, pRefValue);
	  
		if(nonOrth == nNonOrthCorr)
		  {
		    pEqn.solve(mesh.solver("pFinal"));
		  }
		else
		  {
		    pEqn.solve();
		  }
	  
		if (nonOrth == nNonOrthCorr)
		  {
		    phi -= pEqn.flux();
		  }  
	      }

            #include "continuityErrs.H"

	    U -= dt*fvc::grad(p);
	    U.correctBoundaryConditions();
	}       
	
	turbulence->correct();