CFD Online Discussion Forums
Page 1 of 2 1 2
Show 40 post(s) from this thread on one page

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   How to add a body force to a momentum equation (https://www.cfd-online.com/Forums/openfoam-programming-development/209691-how-add-body-force-momentum-equation.html)

anon_q October 22, 2018 19:20
How to add a body force to a momentum equation
 
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. :mad:

2) if it must involve programming, then where in my code (I mean in the pimple loop) should I put the body force?

msaravia October 22, 2018 19:58
Cant you include this force by using a modified gravity vector?

anon_q October 23, 2018 08:11
Quote:
Originally Posted by msaravia (Post 711983)
Cant you include this force by using a modified gravity vector?
Please please, if you have no answer do not comment because I don't know what you are talking about.

Tobi October 23, 2018 08:25
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)
);
As you can see, the fvOptions is included as source term. However, adding a new source is as follows:
Code:
tmp<fvVectorMatrix> 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<fvVectorMatrix> 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

anon_q October 23, 2018 09:49
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)
);
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?

msaravia October 23, 2018 09:58
Quote:
Originally Posted by Evren Linda (Post 712172)
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.

Tobi October 23, 2018 10:52
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 :)

anon_q October 23, 2018 11:36
Quote:
Originally Posted by msaravia (Post 712193)
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.

wyldckat November 3, 2018 16:26
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 (Post 711982)
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. :mad:
... 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 (Post 711982)
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

SvenBent November 29, 2018 09:03
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
    #};
Hope this will be useful in some way.

Best regards
Mathias

I7aniel May 16, 2020 08:45
Quote:
Originally Posted by msaravia (Post 712193)
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

msaravia May 16, 2020 11:32
Quote:
Originally Posted by I7aniel (Post 770771)
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

I7aniel May 18, 2020 03:00
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

msaravia May 18, 2020 16:23
Quote:
Originally Posted by I7aniel (Post 770907)
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.

joshmccraney May 19, 2020 16:56
Quote:
Originally Posted by msaravia (Post 771013)
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!

msaravia May 19, 2020 18:21
Quote:
Originally Posted by joshmccraney (Post 771178)
Hi Martins
Code:
tmp<fvVectorMatrix> 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?

joshmccraney May 19, 2020 19:26
Quote:
Originally Posted by msaravia (Post 771187)
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)?

msaravia May 19, 2020 22:55
Quote:
Originally Posted by joshmccraney (Post 771191)
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 !

joshmccraney May 20, 2020 01:05
Quote:
Originally Posted by msaravia (Post 771203)
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?

joshmccraney May 20, 2020 11:27
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?


All times are GMT -4. The time now is 02:51.
Page 1 of 2 1 2
Show 40 post(s) from this thread on one page