CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   pimpleFoam Body Force (https://www.cfd-online.com/Forums/openfoam-solving/77650-pimplefoam-body-force.html)

mdjames June 29, 2010 16:42

pimpleFoam Body Force
 
Hoping someone has been down this dark and winding trail before. I've seen some modifications to other codes, for instance icoFoam and channelFoam but not this guy...

I'm a Fortran-convert new FOAM user and am having a good amount of trouble, trying to add a source term to the pimpleFoam solver.

I take it it's not as easy as adding a constant to:

tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);

in the UEqn.H file at $FOAM_SOLVERS/incompressible/pimpleFoam

I'd like to have a driven half-channel flow with cyclic b.c's all around. I'm trying to get away from using inlet/outlet patches. I was using fixedValue for pressure and pressureInletVelocity for U but wasn't getting the behavior I need. I'm looking for a little less assumption on the code's part (hence the desire for the body force)

Thanks a TON for any help!

alberto June 30, 2010 03:52

In principle, you can add the forcing term (it is a constant vector, not a scalar constant, since you are formally imposing a constant mean pressure gradient) to the UEqn. There is nothing that prevents you from doing that.

If you want to do it in a more careful way, you might want to include it in the flux instead than in the UEqn, before solving for the pressure equation. The treatment is identical to the one used for the gravity, for example, in bubbleFoam, and it requires some more modifications to the code.

Best,

mdjames June 30, 2010 18:05

Ah, great; she's off and running. Thanks for the quick reply!

I created an H file (class?) to read a dict file which specifies the body force as a vector and its dimensions. After that, I stuck the force into the UEqn.

If you'll pardon my ignorance, what do you mean by including it in the flux?

mdjames June 30, 2010 18:10

Oh, and why is that preferable to including the term in the UEqn?

alberto June 30, 2010 18:23

Including force terms in the flux is a way to improve the robustness of the algorithm in some cases (for example when the force is related to the density, like for the gravity).

If you write your momentum equation as

ddt(U) + div(UU) = div(tau) - grad(p) + F

where F is your forcing term, and then you write in semi-discrete form you get

A*U = H + F

where I put everything I do not want to carry around in H. Now, the predicted velocity is

U = H/A + F/A

and the flux, obtained interpolating this on the cell faces (indicated by _f) is

phi = (H/A)_f . S + (F/A)_f . S

In your case F is constant, so

phi = (H/A)_f . S + F(1/A)_f . S

In OpenFOAM notation, this means that you should add

surfaceScalarField rUAf = fvc::interpolate(rUAf)

surfaceScalarField phiForcing = rUAf*(F & mesh.Sf())

to your flux, after adjustPhi. At this point however you have to modify the correction step for the velocity too (the flux is fine), consistently with the new velocity predictor:

U += fvc::reconstruct(phiForcing) - rUA*fvc::grad(p);

On a side note, you probably won't need this in your case, but it is better to know ;)

Hopefully I did not introduce errors while typing :)

Best,

alberto June 30, 2010 19:12

Quote:

Originally Posted by mdjames (Post 265205)
Oh, and why is that preferable to including the term in the UEqn?

Numerically this is done to improve the robustness of the algorithm, especially with strongly varying terms (gravity affects very differently a denser zone with respect to a less dense one).
In other words, it is a way to improve stability when enforcing the force balance.

Best,

mdjames July 2, 2010 22:06

Oh, okay. I think the first post makes sense to me, if i understand the new variables. The second one definitely does.

Thanks again...you're right, I would like to know the correct way to do it rather than a way.

Next up is getting my forcing term time-dependent..."here be dragons"

I'll check back if I can't manage to work something out.

alberto July 3, 2010 01:39

Quote:

Originally Posted by mdjames (Post 265550)
Oh, okay. I think the first post makes sense to me, if i understand the new variables. The second one definitely does.

Ask if something is not clear :-)

Quote:

Thanks again...you're right, I would like to know the correct way to do it rather than a way.
Let's say nothing is written in stone. If you work with incompressible single-phase codes, you probably won't see any difference with the two approaches. Things are very different in compressible or multiphase flows.

Quote:

Next up is getting my forcing term time-dependent..."here be dragons"

I'll check back if I can't manage to work something out.
The principle is the same, but of course ask if you have problems.

Best,

mdjames July 8, 2010 20:59

I think I'm hot on the trail, I just can't find to seem to find a way to call out the runtime as a dimensioned scalar. I tried creating a new time variable with:

dimensionedScalar timemdj
(
timemdj [0 0 1 0 0 0 0 ] runTime.Value()
);

to include in my equation for the body force but to no avail...

This is getting fun! :)

alberto July 9, 2010 02:40

It is

dimensionedScalar timemdj
(
"timemdj",
dimensionSet(0, 0, 1, 0, 0),
<yourExpression> // Replace with what you need
);

or, simply

dimensionedScalar timemdj = <yourExpression>;

Best,

mdjames September 9, 2010 12:03

I'm confused...

my forcing is described as

dimensionedVector waveP=amp*cos(-freq*time_); //determining the new forcing

and this yields a forcing in the x-direction as I wanted but I realized:

1. I never included any information to make this an x forcing
2. I'm adding a vector to a matrix (in the UEqn, below)???

tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
- waveP
);

Is there some assumption on OF's part or something?

alberto September 9, 2010 13:09

How do you define Amp? It has to be a vector for the code to compile...

Best,

mdjames September 9, 2010 13:15

Okay, I see now...minor freak-out.

I don't quite understand how the vector waveP is added to higher dimension arrays like U...

alberto September 9, 2010 13:22

Quote:

Originally Posted by mdjames (Post 274651)
Okay, I see now...minor freak-out.

I don't quite understand how the vector waveP is added to higher dimension arrays like U...

waveP is a vector of three components, since amp has to be a vector of three components, multiplied by the time-dependent part.

Each component of waveP is added to the corresponding scalar equation for the component of U. So, say you have

g= (gx, gy, gz)

and you do

UEqn += g;

You add gx to the equation for Ux, gy to the equation for Uy and so on. OpenFOAM syntax just makes it easier for you to do this in only one step.

Best,

mdjames September 9, 2010 15:33

Okay, to check my understanding of the grand scheme of things:

Each node has a vector equation associated with it. So, in effect we have the matrix of equations (UEqn) added to a matrix of the same size, with all elements being the vector waveP

correct?

If so, I suppose creating a spatially variable forcing would just involve directly specifying a forcing matrix.

Am I on the right track here?

Thanks a ton!

alberto September 9, 2010 16:11

Yes. Put into the code, you simply define a volVectorField, initialize it with the values of your force term, and add it to UEqn.

P.S. You amp is (|amp|, 0, 0) right? :)

mdjames September 9, 2010 16:19

Beautiful! Seems like this may be a fairly painless solution. Much simpler than some other schemes I had rolling around in my head.

Yep, just without the commas :) It's read from a wavePressDict like viscosity is.

Thomas Baumann October 7, 2010 11:27

Hi,

I'm trying to implement new anisotropic turbulence models in OpenFOAM.
It's like a normal isotropic turbulence-model, with an explicit nonlinear part too.

So the Reynolds-fluxes are defined as: R = linearStress + nonlinearStress (explicit)

so divDevReff=
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))
+ fvc::div(nonlinearStress_)

Would it make here sense to put the nonlinear-part in the flux-part as descirbed above? Would the results be better than putting all of them in a fvVectorMatrix? Or is this only for the convergence, the robustness of the system better?

There are some nonlinear models in OF-1.7 like the nonlinearshih model, but here the nonlinearStress_-term is in the fvVectorMatrix, too. Using them to solve a fully developed channel flow I have velocities out of the wall. Could it be the reason, because the nonlinearStress-part is described in the fvVectorMatrix?

Thanks and regards,
Thomas

be_inspired July 3, 2014 08:16

Hi Alberto,

In case the force is managed in the first way ( force term in the flux), my question is if Ueqn must be modified or not.

I am working over the rotorDiskSource (OF2.3.0) and some pressure wiggles happen when discreted body forces are included directly over the Ueqn.
I think that your approach could improve this situation.
http://www.openfoam.org/mantisbt/view.php?id=1313

A*U = H + F

I think it should be:
A*U = H - grad[p] + F
U=H/A - 1/A*grad[p] + F/A

Excuse me for using the same topic but there are other topics where others have reported this bug and no answer have been given.

http://www.cfd-online.com/Forums/ope...rce-model.html
http://www.cfd-online.com/Forums/ope...l-wiggles.html

I am trying to implement the body force in this other way and I am not able to do it correctly. Could you please help me?

manoj jaiswal July 10, 2014 08:04

urgent help for body force
 
Hi i am a new Foamer and i am simulating flow over a cylinder(in pimplefoam) . i want to add a body force in the flow , so what are the steps needed for the same . I have gone through a tutorial of how to add temperature in icoFoam (http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam) but still not clear about the addition of body force.

Thanks :)





Quote:

Originally Posted by alberto (Post 265207)
Including force terms in the flux is a way to improve the robustness of the algorithm in some cases (for example when the force is related to the density, like for the gravity).

If you write your momentum equation as

ddt(U) + div(UU) = div(tau) - grad(p) + F

where F is your forcing term, and then you write in semi-discrete form you get

A*U = H + F

where I put everything I do not want to carry around in H. Now, the predicted velocity is

U = H/A + F/A

and the flux, obtained interpolating this on the cell faces (indicated by _f) is

phi = (H/A)_f . S + (F/A)_f . S

In your case F is constant, so

phi = (H/A)_f . S + F(1/A)_f . S

In OpenFOAM notation, this means that you should add

surfaceScalarField rUAf = fvc::interpolate(rUAf)

surfaceScalarField phiForcing = rUAf*(F & mesh.Sf())

to your flux, after adjustPhi. At this point however you have to modify the correction step for the velocity too (the flux is fine), consistently with the new velocity predictor:

U += fvc::reconstruct(phiForcing) - rUA*fvc::grad(p);

On a side note, you probably won't need this in your case, but it is better to know ;)

Hopefully I did not introduce errors while typing :)

Best,


jinheng October 24, 2016 00:50

How should I define a acceleration-terms.dat in FvOptions-tabulatedAccelerationSource
 
hi everyone, I want to define a time-change acceleration in sloshing movement. And maybe Fvoptions is a good way. Similar with tanksloshing3D6dof case in interDyMFoam. I want to define a acceleration term in non-inertial reference frame in interFoam. But I can not read the acceleration-term.dat using Fvoptions due to the formate of file.
So my question is :
1、Is there anywhere tell the formate of tabulatedaccelerationsource.
2、otherthan fvoptions, is there any good method to achieve non-inertial reference frame in interFoam. (I have try to added a "acc" just in the Uequ ,but the alpha.water does not change with "acc")

thanks very much for any advises or help.
Usage
Example usage:
\verbatim
SBM
{
type tabulatedAccelerationSource;
active yes;

tabulatedAccelerationSourceCoeffs
{
timeDataFileName "constant/acceleration-terms.dat";
}
}
\endverbatim

heng Jin
OUC

jinheng October 24, 2016 04:15

Quote:

Originally Posted by jinheng (Post 622657)
hi everyone, I want to define a time-change acceleration in sloshing movement. And maybe Fvoptions is a good way. Similar with tanksloshing3D6dof case in interDyMFoam. I want to define a acceleration term in non-inertial reference frame in interFoam. But I can not read the acceleration-term.dat using Fvoptions due to the formate of file.
So my question is :
1、Is there anywhere tell the formate of tabulatedaccelerationsource.
2、otherthan fvoptions, is there any good method to achieve non-inertial reference frame in interFoam. (I have try to added a "acc" just in the Uequ ,but the alpha.water does not change with "acc")

thanks very much for any advises or help.
Usage
Example usage:
\verbatim
SBM
{
type tabulatedAccelerationSource;
active yes;

tabulatedAccelerationSourceCoeffs
{
timeDataFileName "constant/acceleration-terms.dat";
}
}
\endverbatim

heng Jin
OUC

And I have attempted this thread in CFD-ONLINE(http://www.cfd-online.com/Forums/ope...blefile-h.html). I use the same format of the tablefile. But it still give me the same errors:
1、Expected a '(' while reading Tuple2, found on line 2 the label 0
2、wrong token type - expected Scalar, found on line 2 the punctuation token '('

peyman.davvalo.khongar May 16, 2018 08:40

Moi!

Check this:


https://github.com/OpenFOAM/OpenFOAM...5900ea5eb546af

Moikka!

Peyman

Msure June 17, 2021 06:07

Quote:

Originally Posted by be_inspired (Post 499834)
Hi Alberto,

In case the force is managed in the first way ( force term in the flux), my question is if Ueqn must be modified or not.

I am working over the rotorDiskSource (OF2.3.0) and some pressure wiggles happen when discreted body forces are included directly over the Ueqn.
I think that your approach could improve this situation.
http://www.openfoam.org/mantisbt/view.php?id=1313

A*U = H + F

I think it should be:
A*U = H - grad[p] + F
U=H/A - 1/A*grad[p] + F/A

Excuse me for using the same topic but there are other topics where others have reported this bug and no answer have been given.

http://www.cfd-online.com/Forums/ope...rce-model.html
http://www.cfd-online.com/Forums/ope...l-wiggles.html

I am trying to implement the body force in this other way and I am not able to do it correctly. Could you please help me?

Hi,

Have you solved the problem? I am adding a direct force immersed boundary method into pimpleFoam, if I add the extra force from solid in the momentum predictor step, it can run but the results are wrong. I think I should add the force in pressure iterative step, because the force is chaging in every step. Thanks.

mostanad November 4, 2021 21:59

Quote:

Originally Posted by Msure (Post 806278)
Hi,

Have you solved the problem? I am adding a direct force immersed boundary method into pimpleFoam, if I add the extra force from solid in the momentum predictor step, it can run but the results are wrong. I think I should add the force in pressure iterative step, because the force is chaging in every step. Thanks.


Hi Shuo,
Have you solved the problem? I am doing the same thing for IBM, but the results are so dependent on the Number of Piso loops. I think something is wrong with the way I am imposing my force in pressure piso loop. So can you please share to me your experience?
Cheers,
Mohammad


All times are GMT -4. The time now is 03:45.