CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Add a source term to simpleFoam

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

Like Tree3Likes
  • 1 Post By markc
  • 2 Post By psosnows

Reply
 
LinkBack Thread Tools Display Modes
Old   June 10, 2009, 15:16
Default Add a source term to simpleFoam
  #1
Senior Member
 
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 8
markc is on a distinguished road
Hello All,

I want to add a force vectorField as body force to the equations in simpleFoam. These forces do not change in time, it will simply be a vectorField with most cells set at zero and some will have a non-zero force vector.
Is it just as simple as putting:
>>>
eqnResidual = solve
(
UEqn() == -fvc::grad(p) (new part: + source vector)
).initialResidual();
<<<

Probably there is more involved. I searched this forum and surprisingly, similar questions were posted a couple of times but none of them were really answered. Hopefully this time some expert passes by?

Brgds,

Mark
Tushar@cfd likes this.
markc is offline   Reply With Quote

Old   June 10, 2009, 20:09
Default
  #2
Senior Member
 
Pawel Sosnowski
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 105
Rep Power: 9
psosnows is on a distinguished road
it is as simple as that

you just have to make sure that the units match.

the following lines are just because I may be tired,
so if you find those advice to direct, please forgive me

lest say you want to read this field from "constant" folder.
in createFields.H you just add:

Code:
volVectorField my_field_name
(
IOobject(
"my_field",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);
and now it "should" read the field "my_field" from constant folder and allow to use it in the code as such:
Code:
eqnResidual = solve
    (
        UEqn() == -fvc::grad(p) + my_field
).initialResidual();
the dimensions may not match, you maybe will have to divide the field by rho or something
you can do it in the code
Code:
eqnResidual = solve
     (
         UEqn() == -fvc::grad(p) + my_field/rho
).initialResidual();
or manualy while inputing data.

hope it helps you a bit.

Pawel
Mojtaba.a and Detian Liu like this.
psosnows is offline   Reply With Quote

Old   June 11, 2009, 03:40
Default
  #3
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,594
Rep Power: 24
ngj will become famous soon enoughngj will become famous soon enough
Hi Mark and Pawel

I do not agree that it is as simple as that, since if "momentumPrediction" is set to "no", then the momentum equation is not solved, and hence the body force is not taken into consideration. Thus my experience is that in the pEqn.H one need to modify the momentum corrector as follows:

U -= (fvc::grad(p) - my_field / rho )/AU;

However, you might need to consider how this contribution is added. I have used a constant in space, varying in time body force, thus the divergence of the momentum equation gives a 0 (zero) contribution to the pressure equation. As your bodyforce in varying in space, the term will take all the way to the pressure equation, I am not certain that it is sufficient to do as above.

Best regards,

Niels
ngj is offline   Reply With Quote

Old   June 11, 2009, 05:14
Default
  #4
Senior Member
 
Mark Couwenberg
Join Date: Mar 2009
Location: Netherlands
Posts: 130
Rep Power: 8
markc is on a distinguished road
Hello Pawel and Niels,

First of all thanks for your contributions.
Niels, you explain about momentumcorrector but we are talking about simpleFoam, which always uses a momentum corrector, right? In that case, does it still make sense?
Next: now I modified the simpleFoam solver as described by both of you. I tried it on the pitzDaily example (a backward facing step) and add some force on the upperwall, directing downwards. The liquid is not moving downward at all. While writing this I realise this might be caused by the fact that the fluid is incompressible. Does exerting a force on a wall makes sense? I better take some time to test it with a bodyForce on the internal fluid (what is what I actually want to do anyway).
Finally dimensions: Pawel, you are right. I need to supply force per unit volume and divided by rho, which gives: F/V/rho = (kg m s-2) * (m-3) * (kg m-3)^-1 = m/s2.

Once again thanks for your help,

Mark
markc is offline   Reply With Quote

Old   June 11, 2009, 06:02
Default
  #5
Senior Member
 
Pawel Sosnowski
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 105
Rep Power: 9
psosnows is on a distinguished road
Hello Niles and Mark,

Niles, I assumed that "momentumPredictor" is set to "true".
Now I have a question- if we assume that we solve the momentum equation, and that we add, for example, a gravity mass force, we get the UEqn like:
Code:
    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );
the momentum prodictor is set to "true" so we solve the U equation,
Code:
    eqnResidual = solve
    (
        UEqn() == -fvc::grad(p) + GravForce
    ).initialResidual();
we get the intermediet field and start calculating the pressure filed. To do that we use central coefficients of the UEqn:
Code:
volScalarField AU = UEqn().A();
after that we acquire a pressure filed which should be ok. And now: since the pressure equation was solved using velocity affected by the force, don`t you think it is enaugh to correct the velocity with newly acquired pressure gradient? Could you give a comment about that, I got a bit confised...

Mark,
first- check the boundary type- it should be "fixedValue" (I found myself forgeting about it several times)
As for the second thing- I belive that by adding this additional field you simply introduse a velocity source term to the UEqn. For me it seems that by putting some value on the wall would be just like applying a velocity inlet on that wall. (But it would be better if someone veryfied my statement).

Regards,
Pawel
psosnows is offline   Reply With Quote

Old   June 12, 2009, 04:29
Default
  #6
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,594
Rep Power: 24
ngj will become famous soon enoughngj will become famous soon enough
Hi Pawel

In the pressure equation, U is estimated as

U = UEqn().H()/AU;

hence based on all the terms in the UEqn, thus not on neither grad(p) nor my_fields. This means that the pressure correction to U does not include any of these, and therefore you need to have it in the second term as well.
One of my colleagues has been calculating laminar and turbulent oscillating wave boundary layers using OpenFOAM and he achieve perfect comparison with experimental/theoretical data using the method I have described above.

Best regards,

Niels
ngj is offline   Reply With Quote

Old   June 12, 2009, 05:39
Default
  #7
Senior Member
 
Pawel Sosnowski
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 105
Rep Power: 9
psosnows is on a distinguished road
Thenks Niles for pointing that out.

Fortunately I added my source terms into velocity equation so they were included while pressure calculations. Now I will pay more attention to that matter.

Regards,
Pawel
psosnows is offline   Reply With Quote

Old   June 12, 2009, 05:54
Default
  #8
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,594
Rep Power: 24
ngj will become famous soon enoughngj will become famous soon enough
Ah, that explains why you didn't need to add it in the bottom.

Good luck to both of you and enjoy the up-coming weekend,

Niels
ngj is offline   Reply With Quote

Old   September 9, 2009, 05:50
Default
  #9
New Member
 
nick stoppelkamp
Join Date: Jul 2009
Posts: 4
Rep Power: 7
nowhere is on a distinguished road
Hi,

I am also trying to bring in a source term in my simpleFoam solver, following the code snipped below.

Quote:
Originally Posted by psosnows View Post
Code:
eqnResidual = solve
     (
         UEqn() == -fvc::grad(p) + my_field/rho
).initialResidual();

But unfortunately the compiler complains that rho is not defined. So what do I have to do to get thr rho field?

Thanks and best regards,
Nick
nowhere is offline   Reply With Quote

Old   September 10, 2009, 10:55
Default
  #10
Senior Member
 
Pawel Sosnowski
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 105
Rep Power: 9
psosnows is on a distinguished road
Hello Nick!

The way I see it, there are two ways of including rho in your case:

1) when you create your field, divide it by rho at the beginning, set the dimensions as [myField/rho] and solve just:
Code:
 UEqn() == -fvc::grad(p) + myFieldDivRho
2) simply create rho field in the createFields.H file:
Code:
    volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        <initial_scalar_value_of_rho>
    );
but if I were you, I would use the first approach. On the other had, it all depends on what do you want to acheve.

Hope it helped a bit

best regards,
Pawel
psosnows is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
momentum source term zwdi FLUENT 13 December 5, 2013 18:58
How to add a flux source to SIMPLEfoam ehsan_vaghefi OpenFOAM Running, Solving & CFD 0 September 29, 2008 21:19
DxFoam reader update hjasak OpenFOAM Post-Processing 69 April 24, 2008 01:24
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51
Add source term in species equation zhou1 FLUENT 1 October 21, 2003 06:28


All times are GMT -4. The time now is 09:24.