CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   PISO algorithm for particular NS equation (special convection term) (http://www.cfd-online.com/Forums/openfoam-programming-development/90916-piso-algorithm-particular-ns-equation-special-convection-term.html)

Cyp July 25, 2011 04:19

PISO algorithm for particular NS equation (special convection term)
 
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

santiagomarquezd July 25, 2011 21:48

Cyp, this term is a momentum source, take a look to post #9 from Alberto in this thread.

Regards.

Cyp July 26, 2011 03:37

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

levka May 14, 2012 10:44

convection term
 
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)

Cyp May 16, 2012 08:57

Quote:

Originally Posted by levka (Post 360983)
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

levka May 16, 2012 09:06

...
 
Thanks a lot!
I will give a try and come back.

Quote:

Originally Posted by Cyp (Post 361417)
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 May 20, 2012 08:26

...
 
1 Attachment(s)
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?

Cyp May 20, 2012 16:55

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

levka June 12, 2012 09:36

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



All times are GMT -4. The time now is 02:41.