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
Riccardo Maione | Mon, 10/27/2014 - 13:36
sorry as a matter of fact i
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,
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 | Thu, 10/29/2015 - 13:15
Additional expressions in the UEqn
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?
NTT1508 | Wed, 03/30/2016 - 03:42
Difference in the 2.7.0 and 2.9.0
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 | Wed, 03/30/2016 - 10:43
Hi NTT1508!
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.
NTT1508 | Mon, 04/04/2016 - 00:59
Theoretically make sense
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 | Tue, 11/10/2015 - 15:44
I could not find an
I could not find an explanation for my previous question yet. Any help is still appreciated!
hunger | Tue, 01/26/2016 - 08:36
Explanation found!
I think I have found a nice explanation of these terms on the OpenFOAM website:
http://www.openfoam.org/version2.2.0/numerics.php
cgoniva | Tue, 01/26/2016 - 22:50
momentum equation fvm::Sp(...) term
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
2.7.0 was wrong
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
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 theUEqn
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.