[nikwin]

Hello,

I need to compute the spatial terms of the momentum eq. in rhoPisoFoam. As a check I added a part to rhoPisoFoam which computes all the terms of the momentum eq. to see that it is fulfilled. However accoding to my calculations it's not, even though the residual printed during runTime is below tolerance. Any suggestion of what could be wrong?

Thanks
/NW

int main(int argc, char *argv[])
{
#include "setRootCase.H"

#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{

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

runTime++;

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

#include "rhoEqn.H"
#include "UEqn.H"

// --- PISO loop
for (int corr=1; corr<=nCorr; corr++)
{
#include "hEqn.H"
#include "pEqn.H"
}

turbulence->correct();
rho = thermo.rho();

Info << "START CHECK AFTER rho=thermo.rho()" << endl;

volVectorField UEqn4(IOobject("UEqn4",runTime.timeName(),mesh,IOo bject::NO_READ,IOobject::AUTO_WRITE),mesh,dimensio nedVector("UEqn4",dimless,vector::zero));

fvVectorMatrix UEqn4Matrix
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
+ fvc::grad(p)
);

vectorField UEqn4Residual(UEqn4Matrix.residual());
forAll (mesh.cells(),cellI)
{
UEqn4[cellI] = -UEqn4Residual[cellI]/mesh.V()[cellI];
}

vector mean(sum(UEqn4.internalField()*mesh.V())/sum(mesh.V()).value());

Info << "max: " << max(UEqn4.internalField()) << endl;
Info << "min: " << min(UEqn4.internalField()) << endl;
Info << " mean: " << mean << endl;
Info << "variance.x: " << sum(sqr(mean.x()-UEqn4.internalField().component(0)))/mesh.nCells() << endl;
Info << "variance.y: " << sum(sqr(mean.y()-UEqn4.internalField().component(1)))/mesh.nCells() << endl;
Info << "variance.z: " << sum(sqr(mean.z()-UEqn4.internalField().component(2)))/mesh.nCells() << endl;

Info << "END CHECK AFTER rho=thermo.rho()" << endl;

runTime.write();

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

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

return 0;
}

Results from mesh with 2e6 cells, pipe flow.

max: (6181064.994 6807506.826 6688472.669)
min: (-6694996.468 -7079812.026 -6947027.568)
mean: (-77.21956025 27.10909829 -53.67802351)
variance.x: 3019294370
variance.y: 2737483314
variance.z: 2314573170