[XorT]

Hello,


I'm trying to develop some prototype for a realtime 2d fluid dynamics simulation. The goal is to be fast rather than physically exact as long as it is visually plausible (use-case is computer-graphics). However the simulation needs to deal with high density ratios as it needs to handle water-air-interface dynamics in a convincing way.

The method described in the paper "A fast pressure-correction method for incompressible two-fluid flows" (https://doddm.com/publications/2014-jcp-md-af.pdf) appears very interesting as it can employ a FFT based FastPoisson solver to solve the pressure equation, which can achieve realtime performance on high-resolution grids using GPUs.



However I experience severe stability issues (blow ups).
My simulation setting is fairly simple: it is only 2D and the initial conditions are a drop of water (1000kg/m²) in the center with initial velocity in the y-direction on a background of initially still air (1kg/m²). As quantities I keep only cell-centered densities and staggered velocities. I use the MacCormack scheme for advection.



Without advecting the density- and velocity field the simulation is stable and the solved divergence-free velocitiy field stabalizes to a nice (static) approximation (in which mostly the air gets accelerated to match the water' velocity, while the water is barely slowed by the initially still air, due to the massive density ratio).
However if i do advect, even for really small velocities (far less than a grid cell), some velocities will exponentially grow and blow up over multiple iterations. In contrary this does not happen for a density ratio of 1 - in this case I get a nice fluid motion, but of course it behaves as one unified fluid due to equal densities.


Since the described system uses feedback it is quite hard to analyse what causes this instability, so I wanted to ask you first if you have experienced any stability/blowup issues using this modified pressure-solve-scheme or have worked with this method at all. I checked many times my code and can not spot any error :'(

I really hope you can help me out a little to solve this stability issue as this method seems to be quite ideally suited for my use-case and solutions which permit realtime performance are rare


Best regards