CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

SIMPLE force term isnbt converging what have I missed

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   July 17, 2006, 19:21
Default So, I've been working on tryin
  #1
brooksmoses
Guest
 
Posts: n/a
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!
  Reply With Quote

Old   July 17, 2006, 19:38
Default Where you've added the body fo
  #2
Member
 
Michael Prinkey
Join Date: Mar 2009
Posts: 44
Rep Power: 7
mprinkey is on a distinguished road
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.
mprinkey is offline   Reply With Quote

Old   July 17, 2006, 20:25
Default Thanks! A nearly trivial chan
  #3
brooksmoses
Guest
 
Posts: n/a
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!
  Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
MRFZonesC questions what is the mesh_V and why only Coriolis force no centrifugal force waynezw0618 OpenFOAM Running, Solving & CFD 49 April 8, 2008 04:23
Simple rotational force question Mike Main CFD Forum 1 September 12, 2007 12:08
Modelling a force as a momentum source term - HELP Rick FLUENT 0 March 9, 2006 10:28
Problems with SUPG body force term FEM question Main CFD Forum 0 January 21, 2006 17:51
Source Term for Pressure-SIMPLE R.Sureshkumar Main CFD Forum 6 February 23, 1999 05:33


All times are GMT -4. The time now is 13:19.