
[Sponsors] 
July 17, 2006, 19:21 
So, I've been working on tryin

#1 
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(); ********//*Pressurevelocity*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); ************//*Nonorthogonal*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: 176
Rep Power: 13 
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 
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 endtoend chain, just in time for my meeting for my advisor tomorrow. Thank you! 

Thread Tools  
Display Modes  


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 11:28 
Problems with SUPG body force term  FEM question  Main CFD Forum  0  January 21, 2006 18:51 
Source Term for PressureSIMPLE  R.Sureshkumar  Main CFD Forum  6  February 23, 1999 06:33 