# How to add a body force to a momentum equation

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

 October 22, 2018, 19:20 How to add a body force to a momentum equation #1 Senior Member   Join Date: Mar 2018 Posts: 115 Rep Power: 8 Hello If I have a constant body force f(fx fy 0) which has a dimension of force/volume/density. I want to add it to the momentum equation, 1) is it possible to use fvOptions in this case? if yes how? don't post links because there no useful info there, I searched the whole internet without success. 2) if it must involve programming, then where in my code (I mean in the pimple loop) should I put the body force? lpz456 likes this.

 October 22, 2018, 19:58 #2 Member   Martin Join Date: Dec 2011 Posts: 40 Rep Power: 14 Cant you include this force by using a modified gravity vector? lpz456 likes this.

October 23, 2018, 08:11
#3
Senior Member

Join Date: Mar 2018
Posts: 115
Rep Power: 8
Quote:
 Originally Posted by msaravia Cant you include this force by using a modified gravity vector?

 October 23, 2018, 08:25 #4 Super Moderator     Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51 This depends highly on what you are trying to implement. For example, you can add a body force directly to the momentum equation or you can put the source term to the pressure. For the second one, my former colleague would be better suited for the answer. Thus, let us consider the simple way: Adding the source term to the momentum equation. Here you have at least two options: a) Modifying the momentum equation inside the code (UEqn.H file). b) Using fvOptions and adding the terms of your needs; e.g., using the codedSource To a) - The file UEqn.H is located in the solver folder you are using. It is recommended to build a new derivation of the underlying solver (there is a lot of information on the web). The path is \$FOAM_SOLVERS/ Adding an arbitrary source term to the momentum of pimpleFoam looks as follows. The common UEqn is: Code: ```tmp tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevReff(U) == fvOptions(U) );``` As you can see, the fvOptions is included as source term. However, adding a new source is as follows: Code: ```tmp tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevReff(U) == fvOptions(U) + mySource );``` It is obvious that mySource is not defined right now. So you have to do it previously. In the following, we add a source term that is acting just in the x-direction: Code: ```const dimensionedVector mySource("mySource", dimensionSet(0,1,-1,0,0,0,0), vector(-4,0,0)); tmp tUEqn . . .``` That´s it. To b) please follow this link and scroll down to the Detailed Description -> Usage https://cpp.openfoam.org/v6/classFoa...21f7bbb656d86d Good luck raj kumar saini, msaravia, kiski and 12 others like this. __________________ Keep foaming, Tobias Holzmann

 October 23, 2018, 09:49 #5 Senior Member   Join Date: Mar 2018 Posts: 115 Rep Power: 8 Thank you very much, dear Tobi, for your answer it is very clear. Actually I am trying to simulate the 2D flow over a NACA airfoil but without using the real geometry, instead, I replace it as a source term. (Lift and drag forces projected to x and y directions). In the UEqn.H, I have a question. Are the following syntaxes the same?: Code: ```tmp tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevReff(U) == fvOptions(U) + mySource );``` Code: ```tmp tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevReff(U) - mySource == fvOptions(U) );``` When the source term is added to the RHS or substracted from the LHS, do they have the same effect or they afffect the numerical stability? 2)When using fvOptions, I guess there is something called "vectorCodedSource", right?

October 23, 2018, 09:58
#6
Member

Martin
Join Date: Dec 2011
Posts: 40
Rep Power: 14
Quote:
 Originally Posted by Evren Linda Please please, if you have no answer do not comment because I don't know what you are talking about.
Ok ok, no more comments that do not imply the direct solution to your problem. Sometimes in science, questions are answers... Anyway, what I meant was:

A body force is an action exerted on matter essentially by: i) gravity, ii) magnetic fields, iii) electric fields, iv) centrifugal and Coriolis accelerations. From the above, gravity is often taken as constant. As you mention that your force is constant; then I assumed that, whatever this force is, you can add it to gravity. For some solvers, gravity is added at RUNTIME as a vector, that is read from a dictionary. You can check the squareBump example in the OpenFOAM tutorials subdirectory.

Thus, your new gravity vector in your input dictionary would look like:

Code:
`g            g           [0 1 -2 0 0 0 0]  (fx fy -9.81);`
or something similar.

This approach would avoid programming a new solver. Of course this implies that your solver allows this body force.

Finally, I think that Tobi's suggestion is the best approach.

 October 23, 2018, 10:52 #7 Super Moderator     Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51 A few comments, 1. I never saw the modifing gravity vector but this sounds to be a real hack! That's excellent, easy and without headache. The only problem is related to the fact that it acts on the whole domain. 2. Out of the box I would expect that both forumations are the same (related to the question of the equations). One can work with Sp operators to make it more stable but that is not the matter of question. 3. There should be a vectorCodedSource Finally,I really like Martins way because it is actually the same what I mentioned but much easier By the way, the units of the source in my example should be wrong msaravia and lpz456 like this. __________________ Keep foaming, Tobias Holzmann

October 23, 2018, 11:36
#8
Senior Member

Join Date: Mar 2018
Posts: 115
Rep Power: 8
Quote:
 Originally Posted by msaravia Ok ok, no more comments that do not imply the direct solution to your problem. Sometimes in science, questions are answers... Anyway, what I meant was: A body force is an action exerted on matter essentially by: i) gravity, ii) magnetic fields, iii) electric fields, iv) centrifugal and Coriolis accelerations. From the above, gravity is often taken as constant. As you mention that your force is constant; then I assumed that, whatever this force is, you can add it to gravity. For some solvers, gravity is added at RUNTIME as a vector, that is read from a dictionary. You can check the squareBump example in the OpenFOAM tutorials subdirectory. Thus, your new gravity vector in your input dictionary would look like: Code: `g g [0 1 -2 0 0 0 0] (fx fy -9.81);` or something similar. This approach would avoid programming a new solver. Of course this implies that your solver allows this body force. Finally, I think that Tobi's suggestion is the best approach.
I am so sorry, I don't mean what you might think. Actually sometimes people think that they are answering but it is very confusing. You might expect the OP to be understanding what you are talking about but unfortunately that's not the case in the most cases. In my opinion, always expect the OP to be absolute beginner and does not know anything and give a detailed answer, this will help the community

Thank you for your answer. but it is not useful in my case, because I have to insert the body force in a specified cellZone only not in the entire domain.

November 3, 2018, 16:26
#9
Retired Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,975
Blog Entries: 45
Rep Power: 128
Greetings to all,

It took me two or three passes through this thread to figure out which side is up... sorry, what I mean is that I wasn't clear if a solution had been reached yet or not. I had to double-check with Evren to know if I should take a better look into this... so I'll write my post as I go along with figuring things out...

I'm having a hard time understanding the original request on the type of force to be applied, so...
Quote:
 Originally Posted by Evren Linda 1) is it possible to use fvOptions in this case? if yes how? don't post links because there no useful info there, I searched the whole internet without success.
... I'll first start with a link to an entry point to a possible solution: the fvOption "buoyancyForce": https://cpp.openfoam.org/v6/classFoa...e.html#details - quoting from there:
Quote:
 Calculates and applies the buoyancy force rho*g to the momentum equation corresponding to the specified velocity field.
It's not the direct solution to the original question, but this is mostly a "there is already something similar, which could possibly be used as a starting point".

I got there from the base "fv :: option" class: https://cpp.openfoam.org/v6/classFoa...1_1option.html - from there, we can find more possible "fvOptions" if we click on the "cellSetOption" one: https://cpp.openfoam.org/v6/classFoa...SetOption.html - this is because these are bound to the possibility to select a list of cells to which the "fvOption" is applied to.

OK, so the basis for adding a new source term is roughly outlined.

Now I have to try and understand the original request:
Quote:
 Originally Posted by Evren Linda If I have a constant body force f(fx fy 0) which has a dimension of force/volume/density. I want to add it to the momentum equation,
So the first detail we need to keep in mind that the direct source term for the U equation has the units of Force already... I was confused the first time I saw this myself, but that's due to how the equation is constructed. Which means that we need to get rid of the volume and density terms... wait... volume and density? "m3" and "kg/m3"? Sooo... if we try to sort out the units, it would simply become "force/kg"?...

OK ok... if I step back, we do have access to the volumes of the cells and possibly access to the density field, so that shouldn't be a problem. But that raises the following question: What volume should be used and which density?
What I mean is: should we simply make the calculation the same for each individual cell? And should it be applied to the whole field?

Because if this is the case, then it should be simply a matter of:
1. Creating a new class copied from the "buoyancyForce" class.
2. And change "rho*g" for the new "(fx fy 0)/V/rho".
3. We can get the "V" field by following the example from here: https://github.com/OpenFOAM/OpenFOAM...tyForce.C#L132 - I found it by looking at the sibling fvOptions that looked like "maybe this one needs volume?"
I guess this could be coded directly into the "fvOptions" dictionary entry with a "vectorCodedSource", so it wouldn't be necessary to have to deal with creating a complete new library for just a custom fvOption... the example at https://cpp.openfoam.org/v6/classFoa...e.html#details - seems already ready to be used as a good template for this...

But if this is not the desired math operation that is meant to be done - for example, if the idea is to impose only a force directly on a patch - then it's a whole other ball game...

My last concern is whether this a standard feature that is needed for airfoil studies? Because if it is, that I or Tobi or anyone else can look into creating a definitive fvOption for this and contribute it to the OpenFOAM Foundation, so that it's always available in future versions.

Best regards,
Bruno
__________________

 November 29, 2018, 09:03 #10 New Member   Mathias Poulsen Join Date: Feb 2018 Location: Denmark Posts: 9 Rep Power: 8 Hi. I do not know if this is still relevant but to you but here goes nothing. I am not an expert in this, but you could try using fvOptions as you suggested. The following code is the fvOptions file which should be placed in the constant folder of your case. I am not completely sure about the units and expression for the source, but you can try changing it so that is fits your requirements. Code: ```/*--------------------------------*- C++ -*----------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Version: 6 \\/ M anipulation | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object fvOptions; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // codedSource { type vectorCodedSource; selectionMode cellZone; cellZone BANANA; fields (U); name codedSource; codeAddSup #{ const label cellZoneID = mesh().cellZones().findZoneID("BANANA"); // find ID for the cellZone "BANANA" const cellZone& zone = mesh().cellZones()[cellZoneID]; const cellZoneMesh& zoneMesh = zone.zoneMesh(); const labelList& cellsZone = zoneMesh[cellZoneID]; //list of all cellZone cell ID's const scalarField& V = mesh_.V(); // Cell volume Foam::Field >& USource = eqn.source(); // get source from fvMatrix //scalarField & Udiag = eqn.diag(); // get diagonal of fvMatrix const scalarField& rho = mesh().lookupObject("rho"); // get density field const vector g (0,-9.81,0) // gravitational vector forAll(cellsZone,i) { const label celli = cellsZone[i]; USource[celli] -= rho[celli]*g/V[celli]; } #}; codeCorrect #{ Pout<< "**codeCorrect**" << endl; #}; codeSetValue #{ Pout<< "**codeSetValue**" << endl; #}; code #{ \$codeInclude \$codeCorrect \$codeAddSup \$codeSetValue #};``` Hope this will be useful in some way. Best regards Mathias wyldckat, raj kumar saini, BlnPhoenix and 6 others like this.

May 16, 2020, 08:45
#11
New Member

Join Date: May 2020
Posts: 11
Rep Power: 6
Quote:
 Originally Posted by msaravia Ok ok, no more comments that do not imply the direct solution to your problem. Sometimes in science, questions are answers... Anyway, what I meant was: A body force is an action exerted on matter essentially by: i) gravity, ii) magnetic fields, iii) electric fields, iv) centrifugal and Coriolis accelerations. From the above, gravity is often taken as constant. As you mention that your force is constant; then I assumed that, whatever this force is, you can add it to gravity. For some solvers, gravity is added at RUNTIME as a vector, that is read from a dictionary. You can check the squareBump example in the OpenFOAM tutorials subdirectory. Thus, your new gravity vector in your input dictionary would look like: Code: `g g [0 1 -2 0 0 0 0] (fx fy -9.81);` or something similar. This approach would avoid programming a new solver. Of course this implies that your solver allows this body force. Finally, I think that Tobi's suggestion is the best approach.

Hey Martin,

i am trying to implement coriolis and centrifugal accelerations to a case without changing the solver. Any suggestion how to implement this ?
I tried tabulatedAccelerationSource, but keep getting errors no matter what i change.

Kind Regards

Daniel

May 16, 2020, 11:32
#12
Member

Martin
Join Date: Dec 2011
Posts: 40
Rep Power: 14
Quote:
 Originally Posted by I7aniel Hey Martin, i am trying to implement coriolis and centrifugal accelerations to a case without changing the solver. Any suggestion how to implement this ? Daniel
Dear Daniel, I don't think you can use the approach I proposed with the Coriolis force, since you need to specify a different force for each cell. Modifying the gravity vector applies the same force to each cell (providing the density is constant).

You should adopt the approaches proposed by Tomi or Bruno. I do not think there is an easy hack for Coriolis.

Good luck!

Martin

 May 18, 2020, 03:00 #13 New Member   Join Date: May 2020 Posts: 11 Rep Power: 6 Hey Martin, thanks for your reply. I tried to implement the rotation via FvOptions using the tabulatedAccelerationSource, unfortunatly this did not work for me, I also tried using MRFProperties an i tried to rotate my Mesh with dynamic mesh. None of that worked so back to start i guess... Regards Daniel

May 18, 2020, 16:23
#14
Member

Martin
Join Date: Dec 2011
Posts: 40
Rep Power: 14
Quote:
 Originally Posted by I7aniel Hey Martin, thanks for your reply. I tried to implement the rotation via FvOptions using the tabulatedAccelerationSource, unfortunatly this did not work for me, I also tried using MRFProperties an i tried to rotate my Mesh with dynamic mesh. None of that worked so back to start i guess... Regards Daniel

I think that writing a new solver (and probably a library) to take into account this force should not be extremely difficult. But of course you need to know some details about OpenFOAM classes.

May 19, 2020, 16:56
#15
Senior Member

Josh McCraney
Join Date: Jun 2018
Posts: 220
Rep Power: 9
Quote:
 Originally Posted by msaravia I think that writing a new solver (and probably a library) to take into account this force should not be extremely difficult. But of course you need to know some details about OpenFOAM classes.
Hi Martins

Can you briefly explain how this is done? I ask because I am trying to add a time-dependent source term cos(omega*t), so I believe I would edit as follows

Code:
```tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(U) + fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
+ mySource
);```
as Tobias suggested, and then define the source term as

Code:
```# define omega 3

const dimensionedVector mySource("mySource", dimensionSet(0,1,-1,0,0,0,0), gunits*Foam::sin(runTime.value()*omega)*vector(0,1,0));```
where I think runTime.value() is equivalent to t. Any help on anything I'm missing is so very much appreciated!

May 19, 2020, 18:21
#16
Member

Martin
Join Date: Dec 2011
Posts: 40
Rep Power: 14
Quote:
 Originally Posted by joshmccraney Hi Martins Code: ```tmp tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevReff(U) == fvOptions(U) + mySource );```

Hi Josh, it seems you are using the correct approach. You can access to the current time with runTime.value(). It is not working?

May 19, 2020, 19:26
#17
Senior Member

Josh McCraney
Join Date: Jun 2018
Posts: 220
Rep Power: 9
Quote:
 Originally Posted by msaravia Hi Josh, it seems you are using the correct approach. You can access to the current time with runTime.value(). It is not working?
Haven't yet ran it, but the dimensionSet is wrong, isn't it (this is what Tobius seemed to imply in his earlier post)?

Last edited by joshmccraney; May 19, 2020 at 20:53.

May 19, 2020, 22:55
#18
Member

Martin
Join Date: Dec 2011
Posts: 40
Rep Power: 14
Quote:
 Originally Posted by joshmccraney Haven't yet ran it, but the dimensionSet is wrong, isn't it (this is what Tobius seemed to imply in his earlier post)?
I think the units should be (0,1,-2,0,0,0,0). Please check it !

May 20, 2020, 01:05
#19
Senior Member

Josh McCraney
Join Date: Jun 2018
Posts: 220
Rep Power: 9
Quote:
 Originally Posted by msaravia I think the units should be (0,1,-2,0,0,0,0). Please check it !
Okay, so when I run wmake for this UEqn.H file

Code:
```    MRF.correctBoundaryVelocity(U);

# define omega 0.05

const dimensionedVector mySource("mySource", dimensionSet(0,1,-2,0,0,0,0), 10*Foam::sin(runTime.value()*omega)*vector(0,1,0));

fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(rho, U)
==
fvOptions(rho, U)
+ mysource
);

UEqn.relax();
.
.
.```
(green is what I added) I get the error
Code:
```In file included from bfInterFoam.C:142:0:
UEqn.H: In function ‘int main(int, char**)’:
UEqn.H:14:9: error: ‘mysource’ was not declared in this scope
+ mysource
^~~~~~~~
UEqn.H:14:9: note: suggested alternative: ‘mySource’
+ mysource
^~~~~~~~
mySource
/opt/openfoam6/wmake/rules/General/transform:25: recipe for target '/opt/openfoam6/platforms/linux64GccDPInt32Opt/applications/solvers/multiphase/bfInterFoam/bfInterFoam.o' failed
make: *** [/opt/openfoam6/platforms/linux64GccDPInt32Opt/applications/solvers/multiphase/bfInterFoam/bfInterFoam.o] Error 1```
Any idea what's wrong, or rather how to declare mysource?

 May 20, 2020, 11:27 #20 Senior Member   Josh McCraney Join Date: Jun 2018 Posts: 220 Rep Power: 9 For anyone in the future trying to add a body force, my UEqn.H file looks like this Code: ``` MRF.correctBoundaryVelocity(U); # define omega 0.05 const dimensionedVector mySource("mySource", dimensionSet(1,-2,-2,0,0,0,0), 1000*Foam::sin(runTime.value()*omega)*vector(0,1,0)); fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + MRF.DDt(rho, U) + turbulence->divDevRhoReff(rho, U) == fvOptions(rho, U) + mySource ); UEqn.relax(); fvOptions.constrain(UEqn); if (pimple.momentumPredictor()) { solve ( UEqn == fvc::reconstruct ( ( mixture.surfaceTensionForce() - ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh) ) * mesh.magSf() ) ); fvOptions.correct(U); }``` Works no issues (at least it ran and the results look good). Green is what I added. Again, this is OF6, a copy of the interFoam solver. Thanks for all the help Martin! I should ask: do I need to do anything to the pEqn.H or is this all that's necessary to add a body force? raj kumar saini, msaravia, ruloz and 1 others like this.