CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

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

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   February 24, 2014, 14:03
Default problem in writing a solver for a 3-step runge kutta(incompressible flow)
  #1
ooo
Member
 
Join Date: Feb 2012
Posts: 49
Rep Power: 6
ooo is on a distinguished road
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();
}
}
ooo is offline   Reply With Quote

Old   March 4, 2014, 12:19
Default
  #2
ooo
Member
 
Join Date: Feb 2012
Posts: 49
Rep Power: 6
ooo is on a distinguished road
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();
}
}
ooo is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Extrusion with OpenFoam problem No. Iterations 0 Lord Kelvin OpenFOAM Running, Solving & CFD 8 March 28, 2016 11:08
How to write k and epsilon before the abnormal end xiuying OpenFOAM Running, Solving & CFD 8 August 27, 2013 15:33
Maximum words supported by the flow solver exceeded rae CFX 3 June 6, 2013 06:49
Problem with changing mass flow rate roo FLUENT 3 June 4, 2013 14:53
On Setring Step Size in unsteady Solver of FLUENT lzgwhy FLUENT 0 August 19, 2009 22:44


All times are GMT -4. The time now is 16:28.