[ooo]

I'm implementing a 3-step runge kutta descritization for solving incompressible flow.
The results are completely wrong.
I would appreciate if you have a look at the code below and tell me if you see any problem :

while (runTime.loop())
    {

 double alpha,gamma,zeta;

 volVectorField Unew(U), UOld2(U), UTmp(U) ;
 surfaceScalarField phiOld2(phi), phiTmp(phi);
 volScalarField pPhi(p);

        for (int i = 1 ;  i<=3 ; ++i)
        {
            switch (i)
            {
                case 1 : {UOld2=U ; phiOld2=phi ; alpha=4./15. ; gamma= 8./15. ; zeta =0.; break ; }
                case 2 : {UTmp=U ; phiTmp=phi ; alpha=1./15. ; gamma= 5./12. ; zeta =-17./60.; break ; }
                case 3 : {UOld2=UTmp ; phiOld2=phiTmp ; alpha=1./6. ; gamma= 3./4. ; zeta =-5./12.; break ; }
                default : break ;
            }

        adjustPhi(phi, U, p);

        U = U + runTime.deltaT() *
        (
         + 2*alpha*fvc::laplacian(nu,U)
         - 2*alpha*fvc::grad(p)
         - gamma*fvc::div(phi, U)
         - zeta*fvc::div(phiOld2, UOld2)
        );




      solve(alpha*runTime.deltaT()*fvm::laplacian(nu,Unew) - fvm::Sp(1.,Unew)  == //(alpha*runTime.deltaT())  ==
      (-1.)*(U) + alpha*runTime.deltaT()*fvc::laplacian(nu,U) );

      solve (fvm::laplacian(pPhi) == fvc::div(U)/(2.*alpha*runTime.deltaT()));

      U =  Unew - (2.*alpha*runTime.deltaT()*fvc::grad(pPhi));
      p += pPhi - alpha*runTime.deltaT()*nu*(fvc::laplacian(pPhi));

         adjustPhi(phi, U, p);
         U.correctBoundaryConditions();

}
runTime.write();
}
}

[ooo]

I make my code shorter, maybe someone has any idea regarding any problem on this kind of implementation :

while (runTime.loop())
{
for (int i = 1 ;  i<=3 ; ++i)
{
U = U + runTime.deltaT() *
        (
         + 2*alpha*fvc::laplacian(nu,U)
         - 2*alpha*fvc::grad(p)
         - gamma*fvc::div(phi, U)
         - zeta*fvc::div(phiOld2, UOld2)
        );

      solve(alpha*runTime.deltaT()*fvm::laplacian(nu,Unew) - fvm::Sp(1.,Unew)  == //(alpha*runTime.deltaT())  ==
      (-1.)*(U) + alpha*runTime.deltaT()*fvc::laplacian(nu,U) );

      solve (fvm::laplacian(pPhi) == fvc::div(U)/(2.*alpha*runTime.deltaT()));// pPhi is a pseudo pressure without physical meaning

      U =  Unew - (2.*alpha*runTime.deltaT()*fvc::grad(pPhi));
      p += pPhi - alpha*runTime.deltaT()*nu*(fvc::laplacian(pPhi));

      adjustPhi(phi, U, p);
      U.correctBoundaryConditions();
}
}

[hua1015]

Dear Foamer,
Recently, I am interested in writing a solver for a 3-step runge kutta for incompressible flow and suffered some problem. Did you implemented it? Would you mind give me some suggestions?