ИСПОЛЬЗОВАНИЕ СОСТАВНОГО ЯДРА В МЕТОДЕ SPH ДЛЯ МОДЕЛИРОВАНИЯ ДВИЖЕНИЯ СЫПУЧИХ МАТЕРИАЛОВ Текст научной статьи по специальности «Физика»

Вычислительные технологии, 2021, том 26, № 1, с. 50—61. © ФИЦ ИВТ, 2021 Computational Technologies, 2021, vol. 26, no. 1, pp. 50-61. © FRC ICT, 2021

ISSN 1560-7534 elSSN 2313-691X

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

D01:10.25743/ICT.2021.26.1.004

Использование составного ядра в методе SPH для моделирования движения сыпучих материалов

И. И. Потапов, О. В. Решетникова*

Вычислительный центр ДВО РАН, Хабаровск, 680000, Россия

*

Поступила 27 июля 2020 г., доработана 29 сентября 2020 г., принята в печать 9 октября 2020 г.

В работе для моделирования движения сыпучей среды используется метод сглаженных частиц. Для аппроксимации искомых функций предложено новое составное ядро малой связности. Основой для разработки ядра послужило требование к условию о сохранении плотности единичной БРН-частицы. Выполнение данного условия позволяет правильно моделировать поле плотности на границах расчетной области, а также в случаях структурных изменений каркаса гранулированных частиц сыпучей среды. Из анализа решения задачи гидростатики методом БРН получена оценка значения масштаба сглаживающей длины ядра для двумерного случая. Выполнен расчет процесса обрушения гранулированного "столба" и проведено сравнение полученных численных результатов моделирования с экспериментальными данными.

Ключевые слова: движение гранулированной среды, метод сглаженных частиц, составное ядро, связность ядра.

Цитирование: Потапов И.И., Решетникова О.В. Использование составного ядра в методе БРН для моделирования движения сыпучих материалов. Вычислительные технологии. 2021; 26(1):50-61. Б01:10.25743/1СТ.2021.26.1.004.

Введение

В научной литературе по проблематике метода сглаженных частиц (БРИ) [1-6] не наблюдается устоявшейся терминологии при описании функций ядер сглаживания. Поэтому в настоящей работе при описании БРН-метода предлагается использовать следующую систему терминов. Пусть радиусом сглаживания ядра называется величина Н, определяемая как расстояние от центра сглаживаемой частицы до границы, за которой функция ядра Ш(г,Н) принимает нулевое значение (рис. 1), т.е.

где г — расстояние между парой частиц; Шр(г, Н) — порождающая функция ядра.

Радиус сглаживания ядра Н традиционно [1-3] выражается через характерный масштаб частицы Н, называемый сглаживающей длиной ядра: Н = ЯН. Сглаживающая длина ядра Н определяется через номинальный физический диаметр <1 частиц среды, находящихся в нейтральном, плотно упакованном состоянии, а именно Н = Ы, где к —

Рис. 1. Функция ядра W(г,Н) Fig. 1. The kernel function W(r,H)

безразмерный масштаб длины сглаживания относительно диаметра частицы. В расчетах удобно использовать безразмерное расстояние между частицами q = r/h, 0 < q < R, и R = H/h — безразмерный радиус сглаживания ядра. Таким образом, функция ядра приобретает вид

Wpb(q,R), 0 < г < R, 0, q> R,

где Wbp(q, R) — порождающая функция ядра, получаемая из функции Wp(q,R) после выполнения подстановки г = qh. Здесь и далее, где это не будет вызывать неоднозначности, для определения функции ядра индекс b будет опускаться.

h

допустимо лишь в случаях, когда безразмерный радиус сглаживания ядра R равен

h

соотнести с шагом сетки, тогда параметр R эквивалентен размеру сеточного шаблона в масштабе h. Соответственно, чем больше R, тем выше может быть степень аппроксимации получаемого дискретного аналога |1|, Безразмерный масштаб сглаживающей длины к связан с размерностью задачи и зависит от формы функции ядра W(q, R) [2].

Традиционно при решении задач механики сплошной среды SPH-методом 11—31 используются ядра с безразмерным радиусом сглаживания R больше единицы. Популярны ядра с R =2, например суперядро Гаусса [1], би-сплайн Монагана [2], ядро Лю с коллегами |3|, ядро Дисбрана |4|, Дня получения большой гладкости решения используются также ядра со значением R больше двух, например усеченное ядро Гаусса [5] или ядро Морриса |6|,

Ядра с большим радиусом сглаживания R обладают рядом достоинств: обеспечивают порядок точности аппроксимации искомых функций не ниже второго |1|, дают меньшую чувствительность к численной неустойчивости, возникающей при расчетах функций давления |2|, Однако наряду с достоинствами большой безразмерный радиус сглаживания R при ядре имеет недостатки [1-6], особенно важные при моделировании гранулированных сред. Недостатки связаны с ошибками аппроксимации плотности, которые возникают:

1) па свободной поверхности сыпучей среды, что приводит к пефизичпой подвижности частиц;

2) границах скольжения, при движении по ним топких слоев частиц (двух-четырех частиц), что вызывает эффекты пефизичоского отталкивания потока частиц от границ скольжения;

3) при движении сыпучей среды, приводящем к структурным изменениям ее каркаса, что искажает работу реологических моделей среды.

Wb

Отдельно стоит отметить тот факт, что чем больше безразмерный радиус сглаживания Я при ядре, тем более трудоемкими становятся вычисления [1, 2] при решении конкретной задачи.

Первые два из вышеозначенных недостатков ядер с большим Я преодолеваются введением дополнительных компенсирующих алгоритмов при вычислении плотности на границах расчетной области, например фильтров Шепарда [1], Все указанные ошибки переходят в разряд некомпенсируемых при моделировании гранулированных сред, поскольку ядра с большим Я не позволяют корректно описать движение самостоятельных гранул среды без введения дополнительных поправок в реологическую модель.

Известны попытки реализации в БРН компактных ядер с Я =1и к = 1, например ядра Люси [7] или ядра Мюллера с коллегами [8], Такие компактные ядра позволяют при расчетах полностью устранить ошибки вычисления плотности, связанные с неправильной аппроксимацией функции на границах области. Однако ядра Люси и Мюллера, удовлетворяя условию сохранения массы

J = 1, (1)

у

нарушают условие сохранения плотности для единичной частицы

тШ = ро, (2)

где т — масса частицы; р0 — номинальная плотность частицы.

Так, для ядер Люси и Мюллера характерна избыточная плотность в случае единичной частицы тШ(д = 0) > р0 и требуется увеличение масштаба длины сглаживания & для ее сохранения. Отметим, что при к > 1 происходит нарушение поля плотности на свободной поверхности, характерное для ядер с большим безразмерным радиусом сглаживания (Я > 1),

В данной работе предлагается новое составное ядро малой связности, которое удовлетворяет условию сохранения как массы (1), так и плотности единичной БРН-части-цы (2), Как показывают анализ и проведенные численные исследования, предлагаемое ядро хорошо подходит для моделирования движения гранулированных частиц, поскольку позволяет правильно моделировать поле плотности, в том числе и на границах расчетной области, а также в случаях структурных изменений каркаса гранулированных частиц.

1. Математическая постановка задачи

Рассмотрим плоское движение гранулированной тяжелой сыпучей среды. Расчетная область П с границами Г = Г и Г2 в начальный момент времени представлена па рис, 2, Движение сыпучей среды рассматривается в рамках модели плоской деформации [9] и описывается следующими уравнениями: — сохранения массы

V« (риа) = 0, (3)

— движения

= -5а(*VрР + рда, а, ¡3 =1, 2.

(4)

L

Рис. 2. Схема расчетной области задачи Fig. 2. The area for the computational scheme

Здесь р — мгновенная плотность гранулированной среды; иа — скорость среды; V« -оператор Гамильтона; Р — давление; да = (0,д) — проекции вектора ускорения свободного падения на координатные оси; 5а/3 — единичный тензор; а и @ — индексы соответственно горизонтальной и вертикальной координатных осей.

При решении задачи используется следующая гипотеза о связи между функцией давления Р и мгновенной плотностью среды р:

Р = тах (0, рос^/ч [(р/роУ - 1]) , (5)

где с8 — характерная скорость распространения упругой волны в моделируемой среде; 7 — показатель адиабаты.

Уравнения (3)-(5) замыкаются начальным

иа(х, 0) = 0, х е П (6)

и граничными условиями

иа(х,г) = 0, х е Г1, (7)

Р(х, ¿) = 0, х е Г2, (8)

где х = (ха) = (х, у) — координаты материальной точки.

2. Дискретизация задачи методом ЭРН

Выполним дискретизацию задачи (3)-(8) методом БРН. Будем, считать, что сыпучий материал определяется совокупностью N частиц. При аппроксимации плотности р^ и скорости и" для каждой г-й частицы были использованы следующие БРН-аппроксима-ции |1, 2|:

= ^ АзБз, VaAl = ^(А3 - А^Б,Vawij, (9)

где Аг и VaAi — соответственно аппроксимируемая функция и ее градиент; = Ш(д, Я) — значение функции ядра, определяемое через безразмерное расстояние д между парой частиц г и — двумерный "объем" — площадь частицы Sj = mj/(ру¿0)

(00 — толщина слоя гранулированной среды по неучитываемой координате трехмерного

, ттг ¿Шх™ - х"

пространства); = —-----градиент ядра в направлениях координатных

аг г^

С учетом (9) дискретный аналог уравнений (3)-(5) примет вид

^ = Е Т ш' (10)

3

°и? £ От2-. № + ^ + )+ (п)

dt doPj

Рг = max (0, p0c2 [(p/p0)7 - 1]) , (12)

, « -Uy)(xf -x°¡) < 0,

-'s

Vhc8 (< -U*)(x -x<*)

Cij = ^ 0.5(pi + p3) 4 + eerrh ' ^ ^

0, « — u?)(x? — x'j) > 0,

где Cij — искусственная диссипация по Моногану [2], ^ = 0.01 £err = 0.01, Движение частиц сыпучей среды дополняется уравнением

si /-у^

-ж=<■ <13>

Для решения задачи (10)—(13) выбрана схема Верле [8], устойчивость которой обеспечивалась выполнением условия Куранта -Фридрихса-Леви

At = Cu min ( -d—¡—-

j \cs + max |u"|

где Cu = 0.2 — коэффициент Куранта,

3. Составное ядро

В настоящей работе для аппроксимации искомых полей задачи (полей скорости и плотности) предложено новое составное ядро. Ядро имеет единичный безразмерный радиус сглаживания К = 1. Функция ядра определяется формулой

Ж, = +- [ (1 - ^ 0 ^ 1' (14)

У в 0, д> 1. { '

Производная /¿г для составного ядра определяется формулой

пг> Г О2 /4 20\

=^{? 4-;1 м

0, > 1 .

Здесь И — показатель размерности задачи (1,2 или 3); Ав и Вв — коэффищ^тг ядраи

( 7 л 12 / 105

его производной для задач с различной размерностью, Ав = —, 4--, 3 1 -

Д6' ^ V 32п

Вв = ( -,—,— ), Отметим, что для предлагаемого ядра выполняются условия (1)

3 7п п __

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

4. Калибровка составного ядра по градиенту давления

Дня калибровки составного ядра но градиенту давления в состоянии плотной упаковки, когда расстояние между всеми частицами равно г = с учетом К = Ы и д = 1/к перепишем формулу (15) в более удобном дня дальнейших выкладок виде:

dW dr

1

¥d

J w

I o,

w(k), 0,

Q < 1, q > 1,

(16)

где w(k)

9B

D

1 2 4

----+ -

4k \k2 к + 3

Выполним анализ влияния функции т(к) на точность определения градиента давления в плотно упакованной гранулированной среде дня двумерного случая. Из задачи гидростатики (рис. 3) в случае неподвижной среды получим формулу дня давления и градиента давления в слое:

/ ч

Р = РоЯУ, — = Род. (17)

йу

Здесь у — вертикальная координата, направленная вниз от уровня свободной поверхности среды.

Дня предлагаемого составного ядра при плотной упаковке частиц (рис. 3) выполняются условия К = 1, Н = кик = кй. Аналитическое давление на уровне п-го слоя частиц выражается как

Рп = родЬп, (18)

где Ь = ^л/3/2 — расстояние между слоями плотно упакованных частиц по вертикали.

Пусть центральная г-я частица в рассматриваемой плотной упаковке находится на уровне п, тогда окружающие ее ^'-е частицы находятся на уровнях п — 1, пап +1. Дня получения проекций сил па горизонтальную и вертикальную оси необходимо учи-

тывать, что шесть векторов

между центральной г-н частицеи и окружающими

ее j-ми частицами имеют следующие компоненты:

<У <У V _ Y

\_2,

V3 2

1, . (19)

Рп = p0gLn

Рис. 3. Схема плотной упаковки частиц в двумерном случае при решении задачи гидростатики Fig. 3. Scheme of dense packing of particles in the two-dimensional case when solving the problem of hydrostatics

Учитывая двумерноеть задачи Б = к2

Аха

Т.^т (р + Р

'-г' Р] ао

тз

Рз ао

из формул (4) и (5) получим

Ц/с)

-81кч род 1

, 1 /3 I

(п + п - ^ ( -2^ ) +

, л/3 .

+ (п + п - 1^2 +

+ (п + п + 1) ( -2, ~ ) +

, 1 ^

+ (п + п + 1) ( 2 , ^

о

или

¿Р1 Бт(к) /3

0ьха

Окончательно запишем

к2 ро9 2

(п + п - 1) (0, /3) + _ +(п + п + 1) (0,-/3) _

3т ( к)род (0,1).

т

аха

(1 + Д)(0,1)р0д, Д = 3Ц к) - 1.

(20)

При к = 1 из формулы (20) следует АР/Ау = 1.023род, т.е. ошибка аппроксимации градиента давления составляет около 2% (Д = 0.023),

Проанализируем характер изменения функции ошибки Д( к) при изменении плот-п

пользуя уравнение состояния (12) и гидростатический закон (18), получим выражение

п

Рп

7РЛ 1/1 (л , 9а У3 \ Н (л . ч1/7

' - '1+7-2 —п = ро (1 + ^арп)' '

(, + 1Рп\1/1 Л+ да /Д\

Ч1 + Юс2) = 1 + 7 ¿2" ~Гп)

где ар

д а л/3 с2 Т"

малый параметр, позволяющий выполнить по нему линеаризацию

п

Рп ~ Ро (1 + аРп).

(21)

Из выражения (21) следует, что при малых значениях а влиянием толщины столба гранулированного материала на решение задачи можно пренебречь,

В расчетах сыпучих материалов и несжимаемых жидкостей сглаживающая длина ядра к не изменяется, однако площадь частицы в силу (21) будет зависеть от глубины

п

Бп

т

т

о о о

(1 - арп) = к2(1 - арп), к

I т V Ао р о

(22)

С учетом (21) и (22) градиент давления в плотной упаковке частиц гранулированной среды примет вид

_ .ги(к)

V«Р = -(1 - арп)—^р0дЬ

(1 - арп)(п + п - ^ ( -2, ^ I +

1

+ (1 - арп)(п + п - 1) -, — +

2' 2

+ (1 - арп)(п + п + 1) ( -2,) +

/Л \/ л/3

+ (1 - арп)(п + п + 1) ( 2, —2"

)

(23)

Выполняя повторную линеаризацию полученного выражения по малому параметру а, получим

V аР « (1 - 4арп)3т(к)род(0,1). (24)

Из уравнения гидростатики АР ¡Ау = рд и тождеств (21), (24) следует

(1 - 4арп)?ш(к)р0д = (1 + арп)р0д.

По аналогии с несжимаемой моделью введем невязку ядра Дп, теперь уже зависящую от количества слоев частиц п в расчетной области:

Д.

3ЦА;)1 - - 1. 1 + арп

Из условия Д„ = 0 можно получить аналитическую зависимость

/

12-ш(к) + 1 ,__л/3 ,

3ш(к) - 1 ^ Уп = Т"^

(25)

связывающую между собой высоту расчетной области уп, выраженную через количество слоев среды п, и скорость распространения упругой волны в среде с3. Формула (25) уточняет оценку Моногана [9] и справедлива для произвольных одноевязных ядер.

с

в

5. Численные исследования

Верификация модели движения гранулированной сыпучей среды, описываемой уравнениями (10)—(13) в БРН-реализации, проводилась на модельной задаче об обрушении гранулированной колонны (рис, 4), Для анализа результатов вычислений использовались оригинальные эксперименты по обрушению колонны, составленной из стеклянных шариков с диаметром <1 = 0.016 м, В расчетах принимались следующие физико-механические и геометрические параметры: р0 = 2600 кг/м3 — плотность стекла, с3 = 60 = 1 а = 0.16 м и Ь = 0.138 м — соответственно ширина и высота колонны, сложенной из стеклянных шариков. На рис, 4 приведен один из экспериментов: показаны результаты обрушения колонны из состояния начальной укладки шариков (рис, 4, о) до состояния покоя после удаления стенки и их рассыпания (рис, 4, б).

Рис. 4. Вид гранулированной колонны из стеклянных шариков в состоянии начальной укладки (а) и после обрушения в состоянии покоя (б). Контрастной жирной линией неподвижные шарики основания отделены от тела колоны

Fig. 4. The granulated column of glass balls view in the initial state of laying (a) and at rest after the collapse (6)

Рис. 5. Вид гранулированной колонны в состоянии начальной укладки SPH-чаетиц (а) и после обрушения в состоянии покоя (б). Цветовая гамма частиц отражает ноле давления Fig. 5. View of a granular column in the state of initial packing of SPH particles (a) and in a state of rest after collapse (6)

Результаты расчетов SPH-методом для начальной укладки частиц, аналогичной физическому эксперименту, приведены па рис. 5. Здесь помимо пепосредстешюго расположения частиц визуализовапо получаемое при расчете ноле давлений в состояниях перед открытием перегородки (рис. 5, о) и после обрушения частиц (рис. 5, б). Несмотря па достаточно грубую аппроксимацию в "классическом" смысле — всего десять частиц в каждом из одиннадцати слоев, полученное ноле давлений согласуется с теоретическим гидростатическим давлением и имеет максимальную относительную погрешность в 22 % в угловых зонах и около 5 % па вертикальной оси симметрии колоны. Аномалии в угловых зонах вызваны различием в виде граничных условий па боковых границах, где используются подвижные зеркальные фиктивные частицы, и дне .нотка, где применялись неподвижные регулярно расположенные на расстоянии d частицы, имитирующие дно .нотка. При этом дня донных частиц выполнялась линейная аппроксимация давления из расчетной области.

Из сравнения полученных численных и экспериментальных данных можно сделать вывод об их хорошем согласовании, т. е. использование предложенного составного ядра при моделировании гранулированных сред является обоснованным.

Некоторое расхождение результатов, возможно, связано с несовпадением размерности и формы гранул: реальные сферические частицы в численном эксперименте представлялись цилиндрами, проекции которых на фронтальную плоскость образовывали двумерные SPH-чаетицы,

Заключение

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

Благодарности. Работа выполнена при частичной поддержке РФФИ (грант № 18-05-00530 А).

Список литературы

fl] Liu G.R., Liu М.В. Smoothed particle hydrodynamics: A meshfree particle method. World Scientific; 2003: 449.

[2] Monaghan J.J. Particle methods for hydrodynamics. Computer Physics Reports 3. North Holland, Amsterdam; 1985: 71-124.

[3] Liu M.B., Liu G.R. Smoothed particle hydrodynamics (SPH): An overview and recent developments. Arch Comput Methods Eng. 2010; 17:25-76.

[4] Desbrun M., Cani M.P. Smoothed particles: A new paradigm for animating highly deformable bodies. Eurographics Workshop on Computer Animation and Simulation (EGCAS). 1996: 61-76.

[5] Monaghan J.J., Gingold R.A. Shock simulation by the particle method SPH. Journal of Computational Physics. 1983; 52(2):374-389.

[6] Morris J.P., Fox P.J., Zhu Y. Modelling low Reynolds number incompressible flows using SPH. Computational Physics. 1997; 136:214-226.

[7] Lucy L.B. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal. 1977; 82(12):1013-1024.

[8] Muller M., Charypar D., Gross M. Particle-based fluid simulation for interactive applications. Proceedings of 2003 ACM SIGGRAPH Symposium on Computer Animation. 2003; 1540-159.

[9] Monaghan J.J. Simulating free surface flows with SPH. Journal of Computational Physics. 1994; 110: 399-406.

Вычислительные технологии, 2021, том 26, № 1, с. 50-61. © ФИЦ ИВТ, 2021 ISSN 1560-7534

Computational Technologies, 2021, vol. 26, no. 1, pp. 50-61. © FRC ICT, 2021 elSSN 2313-691X

MATHEMATICAL MODELLING

DOI: 10.25743/ICT.2021.26.1.004

Usage of composite kernel in the SPH method for simulating the movement of granular materials

Potapov Igor I., Reshetnikova Olga V.*

Computing center FEB RAS, Khabarovsk, 680000, Russia

*

Received July 21, 2020, revised September 29, 2020, accepted October 9, 2020

Abstract

The purpose of the study is to improve the practice of the SPH methodology which is applied for modelling of movement in the various media. The basis of the SPH-approximation of the function fields is formed by the forms of the smoothing kernel and its derivatives. Popular forms of smoothing kernels are characterized by the presence of significant fatal approximation errors when modelling granular media.

Methodology. The state of granular medium is described by the classical motion and mass conservation equations. Each granule of the medium corresponds to a separate SPH particle. To approximate the density and pressure fields in the SPH particle, a new combination of the smoothing core and its first derivative forms is proposed.

Results. The proposed new composite core fulfills the conditions of mass conservation and density recovery in the particle during SPH modeling. It is shown that the new composite core is characterized by a minimum error of pressure gradient approximation — about 2 %. A new estimate for the velocity of propagation of an elastic wave in a medium, sufficient to obtain a correct numerical solution, is proposed. A comparative analysis of the obtained solutions with experimental data is made.

Findings. The proposed composite shape of the smoothing kernel allows correct simulation of the motion of a granular medium by the SPH method. Its compactness (unit smoothing radius and unit smoothing length) makes it possible to correctly reconstruct the density field at the boundaries of the computational domain and in cases of structural changes in the framework of the granular medium. The numerical solution of the problem of the collapse of a column of granules obtained using the proposed composite core shows good agreement with experimental data.

Keywords: smoothed particle hydrodynamics, kernel approximation, simulation of granular material.

Citation: Potapov I.I., Reshetnikova O.V. Usage of composite kernel in the SPH method for simulating the movement of granular materials. Computational Technologies. 2021; 26(1):50-61. DOI: 10.25743/ICT.2021.26.1.004. (In Russ.)

Acknowledgements. The paper was funded by RFBR according to the research project No. 18-05-00530 A.

References

1. Liu G.R., Liu M.B. Smoothed particle hydrodynamics: A meshfree particle method. World Scientific; 2003: 449.

2. Monaghan J.J. Particle methods for hydrodynamics. Computer Physics Reports 3. North Holland, Amsterdam; 1985: 71-124.

3. Liu M.B., Liu G.R. Smoothed particle hydrodynamics (SPH): An overview and recent developments. Arch Comput Methods Eng. 2010; 17:25-76.

4. Desbrun M., Cani M.P. Smoothed particles: A new paradigm for animating highly deformable bodies. Eurographics Workshop on Computer Animation and Simulation (EGCAS). 1996: 61-76.

5. Monaghan J.J., Gingold R.A. Shock simulation by the particle method SPH. Journal of Computational Physics. 1983; 52(2):374-389.

6. Morris J.P., Fox P.J., Zhu Y. Modelling low Reynolds number incompressible flows using SPH. Computational Physics. 1997; 136:214-226.

7. Lucy L.B. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal. 1977; 82(12):1013-1024.

8. Muller M., Charypar D., Gross M. Particle-based fluid simulation for interactive applications. Proceedings of 2003 ACM SIGGRAPH Symposium on Computer Animation. 2003; 1540-159.

9. Monaghan J.J. Simulating free surface flows with SPH. Journal of Computational Physics. 1994; 110: 399-406.