[markc]

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

[psosnows]

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:

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:

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

eqnResidual = solve
     (
         UEqn() == -fvc::grad(p) + my_field/rho
).initialResidual();

or manualy while inputing data.

hope it helps you a bit.

Pawel

[ngj]

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

[markc]

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

[psosnows]

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:

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

the momentum prodictor is set to "true" so we solve the U equation,

    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:

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

[ngj]

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

[psosnows]

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

[ngj]

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

[nowhere]

Hi,

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

Originally Posted by psosnows

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

[psosnows]

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:

 UEqn() == -fvc::grad(p) + myFieldDivRho

2) simply create rho field in the createFields.H file:

    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