CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

pimpleFoam Body Force

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

Like Tree11Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 29, 2010, 16:42
Default pimpleFoam Body Force
  #1
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
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!
randolph likes this.
mdjames is offline   Reply With Quote

Old   June 30, 2010, 03:52
Default
  #2
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
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,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   June 30, 2010, 18:05
Default
  #3
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
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 is offline   Reply With Quote

Old   June 30, 2010, 18:10
Default
  #4
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
Oh, and why is that preferable to including the term in the UEqn?
mdjames is offline   Reply With Quote

Old   June 30, 2010, 18:23
Default
  #5
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
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 Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   June 30, 2010, 19:12
Default
  #6
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by mdjames View Post
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,
meth likes this.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   July 2, 2010, 22:06
Default
  #7
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
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.
mdjames is offline   Reply With Quote

Old   July 3, 2010, 01:39
Default
  #8
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by mdjames View Post
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,
XJ_Wang likes this.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   July 8, 2010, 20:59
Default
  #9
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
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!
mdjames is offline   Reply With Quote

Old   July 9, 2010, 02:40
Default
  #10
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
It is

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

or, simply

dimensionedScalar timemdj = <yourExpression>;

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 9, 2010, 12:03
Default
  #11
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
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?
mdjames is offline   Reply With Quote

Old   September 9, 2010, 13:09
Default
  #12
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
How do you define Amp? It has to be a vector for the code to compile...

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 9, 2010, 13:15
Default
  #13
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
Okay, I see now...minor freak-out.

I don't quite understand how the vector waveP is added to higher dimension arrays like U...
mdjames is offline   Reply With Quote

Old   September 9, 2010, 13:22
Default
  #14
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by mdjames View Post
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,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 9, 2010, 15:33
Default
  #15
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
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!
mdjames is offline   Reply With Quote

Old   September 9, 2010, 16:11
Default
  #16
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
alberto will become famous soon enoughalberto will become famous soon enough
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?
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.
alberto is offline   Reply With Quote

Old   September 9, 2010, 16:19
Default
  #17
New Member
 
Matt James
Join Date: Jun 2010
Location: Marinette,WI, USA
Posts: 25
Rep Power: 15
mdjames is on a distinguished road
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.
mdjames is offline   Reply With Quote

Old   October 7, 2010, 11:27
Default
  #18
Senior Member
 
Join Date: Apr 2009
Location: Karlsruhe, Germany
Posts: 103
Rep Power: 16
Thomas Baumann is on a distinguished road
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
Thomas Baumann is offline   Reply With Quote

Old   July 3, 2014, 08:16
Default
  #19
Senior Member
 
M. Montero
Join Date: Mar 2009
Location: Madrid
Posts: 153
Rep Power: 17
be_inspired is on a distinguished road
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?

Last edited by be_inspired; July 3, 2014 at 11:20.
be_inspired is offline   Reply With Quote

Old   July 10, 2014, 08:04
Smile urgent help for body force
  #20
New Member
 
Manoj Jaiswal
Join Date: Jun 2014
Posts: 3
Rep Power: 11
manoj jaiswal is on a distinguished road
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 View Post
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,
manoj jaiswal is offline   Reply With Quote

Reply

Tags
body force, pimplefoam, source term

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Body force at the cell face Souviktor Fluent UDF and Scheme Programming 0 March 31, 2009 09:54
DPM Body Force Sandilya Garimella FLUENT 1 April 8, 2008 04:30
Problems with SUPG body force term FEM question Main CFD Forum 0 January 21, 2006 18:51
how to include body force in cfx selvam R.P CFX 4 November 25, 2005 05:01
UDF for lift force on a bluff body sawa FLUENT 2 April 11, 2005 04:06


All times are GMT -4. The time now is 00:50.