[szynka12]

Hello everyone,

So I can't find an explanation for the way that momentum equation is coded in OF multiphase solvers. For example DPMFoam:

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
);

or twoPhaseEulerFoam:

U1Eqn =
(
     fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
   - fvm::Sp(contErr1, U1)
   + MRF.DDt(alpha1*rho1 + Vm, U1)
   + phase1.turbulence().divDevRhoReff(U1)
  ==
      ...
 );

What bugs me, is the Sp source in both equations. In the first case it is given as

Here is the issue, I do recognise this as what appears from the decomposition of the conservative form and it would be logical in my opinion if this was added. My reasoning looks as follows: we have an equation in conservative form, we change it to the non-conservative, so that we solve for instead of :


Instead, it is subtracted, and I don't get why. Additionally this is supposed to make the equations conservative (look at this commit by Henry, it describes the change of the UEqns.H in twoPhaseEulerFoam in OpenFOAM 2.3x). So the question is, can someone explain or provide sources that explain what this source term actually does and how does it make the equations conservative?