Steady incompressible laminar flow using simpleFoam

yongshenglian

Someone posted a message about using simpleFoam to solve steady incompressible flow problem. I just cannot find it anymore.

After reading some papers and posted messages here, I made a new solver based on simpleFoam. It took me a while to compile the solver but finally I have the executable file. However, the solution is always wrong. The velocity is more than 100 times of the inlet velocity over a channel.

I checked my code and compare with simpleFoam. It does not seem that the code is wrong. I also check the boundary condition and initial condition. IF I run PISO, it gives good result.

I post the code here. Hopefully someone can correct me.

Thanks.

*************************

#include "fvCFD.H"
//#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
//#include "incompressible/RASModel/RASModel.H"

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

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 "readSIMPLEControls.H"
# include "initConvergenceCheck.H"

p.storePrevIter();

// Solve the Momentum equation

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

UEqn().relax();

eqnResidual = solve
(
UEqn() == -fvc::grad(p)
).initialResidual();

maxResidual = max(eqnResidual, maxResidual);

// pressure correction
p.boundaryField().updateCoeffs();

volScalarField AU = UEqn().A();
U = UEqn().H()/AU;
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

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

pEqn.setReference(pRefCell, pRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pEqn.solve();
}

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

# include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();

// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();


runTime.write();
# include "convergenceCheck.H"
}

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

return(0);
}


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