Calculation of the continuity density

[LasNikas]

So far, we calculated the ContinuityDensity with eq. 2.18

Adami 2012 uses eq. 2.17 and justify this with a stable simulation even when dealing with high density ratios.

In Monaghan 2005: "Both expressions vanish, as they should, when the velocity is constant. However, when the system involves two or more fluids with large density ratios in contact, the expression (2.17) with ρ in the summation is more accurate (Colagrossi 2004)."

I think we should consider this in the refactoring.

[efaulhaber]

We should definitely compare these properly.

Then I see two options:

1. Stick with the best one. Then we can ignore this for the refactoring and implement it afterwards.
2. Implement both and add an option to switch.

[svchb]

Check out this paper: Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations, J. Bonet *, T.-S.L. Lok

Here they claim in section 2.4.1 and 2.4.2. that you should use a slightly different formulation for the internal force for the continuity equation compared to the summation density case.

[efaulhaber]

Also note that there are different formulations for the momentum equation. I'm also planning to compare these soon. For example, DualSPHysics is using (p_a + p_b)/(rho_a rho_b) instead of our p_a / rho_a^2 + p_b / rho_b^2.

[svchb]

yes the paper claims that (p_a + p_b)/(rho_a rho_b) should be used with a continuity formulation andp_a / rho_a^2 + p_b / rho_b^2with summation density

[svchb]

Since they draw the conclusion purely on mathematical arguments it is probably true.

[LasNikas]

Adami et al produces very good results with an inter particle averaged pressure: with the number density formulation. They do not use any artificial viscosity or empirical parameters but impose a advection velocity. Ramachadran 2019 is also using this formulation for the EDAC scheme.

[svchb]

I fear that these different formulations are probably even case dependent.

[LasNikas]

Indeed. The TVF (transport velocity formulation) of Adami is only for internal flows. But Ramachandra promised a solution for every purpose. I'm very curious.

[svchb]

And in "Robustness and accuracy of SPH formulations for viscous flow" by Mihai Basa et al. they claim that this leads to an improvement:

[LasNikas]

Ramachandran uses also this average pressure, but claims that "that this results in significantly improved results that outperform the traditional TVF scheme."

[svchb]

Without a lot of reference cases this cannot be answered. Anyway this could be made into a paper.

[svchb]

I guess we should for now make these different formulations user selectable until we can draw a conclusion.

[efaulhaber]

Verdict of our Slack discussions: According to multiple papers, 2.17 is better for multi-phase flows. Apparently, it doesn't have any disadvantages apart from the negligibly higher computational cost.

The calculations by Bonet & Lok don't show that (p_a + p_b)/(rho_a rho_b) is the better pressure acceleration when using ContinuityDensity. They only show that this is the correct formulation for 2.17. When repeating their calculations with 2.18, they yield p_a^2 / rho_a + p_b^2 / rho_b, just like for SummationDensity. This is consistent with my findings in #294. Total energy is only conserved when using (p_a + p_b)/(rho_a rho_b) and 2.17 or p_a^2 / rho_a + p_b^2 / rho_b and 2.18.

I will therefore change the formulation to 2.17. I found that the order of operations (m_b * rho_a / rho_b vs rho_a / rho_b * m_b) made a 10% difference in the benchmark of interact! in Julia 1.9. I will leave this issue open to re-check this in future Julia versions.