[jaswi]

Dear All

Greetings.

Conceptually both SIMPLE and PISO formulations can be used to solve problems of transient nature. In the list of OpenFOAM solvers

--> icoFoam is the transient solver for incompressible, laminar flow of Newtonian fluids.

--> transientSimpleFoam is the transient solver for incompressible, turbulent flow of on-Newtonian
fluids.

In principle if one modifies the transientSimpleFoam solver by removing the turbulence calculation part then the two solvers should deliver the same result.

As far my as understanding goes the transientSimpleFoam shall take more time to converge as it uses under-relaxation.

The original icoFoam solver in its unmodified form is:
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)

{

Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readPISOControls.H"
# include "CourantNo.H"
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

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

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

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

adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(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);
}

Now here is the transientSimpleFoam with the changes. I have just removed the turbulence calculation part and the fvMatrixMatrix UEqn is same as in icoFoam. Here is the changed version:

int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;

# include "readPISOControls.H"
# include "CourantNo.H"
// Pressure-velocity SIMPLE corrector loop
for (int corr=0; corr<nCorr; corr++)
{
// Momentum predictor

tmp<fvvectormatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ fvm::laplacian(nu, U)
);

UEqn().relax();

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

p.boundaryField().updateCoeffs();
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

// Store pressure for under-relaxation
p.storePrevIter();

// Non-orthogonal pressure corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);

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

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
# include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
// Momentum corrector
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return(0);
}


I ran both solvers on the cavity with same settings in the controlDict, fvSchemes, fvSolution. the fvSolution has the subDict settings for both algorithms:

PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}

// Relaxation factors are used for transient SIMPLE
relaxationFactors
{
p 0.5;
U 0.8;
}

The results I get are completely different. I am posting here the pressure and velocity fields at t=0.5 for both solvers.

IcoFoam results (all the settings are standard by which I mean the settings given in the UserGuide for the lid driven cavity tutorial)






transientSimpleFoam results





Note the difference in the direction of velocity vectors


The difference in the results could be due to several reason. All I am trying here is to understand how to tune these two solvers which are solving the same physics but using different pressure-velocity coupling.

Any inputs will be helpful

With best regards
Jaswi

[ngj]

Hi Jaswi

Looking at the lowermost plot, which I suspect is the velocity field, I cannot see any boundary, where the velocity is not zero!?!
Thus could it be that you have forgotten to set the boundary condition on U?

Best regards,

Niels

[jaswi]

Hi Niels

Thank for the comment.

Please take a notice of the comments in the post. I have mentioned above that the first two pictures, i.e. pressure and velocity plots are for the led driven cavity using icoFOAM

Now icoFoam uses PISO pressure-velocity coupling algorithm. There is also another solver called transientSimpleFoam which uses the SIMPLE pressure-velocity coupling algorithm . The last two picture are pressure and velocity for the same led driven cavity but now using the modifeid transientSimpleFoam given above in the post.

I again emphasize on the fact that the two solvers
differing only in their pressure-velocity coupling algorithm were run on the same case i.e. lid driven cavity.

Hope that clears the doubt. I am waiting eagerly for the comments :-)

With Regards
Jaswi

[mattijs]

Any reason you're adding and not subtracting the viscous contribution in your UEqn?

+ fvm::laplacian(nu, U)

[jaswi]

Hi Mattijs

Thanks a lot for pointing out that mistake.

I am really sorry for not cross checking my stuff and creating unnecessary trouble.

With that correction , both SIMPLE and PISO deliver identical results , at first glance at least.

I am thankful to you for help. If possible please point me to necessary routines on how I can compare the error for the two solvers and make a quantitative analysis. Just some pointers.

In that way I can do the quantitative analysis for these two pres-velocity algorithms solving the same physical model and put the results on Wiki. It might be useful to somebody trying to understand the solution approach

With Best Regards
Jaswi

[ngj]

Hi

Nice spotting Mattjis...
That could easily explain the lowermost picture, because the change in sign makes the diffusion term into an energy source. Thus if you put on a scale bar, you will probably experience that the maximum velocities in your transientSimpleFoam case is significantly larger than the lid velocity.

- Niels

[bernarde]

Hi

I'm looking for the transientSimpleFoam solver, but I cannot find it in the 1.5 distribution. Is there another location for it?

Thanks