
[Sponsors] 
How to add a body force to a momentum equation 

LinkBack  Thread Tools  Search this Thread  Display Modes 
October 22, 2018, 20:20 
How to add a body force to a momentum equation

#1 
Senior Member
Join Date: Mar 2018
Posts: 115
Rep Power: 5 
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? 

October 22, 2018, 20:58 

#2 
Member
Martin
Join Date: Dec 2011
Posts: 38
Rep Power: 12 
Cant you include this force by using a modified gravity vector?


October 23, 2018, 09:11 

#3 
Senior Member
Join Date: Mar 2018
Posts: 115
Rep Power: 5 

October 23, 2018, 09:25 

#4 
Super Moderator

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 <Type>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<fvVectorMatrix> tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence>divDevReff(U) == fvOptions(U) ); Code:
tmp<fvVectorMatrix> tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence>divDevReff(U) == fvOptions(U) + mySource ); Code:
const dimensionedVector mySource("mySource", dimensionSet(0,1,1,0,0,0,0), vector(4,0,0)); tmp<fvVectorMatrix> tUEqn . . . To b) please follow this link and scroll down to the Detailed Description > Usage https://cpp.openfoam.org/v6/classFoa...21f7bbb656d86d Good luck
__________________
Keep foaming, Tobias Holzmann 

October 23, 2018, 10:49 

#5 
Senior Member
Join Date: Mar 2018
Posts: 115
Rep Power: 5 
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<fvVectorMatrix> tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence>divDevReff(U) == fvOptions(U) + mySource ); Code:
tmp<fvVectorMatrix> tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence>divDevReff(U)  mySource == fvOptions(U) ); 2)When using fvOptions, I guess there is something called "vectorCodedSource", right? 

October 23, 2018, 10:58 

#6  
Member
Martin
Join Date: Dec 2011
Posts: 38
Rep Power: 12 
Quote:
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); 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, 11:52 

#7 
Super Moderator

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
__________________
Keep foaming, Tobias Holzmann 

October 23, 2018, 12:36 

#8  
Senior Member
Join Date: Mar 2018
Posts: 115
Rep Power: 5 
Quote:
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, 17:26 

#9  
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,962
Blog Entries: 45
Rep Power: 125 
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 doublecheck 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:
Quote:
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:
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:
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, 10:03 

#10 
New Member
Mathias Poulsen
Join Date: Feb 2018
Location: Denmark
Posts: 9
Rep Power: 5 
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<Foam::Vector<double> >& USource = eqn.source(); // get source from fvMatrix //scalarField & Udiag = eqn.diag(); // get diagonal of fvMatrix const scalarField& rho = mesh().lookupObject<scalarField>("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 #}; Best regards Mathias 

May 16, 2020, 09:45 

#11  
New Member
Join Date: May 2020
Posts: 11
Rep Power: 3 
Quote:
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, 12:32 

#12  
Member
Martin
Join Date: Dec 2011
Posts: 38
Rep Power: 12 
Quote:
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, 04:00 

#13 
New Member
Join Date: May 2020
Posts: 11
Rep Power: 3 
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, 17:23 

#14  
Member
Martin
Join Date: Dec 2011
Posts: 38
Rep Power: 12 
Quote:
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, 17:56 

#15  
Senior Member
Josh McCraney
Join Date: Jun 2018
Posts: 200
Rep Power: 6 
Quote:
Can you briefly explain how this is done? I ask because I am trying to add a timedependent 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 ); 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)); 

May 19, 2020, 19:21 

#16  
Member
Martin
Join Date: Dec 2011
Posts: 38
Rep Power: 12 
Quote:
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, 20:26 

#17 
Senior Member
Josh McCraney
Join Date: Jun 2018
Posts: 200
Rep Power: 6 
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 21:53. 

May 19, 2020, 23:55 

#18 
Member
Martin
Join Date: Dec 2011
Posts: 38
Rep Power: 12 

May 20, 2020, 02:05 

#19  
Senior Member
Josh McCraney
Join Date: Jun 2018
Posts: 200
Rep Power: 6 
Quote:
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(); . . . 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 

May 20, 2020, 12:27 

#20 
Senior Member
Josh McCraney
Join Date: Jun 2018
Posts: 200
Rep Power: 6 
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); } 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? 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Add extra body force (e.g. gravity) in simpleFoam: cavity example  juchess  OpenFOAM Running, Solving & CFD  6  December 17, 2016 08:19 
[PyFoam] and paraview  eelcovv  OpenFOAM Community Contributions  28  May 30, 2016 10:23 
how to add a force in the momentum equation  guillaumem  OpenFOAM  0  June 14, 2010 04:49 
Viscosity and the Energy Equation  Rich  Main CFD Forum  0  December 16, 2009 15:01 
UDF Addidtional Force to momentum equation  sfawal  FLUENT  10  July 15, 2009 10:29 