CFD Online Discussion Forums - problem in writing a solver for a 3-step runge kutta(incompressible flow)

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - problem in writing a solver for a 3-step runge kutta(incompressible flow) (https://www.cfd-online.com/Forums/openfoam-programming-development/130340-problem-writing-solver-3-step-runge-kutta-incompressible-flow.html)

problem in writing a solver for a 3-step runge kutta(incompressible flow)

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 :

Code:

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();

}

}

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

Code:

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();

}

}

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?