So, I've been working on trying to do some simulations of basic steady laminar flow with a body force term. The force term is supplied as a constant field "f", and that part's working well. The rest of the code is derived from simpleFoam, with the viscous term simplified to a Newtonian model.

The problem is that, instead of converging to zero properly, the residuals in the momentum and pressure equations go to a constant value and sit there. If I change the relaxation factor on the velocity equation, the final values of the residuals and the actual velocity and pressure fields change a bit.

I'm pretty sure this means that I didn't include the force term in the right place, and so the system goes to a fixed point where the velocity equation puts in a correction and the pressure equation takes it back out, or something like that.... Unfortunately, I'm not having any luck comparing this to things like boundaryFoam and trying to understand what I need to change in my code to fix it. Can anyone look at this and see what I need to change?

Here's the core of the solver code:

-----------------

****for*(runTime++;*!runTime.end();*runTime++)
****{
********Info<<*"Time*=*"*<<*runTime.timeName()*<<* nl*<<*endl;

#*******include*"readSIMPLEControls.H"

********p.storePrevIter();

********//*Pressure-velocity*SIMPLE*corrector
********{
************//*Momentum*predictor

************fvVectorMatrix*UEqn
************(
****************fvm::div(phi,*U)
**************-*fvm::laplacian(nu,*U)
************);

************UEqn.relax();

************solve(UEqn*==*-fvc::grad(p)*+*f);**//*Note: this is where the force is added.

************volScalarField*AU*=*UEqn.A();
************U*=*UEqn.H()/AU;

************phi*=*fvc::interpolate(U)*&*mesh.Sf();
************adjustPhi(phi,*U,*p);

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

****************fvScalarMatrix::reference*pRef*=
********************pEqn.setReference(pRefCell,*pR efValue);
****************pEqn.solve();
****************pEqn.unsetReference(pRef);

****************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();

********Info<<*"ExecutionTime*=*"
************<<*runTime.elapsedCpuTime()
************<<*"*s\n\n"*<<*endl;
****}

-----------------

Any suggestions would be much appreciated. Thanks very much!

[mprinkey]

Where you've added the body force term, it is not included in the Jacobi update that generates the predicted velocities (U = UEqn.H()/AU;) so it is going to be missing from the pressure equation.

I think that you either need to put the body force term into the LHS of UEqn:

fvVectorMatrix UEqn
(
fvm::div(phi, U)
- fvm::laplacian(nu, U)
- f
);

That will make sure that the UEqn.H() evaluation includes the body force effects and that will make sure that the pressure effects are present in the predictor.

Or if you exclude the body force term from the Jacobi updated velocities, you need to include it in the pressure equation in an appropriate way as is done in some of the other solvers. That way the velocity corrector step will include the body force effect.

Thanks! A nearly trivial change, which would have probably taken me quite a while to figure out -- somehow I had the impression that that it would get taken into account in the solve() operation regardless. Obviously that was wrong.

In any case, I used the first of the two fixes you suggested, and it works beautifully.

And, with that, the toolchain that I've been building over the last six months of my research is finally debugged enough to be a complete end-to-end chain, just in time for my meeting for my advisor tomorrow. Thank you!