Fluid mechanics: continuity equation

In order to overcome the problem in FLIP solver, I looked into fluid mechanics.

It turns out the condition for incompressibility is not correct, or to say, based on an incorrect assumption.

For continuous fluid, the governing equation is called Continuity Equation

I'll skip the proving part of this equation, lots of reading materials could be found online discussing proving of this equation. In symbolic form, this equation could be expressed as:

Compare to the incompressible confinement, we could see, it's just assuming material derivative for density is 0. Yet this assumption is not correct if fluid cells are marked by the particles: temporay incoherence would lead to a non-zero density material derivative.

One big problem associated with adding this term to the confinement is, we'll have to deal with both time derivative and spacial gradient in the Eulerian grid.

So one of the possible solutions would be: combine SPH and FLIP together in this step. The density carried by each particle provides sufficient information for the material derivative. I'll keep on reading and thinking about other possible solutions.

Solution for compressibility problem in Eulerian method

As I mentioned before, FLIP solver suffers from compressibility problem. In fact, all the grid-based method might be suffering from comressibility problem.

You might consider me to be naive, but I'll try to convince you:

The grid based solver calculate the velocity field based on two formula:

1st is the Navier Stokes equation, 2nd is the incompressible constraint.

For grid solver, all the calculations are done based on one simple assumption:
all the fluid cells have the same density as the rest density.

And this is the key why these solvers are compressible. Though, it's just a small difference, which is hard to tell visually. This problem would become apparent in FLIP solver: with the same particle input, using different grid resolution would lead to different result, some shrink the volume(high resolution), some increase the volume(low resolution).

I tried to associate a density coefficient for each cell for the pressure solve, but it does not contribute. If you combine the rhow and p term in the 1st equation as a whole, you'll find out although the density coefficient would affect the value for pressure, when calculating back to velocity field, this term would be eliminated.

So I realize the problem lies in the second formula. There should be something more than just velocity field to conserve the volume.

It turns out I'm right. Well it's a pity that I found out someone has already done research into this before, or maybe I could be the first one.
With a simple search I found 2 papers dealing with this problem:
This one by Nuttapong Chentanez and Matthias Muller,
and this one by Michael Lentine, Mridul Aanjaneya and Ronald Fedkiw.

What is interesting is that what I've been doing always match these guys research. The papers I've been referring are always done by these names.

So next step: reading these papers and integrate into the FLIP solver. This might be useful for SPH as well, but I need to take a further look.