momentum equation

Riccardo Maione's picture
Submitted by Riccardo Maione on Mon, 10/27/2014 - 13:17

Hello all,

I noticed something when I was reading the different solvers of CFDEM package, normaly the momentum equation is:

(
ddt(voidfraction, U)
+ div(phi, U)
+ divDevReff(U)
==
+ voidfraction*grad(p)
+ Source
);

so from what I have read solvers do not have the voidfraction near grad(p), is it true or I did not notice something that would make things clearer?

Thanks as always

Forums:

Riccardo Maione's picture

Riccardo Maione | Mon, 10/27/2014 - 13:36

sorry as a matter of fact i misread the article of zhou et al., for me model A was set III i am sorry for the inconvenience

Bruno | Mon, 03/02/2015 - 21:43

Indeed,

Only form A (set II of Zhou 2010) has a voidfraction * grad(p). Set I (B full) and Set III only have grad(p) and the void fraction action on the pressure is felt via the gradP and viscForce. In model A, the pressure is a lagrange multiplier to enforce conservation of mass of the fluid phase. In model B, pressure is a lagrange multiplier to enforce mass conservation of the suspension.

Cheers,
Bruno

hunger's picture

hunger | Thu, 10/29/2015 - 13:15

Hello!

I had a look at the difference between cfdemSolverPiso and cfdemSolverPisoScalar. Beside the temperature, I have recognized a difference in the momentum equation.

cfdemSolverPiso (here I removed the particleCloud.divVoidfractionTau-expression for better comparison)
======================================================================================
fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U) + fvm::Sp(fvc::ddt(voidfraction),U)
+ fvm::div(phi,U) + fvm::Sp(fvc::div(phi),U)
+ turbulence->divDevReff(U)
==
- fvm::Sp(Ksl/rho,U)
);
======================================================================================

cfdemSolverPisoScalar
======================================================================================
fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
- fvm::Sp(Ksl/rho,U)
);
======================================================================================

I was wondering where the fvm::Sp()-expressions in cfdemSolverPiso come from? When having a look into some papers, the equation of cfdemSolverPisoScalar would appear sufficient to me.
But then I realized that in the new version of CFDEM 2.9.0, this difference has disappeared (my example is from CFDEM 2.7.1). On the contrary, in both versions, the UEqn of cfdemSolverIB looks like the one of cfdemSolverPisoScalar.
The only idea that I could come up with was connected with the resolved and un-resolved approach, but I didn't find a satisfying explanation. Any suggestions?

Best regards
Harald

NTT1508 | Wed, 03/30/2016 - 03:42

Hi Harald,

In fact I found a difference between 2.70 and 2.90. While the 2.7.0 implements the solution as you shown above, the 2.9.0 uses exactly the same numerical implementation but the possitive +fvm::SP(_) has been changed to negative -fvm::SP(_) .
================================================================
fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U) =====> (This is: +fvm::SP(_) in the 2.7.0)
+ fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
// + turbulence->divDevReff(U)
+ particleCloud.divVoidfractionTau(U, voidfraction)
==
- fvm::Sp(Ksl/rho,U)
);
=======================================================================

I have checked the results from running the 2.7.0 and 2.9.0 on the same case, they are quite different. Do you have any ideas for this?

Thank youk

hunger's picture

hunger | Wed, 03/30/2016 - 10:43

Hi NTT1508!

If you follow the link to the OpenFOAM homepage I have posted earlier, you can read about the theoretical background of this numerical trick.

Concerning your question: They say that in principal this term should not be included for transient simulations, but only for steady-state, because in the end it should be 0 anyway. If you cannot ensure a steady-state solution for every single time-step of your transient simulation, this term is not equal to 0 and can influence the result of a time-step solution. In a transient case, these differences sum up and consequently lead to varying results.

However, this is just aguess. I did not investigate the consequences of this term on the transient PISO solver.

Best regards
Harald

NTT1508 | Mon, 04/04/2016 - 00:59

Hi Harald,

Thank you, it sounds reasonable, though my investigations are truly steady-state (low hydraulic gradient through a stable porous medium). I will note on this. You can try by making a simple particulate bed and let a laminar flow through that.

Regards,

hunger's picture

hunger | Tue, 11/10/2015 - 15:44

I could not find an explanation for my previous question yet. Any help is still appreciated!

Best regards
Harald

cgoniva's picture

cgoniva | Tue, 01/26/2016 - 22:50

Hello Harald,

Thank you for digging into the numerics.

I guess you have also found out that the sign of the fvm::Sp(...) was wrong, which is corrected in latest versions.

Best regards,
Chris

NTT1508 | Mon, 04/04/2016 - 01:50

As the numeric method used in OpenFOAM shows and Chris pointed, the negative -fvm::Sp(...) was wrong in 2.7.0.

sbateman | Thu, 02/02/2017 - 22:07

I'm going to resurrect this thread because I don't understand the purpose of the fvm::Sp(...) terms in the LHS of the cfdemSolverPiso UEqn:

fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+ fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
// + turbulence->divDevReff(U)
+ particleCloud.divVoidfractionTau(U, voidfraction)
==
- fvm::Sp(Ksl/rho,U)
+ fvOptions(U)
);

This page explains why the terms should be included for steady-state flows:
http://openfoam.org/release/2-2-0/numerics-boundedness/
But the fvm::Sp(...) terms are not included explicitly in the UEqn of the OpenFOAM solvers. Instead, there are "bounded" ddtSchemes and divSchemes that are specified in the fvSchemes file. For example, here is from pisoFoam:

fvVectorMatrix UEqn
(
fvm::ddt(U) + fvm::div(phi, U)
+ MRF.DDt(U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);

So I wonder why these terms are included in cfdemSolverPiso.
The simulations I am running are for transient flow, so I don't think I need these "bounded" schemes. Although, I haven't tested whether they actually affect my results or not.