[DuanYB]

Hello Fomers, Good days. Help!!!

I have been learning openfoam for one year, and I am still a novice. There is a problem that has troubled me for more than half a year and has not been resolved. I would like to ask the predecessors in the forum about the influence of the mass source term on the momentum equation and the encoding. By the way, the solver I mainly use is simpleFoam now.

In a stendy, inviscid and incompressible problem, if there is not a mass source term in a continuity equation，one has
div(U)=0 (1)
div(U*U)=-grad(p) (2)
If there is a mass source term, one has
div(U)=S (3)
div(U)*U + U*grad(U)=-grad(p) (4)

Because if there is not a mass source term, div(U*U)=div(U)*U+U*grad(U) , so equation(2) and equation(4) is equal. Even with the source term, the above equation should be also true, but I don't know if the above equation is true after the equation is discretized in openfoam, because when I don't change the momentum equation and just add a divU, openfoam always reports an error. So my questions are
(i) If div(U) is not equal to 0 in openfoam, is the formula div(U*U)=div(U)*U+U*grad(U) still true?
(ii) Because the left hand two terms of formula (4) are nonlinear. So how to write a program in openfoam to realize formula (4)?

I really need help from all of you, even a small suggestion.

[Zena27]

Hi I am quite new to OpenFOAM too and I don't have a clear answer, but have you tried adding the mass source through the fvOptions? or you want to hard code it yourself?

[DuanYB]

Originally Posted by Zena27

Hi I am quite new to OpenFOAM too and I don't have a clear answer, but have you tried adding the mass source through the fvOptions? or you want to hard code it yourself?

Because my source term ia a vector, and all cells have a same source term, so I just change a poisson equation in the pEqn.H. Just like

fvScalarMatrix pEqn
(
  fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA) - scalarDivU
);

As for the fvOptions you said, I think it is mainly to add the momentum source term when building the momentum equation.

[Zena27]

Okay. I find this problem interesting. What exactly the kind of error that you encounter when not editing the momentum equation? Can you post the error message here?

[Zena27]

Hello,
I noticed sth I don't know it's right or wrong. Please correct me if it doesn't make sense.

But I revised How simpleFoam solve for the velocity field. So first if we have the Nacier-Stokes equation:

ρ DU/Dt = -∇ p + ∇ . τ + ρ g

. What actually happens is that we assume stress tensor is a function of U, no gravity force, and steady state, so we end Up with this simple equation:

MU = -∇ p

So basically, lets ignore about how it solves for the U exactly. but it first creates the MU Matrix through the tUEqn object in the UEqn file, and then it equates it to the -∇ p in the momentum predictor step.

solve(UEqn == -fvc::grad(p))

so what I am thinking is to try to add the mass source in the momentum predictor step as this

solve(UEqn == -fvc::grad(p) + mass_Source)

, and then add it too in the pEqn file as you already did.
because the edit I added is the addition of the mass source to the momentum equation and in the pEqn is it's effect on the continuity.

Note: be carreful with the units, and you need to enable the momentum predictor option.

Also, I used this source to revise how the simple algorithm work: https://medium.com/@mustafaabbs2/the...d-61fcc54ab24d

[DuanYB]

Originally Posted by Zena27

Hello,
I noticed sth I don't know it's right or wrong. Please correct me if it doesn't make sense.

But I revised How simpleFoam solve for the velocity field. So first if we have the Nacier-Stokes equation:

ρ DU/Dt = -∇ p + ∇ . τ + ρ g

. What actually happens is that we assume stress tensor is a function of U, no gravity force, and steady state, so we end Up with this simple equation:

MU = -∇ p

So basically, lets ignore about how it solves for the U exactly. but it first creates the MU Matrix through the tUEqn object in the UEqn file, and then it equates it to the -∇ p in the momentum predictor step.

solve(UEqn == -fvc::grad(p))

so what I am thinking is to try to add the mass source in the momentum predictor step as this

solve(UEqn == -fvc::grad(p) + mass_Source)

, and then add it too in the pEqn file as you already did.
because the edit I added is the addition of the mass source to the momentum equation and in the pEqn is it's effect on the continuity.

Note: be carreful with the units, and you need to enable the momentum predictor option.

Also, I used this source to revise how the simple algorithm work: https://medium.com/@mustafaabbs2/the...d-61fcc54ab24d

Hi, I read your reply carefully. One thing I don't know is right, I think what should be added to the momentum equation should be the momentum source term, not the mass source term. So if your "mass_source" refers to the momentum source term, it should be added to tUEqn (the fvOption in the original equation also refers to the momentum source term); if it refers to the mass source term, it should be added to pEqn.
The above is my own understanding after checking some information recently, please correct me if I am wrong. In addition, if you have WeChat or twitter or other instant messaging tools, welcome to send it to me, and we can have further discussions.

[Zena27]

yep I am mistaken, the thing added is the momentum source term not the massSource as I written. But be aware of sth that the mass source term adds momentum to the flow and creates momentum source that you need to add. So, in my opinion it's not just the mass source that should be added.

You said that when you added the mass source term you got an error. Also, are you sure that it's not dimensions error or unit error.

Moreover, I thought about adding the momentum source to the momentum predictor if you don't want to add it by fvOptions. My thought about that is that it's not the best practice but it should work because it satisfies the equation. It's basically you have

MU = -∇ p + Forcing term

and the forcing term here is the momentum source due to the mass source.

[DuanYB]

Originally Posted by Zena27

yep I am mistaken, the thing added is the momentum source term not the massSource as I written. But be aware of sth that the mass source term adds momentum to the flow and creates momentum source that you need to add. So, in my opinion it's not just the mass source that should be added.

You said that when you added the mass source term you got an error. Also, are you sure that it's not dimensions error or unit error.

Moreover, I thought about adding the momentum source to the momentum predictor if you don't want to add it by fvOptions. My thought about that is that it's not the best practice but it should work because it satisfies the equation. It's basically you have

MU = -∇ p + Forcing term

and the forcing term here is the momentum source due to the mass source.

Unfortunately I don't have any, I have a linkedIn account Zeinab Abosedira.

Hey, If we just want to add some kind of force, I think direct addition is correct.
In addition, I would like to ask, you said "But be aware of sth that the mass source term adds momentum to the flow and creates momentum source that you need to add." Can you give me some information or some references? I would like to understand how to add momentum due to the mass source term. Thank you very much!

In addition, the error I made may have nothing to do with the unit. I declared a volScalarField constant as div(U) in "createFields.H", and wrote the correct unit in this file. Dimensions should be fine too, since fvc::div(phiHbyA) is a volume scalar, as is div(U).

[Zena27]

Off course I can search and send you some source. I just have a question can you give me a brief description of your problem, and can you post the error you got here?

[DuanYB]

Originally Posted by Zena27

Off course I can search and send you some source. I just have a question can you give me a brief description of your problem, and can you post the error you got here?

My problem is how to solve an incompressible inviscid steady fluid flow problem. I'm pretty sure my continuity equation is div(U)=0, but not so sure about the momentum equation.

Next I paste my modified simpleFoam code。

1. In the createFields.H, I added the following code:

volScalarField scalarDivU
(
  IOobject
  (
    "scalarDivU",
    runTime.timeName(),
    mesh,
    IOobject::MUST_READ,
    IOobject::AUTO_WRITE
  ),
  mesh
);

2. In the pEqn.H, I added the scalarDivU, like this:

fvScalarMatrix pEqn
(
fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA) - scalarDivU
)

3. I also added a file named scalarDivU in the 0 folder. the main content is

dimensions [0 0 -1 0 0 0 0];
internalField     uniform 1.12;
boundaryField
{
  SURF_EXTERIOR_OUTLETS
  {
    type  fixedValue;
    value uniform 1.12;
  }
  SURF_EXTERIOR_INLETS
  {
    type  fixedValue;
    value uniform 1.12;
  }
}

That is all I changed in the simpleFoam. When I use clion debug, The program stopped at the

solve(UEqn == -fvc::grad(p))

of Time=2(Time=1 no error). The program error is reported as

p = expression failed to parse:error: <user expression 36>:1:1: use of undeclared ideentifier 'p'p^

signal = SIGFPE (signal SIGFPE: floating point divide by zero)


By the way, I really donot know which value is zero??? Maybe the UEqn??? But why???

[Tobermory]

Originally Posted by DuanYB

So how to write a program in openfoam to realize formula (4)

Your formula 4, from your first post is in a non-conservative form, and so cannot be discretised simply with the finite volume approach. Apologies, but I haven't read through the rest of the posts in the thread, but it strikes me from your first post that you may be trying to do something impossible.

I think you probably need to read up more about the basic OpenFOAM solver, and its FV approach - search online for Hrv Jasak's thesis, for example - there's a good description in there, I seem to remember. Good luck!

[DuanYB]

Originally Posted by Tobermory

Your formula 4, from your first post is in a non-conservative form, and so cannot be discretised simply with the finite volume approach. Apologies, but I haven't read through the rest of the posts in the thread, but it strikes me from your first post that you may be trying to do something impossible.

I think you probably need to read up more about the basic OpenFOAM solver, and its FV approach - search online for Hrv Jasak's thesis, for example - there's a good description in there, I seem to remember. Good luck!

Thanks for your reply, I will study it seriously. Do you mean that openfoam can only establish equations in conservative form? Now I know that my efforts may be going in a wrong direction. Thank you again.

In addition, I would like to ask a little question, if the mass source term is added to the continuity equation, does the momentum equation need to change some terms?

[Tobermory]

Do you mean that openfoam can only establish equations in conservative form?

Yes, since it is based on a finite volume approach, whereby the fundamental equations are integrated over a control volume, the grid cell. To do this, you need to take advatange of the Gauss divergence theorems etc, and the equations need to be in conservative form.

if the mass source term is added to the continuity equation, does the momentum equation need to change some terms

No - adding (or removing) mass from a cell does not (directly) change the momentum in the cell, it just affects the density. Any mass addition comes in with zero momentum and zero energy (if you are running a solver with a temperature equation) ... so be careful that this is what you want - often, eg for subgrid mass inlet, you need to apply source terms for mass, momentum & energy to reproduce the inlet stream properties properly.

[DuanYB]

Originally Posted by Tobermory

No - adding (or removing) mass from a cell does not (directly) change the momentum in the cell, it just affects the density.

OK! Thanks very very much !!!

I would like to ask one more question. If at the beginning of the system, the whole space is full of mass source items (that is, div(U) is equal to S everywhere in the solution space). And what I am solving is an incompressible steady-state problem, so will the density of the system change?

Thanks again for your thoughtful answer ！

[Tobermory]

Yes! Think about it from a physical perspective - it is the same as pumping air continuously into a container, where the walls of the container are rigid (i.e. the volume remains fixed). The mass is increasing, and so must be the density.

All the best, T

[DuanYB]

Yes! Think about it from a physical perspective - it is the same as pumping air continuously into a container, where the walls of the container are rigid (i.e. the volume remains fixed). The mass is increasing, and so must be the density.

Thanks for your reply !
Yes, you say the density does increase in this case. But what I want to solve is an incompressible problem. The detailed description of this problem is "a certain substance flows from the inlet to the outlet, and the container is filled with a mass source, during which the whole system is stable, constant temperature, and incompressible".
I made the density constant as a precondition. What I can think of is that adding the mass source term will definitely change the velocity of the fluid compared to no mass source term, but I don't know if this change will change the momentum equation.

According to your explanation, the momentum equation doesn't change even if the density doesn't change right?

Thanks again

[Zena27]

Originally Posted by Tobermory

No - adding (or removing) mass from a cell does not (directly) change the momentum in the cell, it just affects the density. Any mass addition comes in with zero momentum and zero energy (if you are running a solver with a temperature equation) ... so be careful that this is what you want - often, eg for subgrid mass inlet, you need to apply source terms for mass, momentum & energy to reproduce the inlet stream properties properly.

Hello, I need a question about this?
How the momentum is not changed. Isn't the mass we add has a velocity and this velocity? and then we have the momentum is basically mass*velocity. Also, concerning the energy equation, isn't the mass added has an internal energy thus it changes the energy too. Please, I need more explanation about this.

[Tobermory]

Apologies - I forgot you were dealing with an incompressible system. In that case, density is indeed constant and instead you will get an increase in velocity downstream to accomodate the extra mass flow. This will be accompanied by a pressure rise to counter the increased momentum of the flow, thereby leaving the momentum equation balance unchanged. This back-pressure is the resistance to you pumping the mass source in, and would in a real system affect the upstream conditions ... etc.

[Tobermory]

Originally Posted by Zena27

How the momentum is not changed.

Apologies - I meant the momentum balance (ie momentum plus body forces like pressure), not the momentum itself - Duan was wanting to know how this affects the momentum equation, and that was what I was trying to answer. Sorry for the confusion.

[Zena27]

Okay Now I understand thank you so much, and sorry DuanYB for confusing you earlier. :""