
[Sponsors] 
June 29, 2010, 15:42 
pimpleFoam Body Force

#1 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
Hoping someone has been down this dark and winding trail before. I've seen some modifications to other codes, for instance icoFoam and channelFoam but not this guy...
I'm a Fortranconvert new FOAM user and am having a good amount of trouble, trying to add a source term to the pimpleFoam solver. I take it it's not as easy as adding a constant to: tmp<fvVectorMatrix> UEqn ( fvm::ddt(U) + fvm::div(phi, U) + turbulence>divDevReff(U) ); in the UEqn.H file at $FOAM_SOLVERS/incompressible/pimpleFoam I'd like to have a driven halfchannel flow with cyclic b.c's all around. I'm trying to get away from using inlet/outlet patches. I was using fixedValue for pressure and pressureInletVelocity for U but wasn't getting the behavior I need. I'm looking for a little less assumption on the code's part (hence the desire for the body force) Thanks a TON for any help! 

June 30, 2010, 02:52 

#2 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
In principle, you can add the forcing term (it is a constant vector, not a scalar constant, since you are formally imposing a constant mean pressure gradient) to the UEqn. There is nothing that prevents you from doing that.
If you want to do it in a more careful way, you might want to include it in the flux instead than in the UEqn, before solving for the pressure equation. The treatment is identical to the one used for the gravity, for example, in bubbleFoam, and it requires some more modifications to the code. Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

June 30, 2010, 17:05 

#3 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
Ah, great; she's off and running. Thanks for the quick reply!
I created an H file (class?) to read a dict file which specifies the body force as a vector and its dimensions. After that, I stuck the force into the UEqn. If you'll pardon my ignorance, what do you mean by including it in the flux? 

June 30, 2010, 17:10 

#4 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
Oh, and why is that preferable to including the term in the UEqn?


June 30, 2010, 17:23 

#5 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Including force terms in the flux is a way to improve the robustness of the algorithm in some cases (for example when the force is related to the density, like for the gravity).
If you write your momentum equation as ddt(U) + div(UU) = div(tau)  grad(p) + F where F is your forcing term, and then you write in semidiscrete form you get A*U = H + F where I put everything I do not want to carry around in H. Now, the predicted velocity is U = H/A + F/A and the flux, obtained interpolating this on the cell faces (indicated by _f) is phi = (H/A)_f . S + (F/A)_f . S In your case F is constant, so phi = (H/A)_f . S + F(1/A)_f . S In OpenFOAM notation, this means that you should add surfaceScalarField rUAf = fvc::interpolate(rUAf) surfaceScalarField phiForcing = rUAf*(F & mesh.Sf()) to your flux, after adjustPhi. At this point however you have to modify the correction step for the velocity too (the flux is fine), consistently with the new velocity predictor: U += fvc::reconstruct(phiForcing)  rUA*fvc::grad(p); On a side note, you probably won't need this in your case, but it is better to know Hopefully I did not introduce errors while typing Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

June 30, 2010, 18:12 

#6  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Quote:
In other words, it is a way to improve stability when enforcing the force balance. Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

July 2, 2010, 21:06 

#7 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
Oh, okay. I think the first post makes sense to me, if i understand the new variables. The second one definitely does.
Thanks again...you're right, I would like to know the correct way to do it rather than a way. Next up is getting my forcing term timedependent..."here be dragons" I'll check back if I can't manage to work something out. 

July 3, 2010, 00:39 

#8  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Quote:
Quote:
Quote:
Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

July 8, 2010, 19:59 

#9 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
I think I'm hot on the trail, I just can't find to seem to find a way to call out the runtime as a dimensioned scalar. I tried creating a new time variable with:
dimensionedScalar timemdj ( timemdj [0 0 1 0 0 0 0 ] runTime.Value() ); to include in my equation for the body force but to no avail... This is getting fun! 

July 9, 2010, 01:40 

#10 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
It is
dimensionedScalar timemdj ( "timemdj", dimensionSet(0, 0, 1, 0, 0), <yourExpression> // Replace with what you need ); or, simply dimensionedScalar timemdj = <yourExpression>; Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 9, 2010, 11:03 

#11 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
I'm confused...
my forcing is described as dimensionedVector waveP=amp*cos(freq*time_); //determining the new forcing and this yields a forcing in the xdirection as I wanted but I realized: 1. I never included any information to make this an x forcing 2. I'm adding a vector to a matrix (in the UEqn, below)??? tmp<fvVectorMatrix> UEqn ( fvm::ddt(U) + fvm::div(phi, U) + turbulence>divDevReff(U)  waveP ); Is there some assumption on OF's part or something? 

September 9, 2010, 12:09 

#12 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
How do you define Amp? It has to be a vector for the code to compile...
Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 9, 2010, 12:15 

#13 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
Okay, I see now...minor freakout.
I don't quite understand how the vector waveP is added to higher dimension arrays like U... 

September 9, 2010, 12:22 

#14  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Quote:
Each component of waveP is added to the corresponding scalar equation for the component of U. So, say you have g= (gx, gy, gz) and you do UEqn += g; You add gx to the equation for Ux, gy to the equation for Uy and so on. OpenFOAM syntax just makes it easier for you to do this in only one step. Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 9, 2010, 14:33 

#15 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
Okay, to check my understanding of the grand scheme of things:
Each node has a vector equation associated with it. So, in effect we have the matrix of equations (UEqn) added to a matrix of the same size, with all elements being the vector waveP correct? If so, I suppose creating a spatially variable forcing would just involve directly specifying a forcing matrix. Am I on the right track here? Thanks a ton! 

September 9, 2010, 15:11 

#16 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,911
Rep Power: 28 
Yes. Put into the code, you simply define a volVectorField, initialize it with the values of your force term, and add it to UEqn.
P.S. You amp is (amp, 0, 0) right?
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541) OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 9, 2010, 15:19 

#17 
New Member
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 9 
Beautiful! Seems like this may be a fairly painless solution. Much simpler than some other schemes I had rolling around in my head.
Yep, just without the commas It's read from a wavePressDict like viscosity is. 

October 7, 2010, 10:27 

#18 
Member
Join Date: Apr 2009
Location: Karlsruhe, Germany
Posts: 96
Rep Power: 10 
Hi,
I'm trying to implement new anisotropic turbulence models in OpenFOAM. It's like a normal isotropic turbulencemodel, with an explicit nonlinear part too. So the Reynoldsfluxes are defined as: R = linearStress + nonlinearStress (explicit) so divDevReff=  fvm::laplacian(nuEff(), U)  fvc::div(nuEff()*dev(fvc::grad(U)().T())) + fvc::div(nonlinearStress_) Would it make here sense to put the nonlinearpart in the fluxpart as descirbed above? Would the results be better than putting all of them in a fvVectorMatrix? Or is this only for the convergence, the robustness of the system better? There are some nonlinear models in OF1.7 like the nonlinearshih model, but here the nonlinearStress_term is in the fvVectorMatrix, too. Using them to solve a fully developed channel flow I have velocities out of the wall. Could it be the reason, because the nonlinearStresspart is described in the fvVectorMatrix? Thanks and regards, Thomas 

July 3, 2014, 07:16 

#19 
Senior Member
M. Montero
Join Date: Mar 2009
Location: Madrid
Posts: 125
Rep Power: 10 
Hi Alberto,
In case the force is managed in the first way ( force term in the flux), my question is if Ueqn must be modified or not. I am working over the rotorDiskSource (OF2.3.0) and some pressure wiggles happen when discreted body forces are included directly over the Ueqn. I think that your approach could improve this situation. http://www.openfoam.org/mantisbt/view.php?id=1313 A*U = H + F I think it should be: A*U = H  grad[p] + F U=H/A  1/A*grad[p] + F/A Excuse me for using the same topic but there are other topics where others have reported this bug and no answer have been given. http://www.cfdonline.com/Forums/ope...rcemodel.html http://www.cfdonline.com/Forums/ope...lwiggles.html I am trying to implement the body force in this other way and I am not able to do it correctly. Could you please help me? Last edited by be_inspired; July 3, 2014 at 10:20. 

July 10, 2014, 07:04 
urgent help for body force

#20  
New Member
Manoj Jaiswal
Join Date: Jun 2014
Posts: 3
Rep Power: 5 
Hi i am a new Foamer and i am simulating flow over a cylinder(in pimplefoam) . i want to add a body force in the flow , so what are the steps needed for the same . I have gone through a tutorial of how to add temperature in icoFoam (http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam) but still not clear about the addition of body force.
Thanks Quote:


Tags 
body force, pimplefoam, source term 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Body force at the cell face  Souviktor  Fluent UDF and Scheme Programming  0  March 31, 2009 08:54 
DPM Body Force  Sandilya Garimella  FLUENT  1  April 8, 2008 03:30 
Problems with SUPG body force term  FEM question  Main CFD Forum  0  January 21, 2006 18:51 
how to include body force in cfx  selvam R.P  CFX  4  November 25, 2005 05:01 
UDF for lift force on a bluff body  sawa  FLUENT  2  April 11, 2005 03:06 