# SIMPLE force term isnbt converging what have I missed

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 LinkBack Thread Tools Display Modes
 July 17, 2006, 19:21 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!

 July 17, 2006, 19:38 Where you've added the body fo #2 Senior Member   Michael Prinkey Join Date: Mar 2009 Location: Pittsburgh PA Posts: 352 Rep Power: 17 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.

 July 17, 2006, 20:25 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!

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post waynezw0618 OpenFOAM Running, Solving & CFD 49 April 8, 2008 04:23 Mike Main CFD Forum 1 September 12, 2007 12:08 Rick FLUENT 0 March 9, 2006 11:28 FEM question Main CFD Forum 0 January 21, 2006 18:51 R.Sureshkumar Main CFD Forum 6 February 23, 1999 06:33

All times are GMT -4. The time now is 16:05.

 Contact Us - CFD Online - Privacy Statement - Top