Variation of gravity with time
 

Dear OpenFoam users,

i'm trying to simulate the effects of the acceleration's variation of a rocket on a particular experiment. My problem is that the acceleration value varies with the mission time.

Does anyone know how to set the variation of the gravity with the time?

Thank you
Bye

Whyman thanks for posting this thread ,, I'm too have a simulation that need variation of gravity with time ,, also i want to know how to implement centrifugal force in openFoam


regards

Had you guys a look at the sloshing tank tutorial:

tutorials/multiphase/interDyMFoam/ras/sloshingTank2D

This will properly match you expectation.

kinda regards

hi idrama

thank you ,, it is a rotating wall boundary condtion
do you know how to implement a given body force ( and also how to make it vary with time ) in to the solver governing equations ' interFoam ' for example

Hi,

look at how the gravity is implemented in interFoam, by means of its contribution to the flux.

In the momentum predictor, if solved, the gravity is not included in the matrix but added on the RHS of the equation as (UEqn.H)

fvc::reconstruct(fvc::interpolate(rho)*(g & mesh.Sf()))

and its contribution to the flux is accounted for as (pEqn.H)

phi = phiU + (fvc::interpolate(rho)*(g & mesh.Sf()))*rUAf;

You can do the same for other body force terms, starting from the semi-discrete form of the momentum equation and deriving the expression for the flux contribution of the additional term.

If you simply want the gravity to change with time, simply remove

• Remove #include "readGravitationalAcceleration.H"
• Define a volVectorField that contains your gravitational acceleration, depending on the current time
  • volVectorField g = <put your expression here>;

Best,

Alberto,

thanks for your useful reply.
I just have only one more question about this topic.

You suggested to write a particular expression for the gravity. What am i supposed to do in case i have a list of gravity values for each flowtime, which do not respond to a specific expression?

For example, if i have this type of gravity data file:

flowtime gx gy gz
0.1 0.1 0.1 1
0.15 0.2 0.1 1.5
0.2 0.8 -0.5 1.5
.........


Finally, is it possible to apply this variation to OpenFoam version 1.5?

Thank you
Regards

Hello Alberto,

in my case the gravity only turns about the y-axis in 4 sec. So I changed the code according to your last post as follows:

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
//#include "readGravitationalAcceleration.H"
volVectorField g=-9.81*vector(0,0,1); // There may be dependencies in the next includes.
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "readTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"

runTime++;

g=9.81*sin(pi/4*runTime)*vector(1,0,0)-9.81*cos(pi/4*runTime)*vector(0,0,1);

Info<< "Time = " << runTime.timeName() << nl << endl;


But when I am compiling it the compiler ends up with an error message. What did I wrong?

kinda regards

Claus

Well, it might not be the most elegant way, but a quick one could be to write a function that returns the gravity value checking the time range.

In pseudo-code

if (t <= 0.1)
g = (0.1 0.1 1)
else if (0.1 < t < 0.15)
g = (0.1 0.2 1.5)
else if (...)
...

[QUOTE
Finally, is it possible to apply this variation to OpenFoam version 1.5?

Yes. In OF 1.6 the gravity field has been redefined in a slightly different way than in 1.5, but this should not make a big difference.

Best,

Without the error message, it is a bit hard :-)

However, here

g=9.81*sin(pi/4*runTime)*vector(1,0,0)-9.81*cos(pi/4*runTime)*vector(0,0,1);

you use runTime in the wrong way. It is a Time object, but what you need is your current time, which is provided by runTime.timeName().

In addition, you do not provide the correct units to g, so the solver will complain when running. You can define a multiplier with the correct units as

const dimensionedScalar gunits("gunits", dimensionSet(0,1,-2,0,0,0,0), 1.0);

and multiply your value by it.

Best,

Hello Alberto,

for one week I am trying it to get the code going. I modified the code as follows:

vector g0(0,0,9.81);

#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
//#include "readGravitationalAcceleration.H"
dimensionedVector g("g",dimensionSet(0,1,-2,0,0,0,0),g0);
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "readTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"

runTime++;

g0=9.81*sin(runTime.timeName()*pi/4.0)*vector(1,0,0)-9.81*cos(runTime.timeName()t*pi/4.0)*vector(0,0,1);

dimensionedVector g("g",dimensionSet(0,1,-2,0,0,0,0),g0);

Info<< "Time = " << runTime.timeName() << nl << endl;

By compinling I get the error message:

MyInterFoam.C:87: error: no match for ‘operator*’ in ‘Foam::Time::timeName() const() * 3.141592653589793115997963468544185161590576171875 e+0’
/home/idrama/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/dimensionSet.H:266: note: candidates are: Foam::dimensionSet Foam::operator*(const Foam::dimensionSet&, const Foam::dimensionSet&)

As far as I figured out the problem is that I have double operation * but the data type of tunTime.timeName() is word. How can I get a double value from runTime? Or, in general, what I have to to else to get the code running?

Cheers,

Claus

The problem is in how you set the dimensions. Replace
with

and remove the definition of vector g0 at the beginning.

In other words you define a constant, of magnitude 1 with the dimensional units of an acceleration, and you use it as a multiplier for your vector, to set the appropriate unit.

Best,

Thanks All!

PROBLEM SOLVED, MISSION ACCOMPLISHED:

#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"

#define pi 3.141592653589793238


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{

const dimensionedScalar gunits("gunits", dimensionSet(0,1,-2,0,0,0,0), 9.81);


#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "readTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"

runTime++;

g=gunits*Foam::sin(runTime.value()*pi/4.0)*vector(1,0,0)-gunits*Foam::cos(runTime.value()*pi/4.0)*vector(0,0,1);


Info<< "Time = " << runTime.timeName() << nl << endl;

Kinda regards!

Claus

Good! And you found part of the answer yourself, since I misled you with timeName(). Sorry about that :(

Best,

Hi alberto ;

i want to make the gravity vector constant in magnitude but rotates with a known omega ,, and in the same situation there exist a constant centrifugal force (constant magnitude and constant direction).. both in one simulation

gravity vector rotates (omega) + centrifugal force (constant vector)


does it possible to do that in the runtime as in the previous case ?

best regards.

Hi alberto ;

i want to make the gravity vector constant in magnitude but rotates with a known omega ,, and in the same situation there exist a constant centrifugal force (constant magnitude and constant direction).. both in one simulation

gravity vector rotates (omega) + centrifugal force (constant vector)


does it possible to do that in the runtime as in the previous case ?

best regards.

further problems
 

Hi guys,

i would like to replace my problem, beacause i can't solve it yet.

Alberto wrote about a pseudo-code for the gravity in case i had random values.

My problems is: what happens if i have a list of thousands values of g?

My idea is the following.

1) To write a gravity file as a text file with all the x,y,z values

2) To write the time file with the list of the time values

3) Every time-step, the solver will read these two files (g and t) and change the g-value in the "g" file of the case contained in the "constant" directory

4) The solver keep calculating and cycle it.


Alternatively i received this comment:

"
Obviously the standard OpenFOAM solvers don't include the variation in
gravity - most of them apply to flows in a stationary reference frame
on Earth!

If you want to use data from a file whose values you wish to
interpolate between, there is a class set up for this called
interpolationTable<Type>. You could use this for a vector quantity,
by constructing an interpolationTable<vector>. There are a choice of
constructors that take a file name, or dictionary entry for the file
name as arguments, see:
$FOAM_SRC/OpenFOAM/interpolations/interpolationTable/interpolationTable.C

Some time varying boundary conditions use the same class to handle
time-value data. Examples are:
$FOAM_SRC/finiteVolume/fields/fvPatchFields/derived/timeVaryingFlowRateInletVelocity/
$FOAM_SRC/finiteVolume/fields/fvPatchFields/derived/timeVaryingUniformFixedValue/
$FOAM_SRC/finiteVolume/fields/fvPatchFields/derived/timeVaryingUniformTotalPressure/

The gravity vector can be overridden with a vector from your
interpolationTable<vector> in whichever solver you need to use. If
the interpoltationTable<vector> you construct is called myInTable, the
cose would look something like:

g.value() = myIntTable(runTime.value());
"
Unfortunately, i don't know how to implement these ideas to the code.

What do you think about it? Do you have any suggestion?

Regards

Stefano

Thousands of values as a function of time, or does your "gravity vector" change as a function of the position too?

Best,
Alberto

Thousands of values as a function of time. It is uniform in space (i mean constant in position).

Imagine the small example that i wrote before but multiplied by thounsands values.


flowtime gx gy gz
0.1 0.1 0.1 1
0.15 0.2 0.1 1.5
0.2 0.8 -0.5 1.5
.........(thousands rows)....


Stefano

Hi,

I agree with the suggestion of using an interpolationTable. You can take a look at this:

http://foam.sourceforge.net/doc/Doxy....html#_details

and use the constructor from the name of the file containing the data table, for example, or, if you prefer, you can set the file name containing the data in a specified dictionary and construct using that dictionary (as done in timeVaryingUniformFixedValue).

What problems did you meet when trying to implement this?

Best,

My only problem is that i'm not so expert yet to change the code. I'm still studying it, but it's still quite complicate for me.

Moreover the interpolation file is used (i think) only to interpolate between two values: but how does the code read the file, select the right value for each flowtime and change the value in its calculation?



Regards

Stefano

OpenFoam 1.7.x, Gravity rotation for interFoam
 

Dear all
For OF 1.7.x, my "interFoam.C" looks like this with comments for my self (had to change it sligtly relative to the above discussion) Hope this is of help to someone! :)
NOTE: // a: //b: etc. marks the changes.

---------------------------------------------------------
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"
#include "interpolationTable.H"

// a:
#define pi 3.141592653589793238

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{

// b:
const dimensionedScalar gunits("gunits", dimensionSet(0,1,-2,0,0,0,0), 9.81);

#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"


#include "readTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"

runTime++;

Info<< "Time = " << runTime.timeName() << nl << endl;

// ------------------
// kemur fra creatFields.H, virdist breyta litlu ad hafa thetta a:
// dimensionedVector g0(g);
// ------------------

// c:
// the file ./constants/g seems to be overwritten or ignored.
// Have tested by setting gravity to -0.05 m/s2 as well as -99.8 m/s2 in
// ./constants/g => same result is produced regardless! I.e. the line
// below seems to dominate over anything that is written in ./constants/g.
g=gunits*Foam::sin(runTime.value()*pi/2.0)*vector(1,0,0)-gunits*Foam::cos(runTime.value()*pi/2.0)*vector(0,0,1);

// ------------------
// Comes from creatFields.H, virdist breyta litlu ad hafa thetta a:
// dimensionedVector g0(g);
// ------------------

// --------- from createFields.H
// d:
// Comes from the file createFields.H. You dont have to delete the lines
// in createFields.H. The field gh is just initilaized there.
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
// -------------------------------------


twoPhaseProperties.correct();

#include "alphaEqnSubCycle.H"

#include "UEqn.H"

// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}

turbulence->correct();

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
}

This thread is a bit old but hopefully someone could help....

i have the same problem, wanna change g and make it as a function of time. I have OF 2.1 and in interFoam solver file i cant find the file : readGravitationalAcceleration.C ... I tried to change the interFoam.C file like diescribed above but it didnt work (i ran interfoam as normal on a test case and it just takes the constant value of g in "constant" file)... i tried to write the codes in creatFields.H but it didnt work either (again interFoam runs normal with no changes as if i didnt change anything in the code).

I would really appreciate some help or advice

Thanks

interFoam 2.1.x Rotation
 

Hi Callahance

I just compiled gravity rotation into interFoam 2.1.x. The files needed to be changed/modified are the interFoam.C and the UEqn.H

here are my modifications (search for VVPF in interFoam.C and UEqn.H) to see the changes.

Note that this compile without any problems, but I haven't had time to test the binaries. Hope this is enough to get you started:)

cheers
JonW

JonW... all what i could say is : THANKS... ill try it out and give a feedback... thanks very much again

Hey Claus,

Hello. This is a very old post. Its feels quite strange to ask a question now. Please bear with me. I want to know how the code worked since you never defined the variable 'g'. I'm trying to do something similar for a plunging airfoil. I get an error " In function ‘int main(int, char**)’: , ‘g’ was not declared in this scope "

If you no longer remember thats OK.

Thanks,
Pru.

you have to define g in createFields.H with the #include "readGravitationalAccelation.H" - thingi (see the createFields.H in the interFoam solver).

If you are going to use gravity in a single phase fluid, check out
http://www.cfd-online.com/Forums/ope...interfoam.html

J.

Heaving reference frame
 

Hey john thanks for the reply!

I'm simulating a flapping airfoil and I want to use a heaving reference frame instead of using a deforming mesh. This means that the fluid should accelerate up and down in a sinusoidal motion. So I changed the Ueqn as follows

I defined max_ acceleration and frequency in createFields.H as follows

When I do wmake I get the following error.
However Claus managed to run his code.

I want to know where I went wrong in the process. I will check out the link you mentioned.

i made it whyman........ simply enter following code in your creatFields.H


// Read the data file and initialise the interpolation table
interpolationTable<vector> timeSeriesAcceleration
(
runTime.path()/runTime.caseConstant()/"acceleration.dat"
);
g.value() = timeSeriesAcceleration(runTime.value());

compile it

and then make an acceleration.dat file in your constants case directory in following format

4
(
(0 (1000 0 0))
(3 (1000 0 0))
(4 (1000 0 0))
(5 (12000 0 0))
)

run the case and enjoy

I tried to specify acceleration from file acceleration.dat file. but there is a problem. solver gets only first value of acceleration and after that it dont read acceleration.dat file. i added following code in interfoam creatFields.h

// Read the data file and initialise the interpolation table
interpolationTable<vector> timeSeriesAcceleration
(
runTime.path()/runTime.caseSystem()/"acceleration.dat"
);
g.value() = timeSeriesAcceleration(runTime.value());

please try and help me too......

createFields.H is just read at the beginning of the simulation. You should update g.value() at each time step

ya i know problem is creatfields.h is read once but......
 

yes boss i know creatFields.H is read just once. and i have to update g.value at each time step.... but the problem is how to do it..... i have tried by adding while loop in creatFields for g.value() but it doesnt work. as you said creatFields is read only once. if i take out g.value form creatFields.H and place it in interFoam g.value does not replace g vector specified in constants directory of the case. please suggest if you have an idea.....

While loops in createFields.H doesn't make sense, since the timeStep isn't updated. Update g.value() in runTime while loop (interFoam.C).

For example:

sir i tried this but it doesnt solve for acceleration from acceleration.dat instead it solves from g file. i have already tried this.....

i have tried it by keeping g.value() = timeSeriesAcceleration(runTime.value()); at the same place as you kept... but doing this solver stops taking value from acceleration.dat it takes g value from g file. as it solves g in creatFields for gh and p_gh both gets g value.....

any other suggestion sir???

then you should update gh and ghf, too.

thank you sir it works

Hi，clause. Your answer to changable acceleration seems correct, but how can I handle the g file in const folder. Thanks in advance.

no match for operator error in implementing variable gravity
 

Based on the direction given in this thread,

I tried using the line
g=gunits*vector(0,Foam::sin(runTime.value()*pi*2.0 *freq),1) for sinusoidal lateral acceleration


I get an error

interFoam2.C:114:3: error: no match for âoperator=â (operand types are âconst Foam::meshObjects::gravityâ and âFoam::dimensioned<Foam::Vector<double> >â)
g=gunits*vector(0,Foam::sin(runTime.value()*pi*2.0 *freq),1);

I also tried the line given in the thread as such to test. That too gave a similar error. Am I missing some step? I do hope I get a response considering that this is an very old thread.

I am using v1812

I dont think you can put the "Foam::"-thingi into vector()
Rather try this,...
g=gunits*Foam::sin(runTime.value()*pi*2.0*freq)*ve ctor(0,1,0) + gunits*vector(0,0,1);


If this is indeed gravity, make sure that abs(g) = 9.81. Like this vector is now, then its not the case.
Hope this helps.