
[Sponsors] 
June 10, 2009, 16:16 
Add a source term to simpleFoam

#1 
Senior Member
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 13 
Hello All,
I want to add a force vectorField as body force to the equations in simpleFoam. These forces do not change in time, it will simply be a vectorField with most cells set at zero and some will have a nonzero force vector. Is it just as simple as putting: >>> eqnResidual = solve ( UEqn() == fvc::grad(p) (new part: + source vector) ).initialResidual(); <<< Probably there is more involved. I searched this forum and surprisingly, similar questions were posted a couple of times but none of them were really answered. Hopefully this time some expert passes by? Brgds, Mark 

June 10, 2009, 21:09 

#2 
Senior Member
Pawel Sosnowski
Join Date: Mar 2009
Location: Munich, Germany
Posts: 105
Rep Power: 14 
it is as simple as that
you just have to make sure that the units match. the following lines are just because I may be tired, so if you find those advice to direct, please forgive me lest say you want to read this field from "constant" folder. in createFields.H you just add: Code:
volVectorField my_field_name ( IOobject( "my_field", runTime.constant(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE ), mesh ); Code:
eqnResidual = solve ( UEqn() == fvc::grad(p) + my_field ).initialResidual(); you can do it in the code Code:
eqnResidual = solve ( UEqn() == fvc::grad(p) + my_field/rho ).initialResidual(); hope it helps you a bit. Pawel 

June 11, 2009, 04:40 

#3 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,869
Rep Power: 33 
Hi Mark and Pawel
I do not agree that it is as simple as that, since if "momentumPrediction" is set to "no", then the momentum equation is not solved, and hence the body force is not taken into consideration. Thus my experience is that in the pEqn.H one need to modify the momentum corrector as follows: U = (fvc::grad(p)  my_field / rho )/AU; However, you might need to consider how this contribution is added. I have used a constant in space, varying in time body force, thus the divergence of the momentum equation gives a 0 (zero) contribution to the pressure equation. As your bodyforce in varying in space, the term will take all the way to the pressure equation, I am not certain that it is sufficient to do as above. Best regards, Niels 

June 11, 2009, 06:14 

#4 
Senior Member
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 13 
Hello Pawel and Niels,
First of all thanks for your contributions. Niels, you explain about momentumcorrector but we are talking about simpleFoam, which always uses a momentum corrector, right? In that case, does it still make sense? Next: now I modified the simpleFoam solver as described by both of you. I tried it on the pitzDaily example (a backward facing step) and add some force on the upperwall, directing downwards. The liquid is not moving downward at all. While writing this I realise this might be caused by the fact that the fluid is incompressible. Does exerting a force on a wall makes sense? I better take some time to test it with a bodyForce on the internal fluid (what is what I actually want to do anyway). Finally dimensions: Pawel, you are right. I need to supply force per unit volume and divided by rho, which gives: F/V/rho = (kg m s2) * (m3) * (kg m3)^1 = m/s2. Once again thanks for your help, Mark 

June 11, 2009, 07:02 

#5 
Senior Member
Pawel Sosnowski
Join Date: Mar 2009
Location: Munich, Germany
Posts: 105
Rep Power: 14 
Hello Niles and Mark,
Niles, I assumed that "momentumPredictor" is set to "true". Now I have a question if we assume that we solve the momentum equation, and that we add, for example, a gravity mass force, we get the UEqn like: Code:
tmp<fvVectorMatrix> UEqn ( fvm::div(phi, U) + turbulence>divDevReff(U) ); Code:
eqnResidual = solve ( UEqn() == fvc::grad(p) + GravForce ).initialResidual(); Code:
volScalarField AU = UEqn().A(); Mark, first check the boundary type it should be "fixedValue" (I found myself forgeting about it several times) As for the second thing I belive that by adding this additional field you simply introduse a velocity source term to the UEqn. For me it seems that by putting some value on the wall would be just like applying a velocity inlet on that wall. (But it would be better if someone veryfied my statement). Regards, Pawel 

June 12, 2009, 05:29 

#6 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,869
Rep Power: 33 
Hi Pawel
In the pressure equation, U is estimated as U = UEqn().H()/AU; hence based on all the terms in the UEqn, thus not on neither grad(p) nor my_fields. This means that the pressure correction to U does not include any of these, and therefore you need to have it in the second term as well. One of my colleagues has been calculating laminar and turbulent oscillating wave boundary layers using OpenFOAM and he achieve perfect comparison with experimental/theoretical data using the method I have described above. Best regards, Niels 

June 12, 2009, 06:39 

#7 
Senior Member
Pawel Sosnowski
Join Date: Mar 2009
Location: Munich, Germany
Posts: 105
Rep Power: 14 
Thenks Niles for pointing that out.
Fortunately I added my source terms into velocity equation so they were included while pressure calculations. Now I will pay more attention to that matter. Regards, Pawel 

June 12, 2009, 06:54 

#8 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,869
Rep Power: 33 
Ah, that explains why you didn't need to add it in the bottom.
Good luck to both of you and enjoy the upcoming weekend, Niels 

September 9, 2009, 06:50 

#9  
New Member
nick stoppelkamp
Join Date: Jul 2009
Posts: 4
Rep Power: 13 
Hi,
I am also trying to bring in a source term in my simpleFoam solver, following the code snipped below. Quote:
But unfortunately the compiler complains that rho is not defined. So what do I have to do to get thr rho field? Thanks and best regards, Nick 

September 10, 2009, 11:55 

#10 
Senior Member
Pawel Sosnowski
Join Date: Mar 2009
Location: Munich, Germany
Posts: 105
Rep Power: 14 
Hello Nick!
The way I see it, there are two ways of including rho in your case: 1) when you create your field, divide it by rho at the beginning, set the dimensions as [myField/rho] and solve just: Code:
UEqn() == fvc::grad(p) + myFieldDivRho Code:
volScalarField rho ( IOobject ( "rho", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), <initial_scalar_value_of_rho> ); Hope it helped a bit best regards, Pawel 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
momentum source term  zwdi  FLUENT  14  June 27, 2017 16:40 
How to add a flux source to SIMPLEfoam  ehsan_vaghefi  OpenFOAM Running, Solving & CFD  0  September 29, 2008 22:19 
DxFoam reader update  hjasak  OpenFOAM PostProcessing  69  April 24, 2008 02:24 
DecomposePar links against liblamso0 with OpenMPI  jens_klostermann  OpenFOAM Bugs  11  June 28, 2007 18:51 
Add source term in species equation  zhou1  FLUENT  1  October 21, 2003 07:28 