Моделирование жидкости несжимаемой жидкости

В моделировании жидкости на основе статья Стама, жидкость моделируется как сетка плотностей. Плотности «обычно принимают значение от нуля до единицы», но могут быть больше. Границы реализуются в основном так, как описано в статье:

Простой способ реализации внутренних границ — выделить булеву сетку.
который указывает, какие клетки заняты объектом или нет. Тогда мы просто должны добавить
некоторый код подпрограммы set_bnd () для заполнения значений для занятых ячеек из значений
их прямые соседи.

int surround = !bound[IX(i+1,j)] + !bound[IX(i-1,j)] + !bound[IX(i,j+1)] + !bound[IX(i,j-1)];
if (!surround) x[IX(i,j)] = 0;
else
x[IX(i,j)] = ((bound[IX(i+1,j)] ? 0 : x[IX(i+1,j)]) +
(bound[IX(i-1,j)] ? 0 : x[IX(i-1,j)]) +
(bound[IX(i,j+1)] ? 0 : x[IX(i,j+1)]) +
(bound[IX(i,j-1)] ? 0 : x[IX(i,j-1)])) / surround;

Метод плотности хорошо работает для сжимаемых жидкостей, таких как воздух, огонь или дым. Есть ли способ изменить граничные процедуры, чтобы плотность (ограниченная одной жидкостью) была ограничена значением, например 1? Это будет, скажем, клетка, полностью заполненная частицами воды Плотность, которая была бы больше единицы, должна была быть отброшена к соседним ячейкам. Stam перечисляет идеи для расширения, но не включает как:

Еще одним расширением является использование этого решателя в качестве основы для оживления потоков воды. В этом случае есть
две жидкости с различной плотностью: вода и воздух. Воздух обычно не моделируется, а решатель
труднее реализовать по следующей причине: домен водной жидкости меняется
время и нужно как-то отслеживать, и должны применяться правильные граничные условия
на интерфейсе. Водную область можно отслеживать с помощью частиц, которые просто перемещаются через
жидкость, как это сделано Фостером и Метаксасом [Foster96] или может отслеживаться с помощью комбинации
частицы и наборы уровней [Foster01, Enright02].