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

PISO algorithm for particular NS equation (special convection term)

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

Reply
 
LinkBack Thread Tools Display Modes
Old   July 25, 2011, 04:19
Default PISO algorithm for particular NS equation (special convection term)
  #1
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 229
Rep Power: 9
Cyp is on a distinguished road
Hi !

I would like to solve a special "Naviers-Stokes" equation. It consists in a classical NS equation with an additionnal source term and a special convection term. Indeed, in the latter the velocity we look for solving is transported by a constant velocity.

The final system can be formulated as follow:

\frac{\partial \rho \textbf{U}}{\partial t} + \nabla \cdot (\rho \textbf{U}_{0} \textbf{U}) = -\nabla p + \nabla \cdot (\nu \nabla U) - \nu \textbf{e}_{0}

where e0 is a constant vector and U0 is a constant velocity field

and of course the continuity equation for incompressible fluid:
\nabla \cdot (\rho \textbf{U}) = 0

Moreover, the density and the viscosity are supposed to be constant.

I tried the following PISO loop, but I not sure the PISO algorithm described in Jasak thesis is suitable to my equations :

Quote:

surfaceScalarField phi0 = rho*(fvc::interpolate(U0) & mesh.Sf());

while (runTime.run())
{

#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"

runTime++;

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

fvVectorMatrix UEqn
(
fvm::ddt(rho,U)
+ fvm::div(phi0, U)
- fvm::laplacian(mu, U)
+ UnitE*mu*vector(0,1,0)
);

solve(UEqn == -fvc::grad(p));

// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A();

U = rUA*UEqn.H();
surfaceScalarField phiRho = phi/rho;

phi = rho*
(
(fvc::interpolate(U) & mesh.Sf())) + rho*rho*(fvc::ddtPhiCorr(rUA, U, phiRho)
);

adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rho*rUA, p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

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

#include "continuityErrs.H"

U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
Specially I would like confirmation of the treatment of the pressure equation and the ddtPhiCorr correction.

Thank you for your help,
Cyp
Cyp is offline   Reply With Quote

Old   July 25, 2011, 21:48
Default
  #2
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Cyp, this term is a momentum source, take a look to post #9 from Alberto in this thread.

Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   July 26, 2011, 03:37
Default
  #3
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 229
Rep Power: 9
Cyp is on a distinguished road
Thank you for your answer.

Indeed the treatment of a momentum source term is a topic widely discussed. My main question deals with the convection term \nabla \cdot \rho\textbf{U}_{0}\textbf{U} when \textbf{U}_{0} is a constant velocity different than the variable \textbf{U} we looked to solve.

I am not sure the classical process with phi0 instead of phi in the convection term as I mentioned in the quoted code above is correct ?

Regards,
Cyp
Cyp is offline   Reply With Quote

Old   May 14, 2012, 10:44
Default convection term
  #4
Member
 
Lev
Join Date: Dec 2010
Posts: 31
Rep Power: 6
levka is on a distinguished road
Hello, Cyp
did you find any confirmation on topic how to treat convection term with some known U0 ?

I have similar task, but without source term. Do i have to modify standard PISO loop from IcoFoam?

In my case i have:
+ fvm::div(phi, U)
+ fvm::div(phi0, U)
+ fvm::div(phi, U0)
levka is offline   Reply With Quote

Old   May 16, 2012, 08:57
Default
  #5
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 229
Rep Power: 9
Cyp is on a distinguished road
Quote:
Originally Posted by levka View Post
Hello, Cyp
did you find any confirmation on topic how to treat convection term with some known U0 ?

I have similar task, but without source term. Do i have to modify standard PISO loop from IcoFoam?

In my case i have:
+ fvm::div(phi, U)
+ fvm::div(phi0, U)
+ fvm::div(phi, U0)

Hi levka,

If you want to add fvm::div(phi0, U) in the Naviers-Stokes equations you do not have to change the PISO loop.

However, the last term you mentioned cannot be treated implicitly (I mean, using fvm::div(...)). Indeed the unknown of this equation is U. When you use fvm::blahblah , you build the linear matrix system where U is the unknown. In your equation, fvm::div(phi,U0) means that you want to build the matrix for U0 within the matrix for U... It doesn't make sens. I suggest you consider this term as an explicit source term using fvc::div(phi,U0). Once again, in this case, the PISO loop doesn't change.

Best,
Cyp
Cyp is offline   Reply With Quote

Old   May 16, 2012, 09:06
Default ...
  #6
Member
 
Lev
Join Date: Dec 2010
Posts: 31
Rep Power: 6
levka is on a distinguished road
Thanks a lot!
I will give a try and come back.

Quote:
Originally Posted by Cyp View Post
Hi levka,

If you want to add fvm::div(phi0, U) in the Naviers-Stokes equations you do not have to change the PISO loop.

However, the last term you mentioned cannot be treated implicitly (I mean, using fvm::div(...)). Indeed the unknown of this equation is U. When you use fvm::blahblah , you build the linear matrix system where U is the unknown. In your equation, fvm::div(phi,U0) means that you want to build the matrix for U0 within the matrix for U... It doesn't make sens. I suggest you consider this term as an explicit source term using fvc::div(phi,U0). Once again, in this case, the PISO loop doesn't change.

Best,
Cyp
levka is offline   Reply With Quote

Old   May 20, 2012, 08:26
Default ...
  #7
Member
 
Lev
Join Date: Dec 2010
Posts: 31
Rep Power: 6
levka is on a distinguished road
Hello, Cyp
Have a look at the attached picture (that is a probe of U in time). I computed the following:


fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
+ fvm::div(phi0, U)
);

solve(UEqn == -fvc::grad(p)- fvc::div(phi, U0));


Where U0- constant velocity profile, phi0- flux based on that U0. U(initially)- profile of full turbulent flow at Re=5000.

The result i obtained is not expected: flow became absolutely laminar...
I suppose to see turbulent flow, according to the math equations

Do you have any idea?
Attached Images
File Type: jpg look2.jpg (55.7 KB, 21 views)
levka is offline   Reply With Quote

Old   May 20, 2012, 16:55
Default
  #8
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 229
Rep Power: 9
Cyp is on a distinguished road
If you consider the additional source term as an argument of the solve function, you have to modified the reconstruction of the velocity (see the link mentioned in the second post of this topic).

You should try :

Code:
fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          + fvm::div(phi0, U)
 + fvc::div(phi, U0)
          );

        solve(UEqn == -fvc::grad(p));
Regards,
Cyp
Cyp is offline   Reply With Quote

Old   June 12, 2012, 09:36
Default ...
  #9
Member
 
Lev
Join Date: Dec 2010
Posts: 31
Rep Power: 6
levka is on a distinguished road
Cyp,
finally i found that your idea is working, thanks!

Code:
fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          + fvm::div(phi0, U)
          + fvc::div(phi, U0)
          );

        solve(UEqn == -fvc::grad(p));
levka 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
source term units - poission equation johnwinter Fluent UDF and Scheme Programming 0 June 6, 2011 00:53
Estimate convection and diffusion terms in transport equation Lance CFX 0 April 4, 2011 11:16
Pressure work term in turbulent K.E. equation lost.identity Main CFD Forum 0 March 8, 2011 13:21
Modified Equation for CFX algorithm Craig Johansen CFX 0 August 27, 2004 23:02
bouyancy term in epsilon equation Michael Main CFD Forum 1 June 25, 1999 10:20


All times are GMT -4. The time now is 10:20.