МОДЕЛИРОВАНИЕ СПЛОШНОЙ СРЕДЫ ПРИ ПОМОЩИ МЕТОДА SPH Текст научной статьи по специальности «Математика»

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

1. Хоар Ч. Взаимодействующие последовательные процессы: пер. с англ. - М., Мир, 1989. - 264 с., ил.

Волгоградский государственный университет, г. Волгоград

В настоящей работе рассмотрен метод гидродинамики гладких частиц в приложении к моделированию сплошных сред. Реализована и протестирована программная модель на основе метода SPH с визуализацией при помощи OpenGL.

Ключевые слова SPH, гидродинамика, моделирование.

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

2. Аппроксимация значений физических величин при помощи SPH 2.1. Вычисление значения атрибута и сглаживающая функция

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

Значение атрибута интерполируется посредством взвешенной суммы воздействий всех частиц и вычисляется по формуле

МОДЕЛИРОВАНИЕ СПЛОШНОЙ СРЕДЫ ПРИ ПОМОЩИ МЕТОДА SPH

© Кибко М.О.1

1. Введение

1 Аспирант очного обучения.

где AS - вычисляемое сглаженное значение атрибута;

r - позиция частицы в пространстве, м;

Aj - значение атрибута для j-ой частицы;

W(r - rj, h) - сглаживающая функция (функция ядра);

h - длина сглаживания, м.

Функция W(r, h) называется сглаживающей функцией ядра. Данная функция выбирается при моделировании и должна быть гладкой. При этом должны также выполняться следующие условия:

lim W (r - r',h)- dr'= ö (r — r') dr', (2)

jw (r',h) — dr' = 1, (3)

где S(r - r') - функция Дирака в точке r.

Фактически сглаживающая функция определяет вес влияния каждой частицы на вычисляемую величину. Чем ближе соседняя частица к той частице, для которой вычисляется величина атрибута - тем, соответственно, больше её влияние на результат и тем больше сглаживающая функция в этой позиции. При этом сглаживающая функция имеет такой параметр, как длина сглаживания h. Этот параметр определяет масштаб функции и соответственно разрешение, с которым происходит моделирование, так как основной вклад в вычисление величины вносят частицы на расстоянии не более чем 2h от частицы, для которой производятся вычисления. Другие частицы либо вносят незначительный вклад в вычисления либо не вносят его вовсе, в зависимости от выбранной сглаживающей функции.

Стоит отметить, что интерполяция значения атрибута может иметь второй порядок точности [1] только если сглаживающая функция нормализована и удовлетворяет условию

W (r, h) = W (—r, h). (4)

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

dx j р dx '

где х - величина, по которой вычисляется производная (как правило, пространственная координата).

Таким образом, градиент величины A можно представить как

A

VAS (r) = Yjmj-LVW(r -rj,h), (6)

J PL

а лапласиан величины А принимает форму

А

V2 А ( г ) = Х у2 ^ (Г - г, к). (7)

] Р]

2.2. Вычисление плотности и закон сохранения массы в 8РИ

Подставив плотность как вычисляемую величину в уравнение (1), получаем стандартный подход для вычисления плотности в методе 8РИ

Р

= X(Г - Г,к).

(8)

Данный метод позволяет гарантировать сохранение массы в модели, но создает ряд проблем, связанных с поверхностью жидкости или газа. Частицы, которые оказываются на поверхности или близко к ней имеют значительно меньшую плотность, чем частицы внутри жидкости (или газа) вследствие того, что они имеют значительно меньше частиц соседей. В последствии неверное вычисление плотности также приводит к неверному вычислению давления. Данное свойство производит артефакты поверхностного натяжения [2].

2.3. Закон сохранения импульса и вычисление ускорения частицы

Уравнение Навье-Стокса определяет закон сохранения импульса для метода 8РИ и позволяет вычислять изменение ускорения частиц во времени.

Уравнение формулируется следующим образом

+ ^ ) = -^ + р§ + (9)

где р- плотность частицы, кг/м3;

V - скорость частицы, км/ч;

/ - время, с;

р - давление в данной точке, Па;

g - ускорение свободного падения, ~9,8 км/с2;

[ - динамический коэффициент вязкости, Па-с.

Соблюдение закона сохранения импульса необходимо, чтобы моделируемая жидкость была несжимаемой, наряду с законом сохранения массы.

Использование лагранжевой гидродинамики (представление жидкости частицами) вместо эйлеровой (представление в виде сетки) упрощает уравнение (9). В левой части уравнения (9) сумма частной производной скорости по времени и члена, описывающего конвекцию, заменяется производной Лагранжа)

-= —(10)

О/

А так, как частицы движутся вместе с жидкостью, то производная Ла-гранжа может быть заменена производной скорости частицы по времени, что означает, что конвективный член уУу не нужен для моделей с частицами.

На правой стороне уравнения (9) имеются три силовых поля, моделирующие давление -Ур, внешние силы р и вязкость /1У2у. Сумма этих трёх силовых полей

/ = -Ур + pg + (11)

определяет изменение импульса рна левой стороне уравнения (9). Соответственно, для ускорения 1-ой частицы, получаем

а =-Т7 = — > (12)

Ш р

где /, - сумма силовых полей для данной частицы.

3. Реализация численной модели при помощи метода 8РИ

3.1. Моделирование сил, влияющих на ускорение частиц

Применение правила 8РИ, описанного в уравнении (1) к первому члену правой части уравнения (9) -Ур, описывающему давление, приводит к формуле, позволяющей вычислить силу давления для конкретной частицы

ГГШГе =-Ур (Т ) = -£ тур УЖ (п - Ту, к). (13)

у ру

Сила давления, вычисляемая таким образом - несимметрична, что влияет на качество моделирования, делая его менее точным, а модель - более неустойчивой. Наиболее простым способом симметризации уравнения (13) может служить следующий способ, который и был использован в данной работе

ГГШ" = -Е ту^^ УЖ (г - Ту, к), (14)

у 2ру

где р/ - значение давления для /-ой частицы.

Данный способ является подходящим, обеспечивая устойчивость модели и совсем незначительно увеличивая требования к вычислительным ресурсам. Сила давления, вычисляемая подобным образом, является симметричной, так как использует среднее арифметическое давлений взаимодействующих частиц.

Так как частицы имеют только три свойства, хранимых в памяти непосредственно: масса, позиция в пространстве и скорость, то давление в месте

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

В данной работе использовалось уравнение состояния, предложенное Мо-наганом [3], как наиболее подходящее для моделирования потоков жидкости

где В - постоянная давления, параметр моделирования, Па; у- показатель адиабаты.

Применение правила SPH (1) к члену правой части уравнения (9) /иУ2у, описывающему вязкость, приводит к формуле, применимой для вычисления асимметричной силы вязкости

где Vу - скоростьу-ой частицы.

Данная сила асимметрична, так как поле скоростей различно для разных частиц.

Так как вязкость зависит от разности скоростей и не от абсолютных скоростей, то существует достаточно простой и естественный способ сделать эту силу симметричной. Достаточно просто использовать разности скоростей вместо абсолютной скорости

где vi - скорость /-ой частицы.

Уравнение (17) можно интерпретировать следующим образом: частица ускоряется в направлении движения окружающей её среды (других частиц).

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

Молекулы жидкости подвержены притяжению со стороны соседних молекул. Внутри жидкости эти межмолекулярные силы одинаковы во всех направлениях и балансируют друг друга. Силы же, воздействующие на молекулы на поверхности жидкости - не сбалансированы. Сила поверхностного натяжения действует по направлению нормали к поверхности жидкости. Также эта сила способствует минимизации искривления поверхности. Чем

(15)

(16)

(17)

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

Поверхность жидкости находится с использованием дополнительного поля, имеющего значения 1 в местах нахождения частиц и 0 в других позициях. Это поле называется «цветовым полем». При применении (1) к данному полю, получаем

( Г ) = Е ту—(г - Г'к )'

(18)

т

/ 'Р.

где с8 - значение сглаженного цветового поля для частицы /. Градиент сглаженного цветового поля

п = Ус3, (19)

даёт поле нормалей, направленных к поверхности жидкости, а дивергенция п измеряет искривление поверхности

-У2сч

к = . (20)

Iп I

Минус в формуле (20) необходим для того, чтобы получить положительную кривизну для выпуклых объемов жидкости.

Сила поверхностного натяжения формулируется следующим образом

ГФсе= акп (21)

I п |

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

¡шфсе= <жп. (22)

Вычисление п / |п| в тех позициях, где |п| мало, может вызывать проблемы с вычислениями. Сила (22) вычисляется только если |п| превышает некоторое заданное минимальное значение.

Внешние силы в уравнении (9) представлены гравитацией pg. Сила гравитации применяется к частицам непосредственно, без использования 8РИ. Аналогичным образом к частицам могут применяться и другие внешние силы. В частности, в разработанной модели существуют твёрдые стенки,

которые, при прохождении частицы через них, отражают зеркально её скорость и выталкивают частицу обратно в моделируемое пространство. Значение ускорения свободного падения берётся равным 9,8 м/с2.

3.2. Ядра сглаживания, используемые в численной модели

Устойчивость, точность и требовательность к ресурсам метода SPH сильно зависит от выбора сглаживающих функций (ядер сглаживания). Используемые в данной работе ядра сглаживания имеют второй порядок точности, так как они удовлетворяют условиям (2) и (3) и (4). В дополнение к этому, ядра такого типа способствуют устойчивости системы [1]. В данной работе используется три различных ядра сглаживания. Первое из них используется во всех случаях, кроме двух, указанных далее и формулируется следующим образом

График этого ядра сглаживания и двух других ядер, используемых в данной работе изображен на рисунке 1. Толстой чертой изображены собственно графики сглаживающих функций, тонкой чертой - их градиенты, и пунктирной - лапласианы.

Важная черта ядра сглаживания (23) состоит в том, что г в нём появляется только в возведённом в квадрат виде, что означает, что отсутствует необходимость вычислять квадратный корень для вычисления скалярного расстояния, что значительно уменьшает сложность вычислений и, соответственно, требования к машинным ресурсам. Тем не менее, при применении данного ядра для подсчёта силового поля давления, частицы часто собираются в отдельные кластеры под высоким давлением. Как только частицы

(h - r j , если 0 < r < h; 0, в другом случае.

(23)

Рис. 1. Графики трёх ядер сглаживания (poly6, spiky, viscosity - слева направо)

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

Данное ядро производит необходимые отталкивающие силы, так как его градиент не стремится к нулю в центре, а следовательно, при его использовании для вычисления сил давления не будет возникать проблем, описанных выше. Также для сглаживающей функции (24) производные первого и второго порядка на границе равны нулю.

Вязкость - это явление, вызываемое трением, и следовательно, вязкость уменьшает кинетическую энергию, преобразовывая её в тепловую. Таким образом, вязкость должна иметь лишь только сглаживающий эффект на поле скоростей и не должна ускорять частицы. Если использовать в качестве ядра сглаживания для вязкости (23) или (24), то они не дадут необходимых качеств полученной в результате вычислений силе. Как видно на рисунке 1, для двух частиц, находящихся рядом друг с другом лапласиан сглаженного поля скоростей (от которого зависит вязкость) может становиться отрицательным, что приводит к вычислению таких сил, которые будут ускорять относительную скорость частицы, вместо того, чтобы уменьшать её. При вычислениях в реальном времени и сравнительно небольшом числе частиц данный артефакт может вызывать проблемы с устойчивостью модели. Таким образом, в данной работе для сил вязкости используется следующее ядро сглаживания

Лапласиан данной функции принимает положительные значения на всём промежутке действия ядра. Также данная функция имеет следующие дополнительные свойства

(25)

0, в другом случае.

VWU ( Гh (h - r ),

Wpy (| Г |= h, h) = 0, VWp&y (| r |= h,h) = 0.

45

spiky

(26)

(27)

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

3.4. Отслеживание и визуализация поверхности жидкости

Цветовое поле (18) и его градиент (19) могут быть использованы для того, чтобы идентифицировать частицы, находящиеся на поверхности жидкости (или газа) и для того, чтобы вычислять нормали к поверхности.

Частица под номером / считается частицей, находящейся на поверхности, если

где I - параметр моделирования, обозначающий некоторую границу между частицами на поверхности жидкости и частицами внутри жидкости.

Направление нормали к поверхности в позиции ьой частицы можно определяется выражением

4. Результаты тестирования программного приложения

Тестовый случай прорыва плотины позволяет проследить воздействие различных силовых полей на частицы и наглядно оценить соответствие модели экспериментальным данным. Данный тестовый случай был взят с сайта организации SPHERIC (SPH European Research Interest Community), где он был описан со ссылкой на статью Имре Яноши [10].

Исходная ситуация описана на схеме, изображенной на рис. 2. Параметр do задаёт глубину воды по левую сторону от плотины и в данном эксперименте равен 15 сантиметрам. Параметр d устанавливает глубину воды по правую сторону от плотины и равен 38 миллиметрам. На расстоянии 38 сантиметров вправо от левой стенки находится барьер (плотина), разделяющий экспериментальную установку на левую часть с глубиной 15 сантиметров и правую с глубиной воды 38 миллиметров. Общая длина экспериментальной установки - 1 метр 4 сантиметра.

I n ( rt )!> I,

(29)

Nsurface (Г, ) = ~П (Г )•

(30)

do

0.38

Рис. 2. Начальное состояние экспериментальной установки

Эксперимент записывался на две камеры одновременно, запись была проведена на участке оси X от 0,38 м до 1,04 м. На рис. 3 изображены результаты эксперимента, полученные Яноши [10].

Рис. 3. Кадры с двух камер, демонстрирующие результаты эксперимента

Результаты проведения эксперимента по прорыву плотины на реализованной численной модели приведены на рисунке 4.

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

Наблюдается качественное соответствие между экспериментом Яноши и экспериментом, проведённым на реализованной численной модели. Добиться количественного соответствия возможно, увеличив количество частиц, используемых в модели, установив движущийся вверх твёрдый барьер (плотину), а так же продолжая подбор параметров модели до получения необходимых результатов.

Рис. 4. Результаты эксперимента на численной модели

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

1. Muller M. Particle-Based Fluid Simulation for Interactive Applications / M. Muller, D. Charypar, M. Gross // Eurographics/SIGGGRAPH Symposium on Computer Animation / Edited by D. Breen and M. Lin. - Zurich: Federal Institute of Technology Zurich (ETHZ), 2003. - p. 7.

2. Schechter H. Ghost SPH for Animating Water / H. Schechter, R. Bridson // Eurographics/SIGGGRAPH Symposium on Computer Animation / Edited by N. Thuerey. - Vancouver: University of British Columbia (UBC), 2012. - p. 8.

3. Monaghan J.J. Simulating Free Surface Flows with SPH / J.J. Monaghan // Journal of Computational Physics. - 1994. - № 110 (2). - p. 399-406.

4. Lucy L.B. A numerical approach to the testing of the fission hypothesis / L.B. Lucy // The Astronomical Journal. - 1977. - № 82. - p. 1013-1024.

5. Desbrun M. Smoothed particles: A new paradigm for animating highly de-formable bodies / M. Desbrun, M.P. Cani // Computer Animation and Simulation

'96 (Proceedings of EG Workshop on Animation and Simulation) / Edited by Edited by D. Thalman and M. van de Panne. - Wien: Springer-Verlag, 1996. - p. 61-76.

6. Braune, L. An initiation to SPH / L. Braune, T. Lewiner - Rio de Janeiro: PUC-Rio, 2013. - p. 1-7.

7. Sun F. Investigations of Boundary Treatments in Incompressible Smoothed Particle Hydrodynamics for Fluid-Structural Interactions / F. Sun, M. Tan, J. Xing // Recent Researches in Mechanics / Edited by N. Mastorakis - Corfu Island: WSEAS Press, 2011. - p. 92-97.

8. Muller M. Particle-Based Fluid Interaction / M. Muller, B. Solenthaler, R. Keiser, M. Gross // Eurographics/SIGGGRAPH Symposium on Computer Animation / Edited by D. Breen and M. Lin. - Zurich: Federal Institute of Technology Zurich (ETHZ), 2005. - p. 1-7.

9. Randles P.W. Smoothed Particle Hydrodynamics: Somer recent improvements and applications / P.W. Randles, L.D. Libersky // Computer Methods in Appliedn Mechanics and Engineering - 1996 - № 139 - p. 375-408.

10. Janosi I.M. Turbulent drag reduction in dam-break flows / I.M. Janosi, D. Jan, K.G. Szabo, T. Tel // Experiments in Fluids - 2004 - № 37 - p. 219-229.

РАЗРАБОТКА POINT-OF-INTEREST РЕКОМЕНДАТЕЛЬНОЙ СИСТЕМЫ ПРИ ИСПОЛЬЗОВАНИИ LBSN СОЦИАЛЬНЫХ СЕТЕЙ

© Пахомова К.И.1

Сибирский федеральный университет, г. Красноярск

В статье рассматривается проблема изучения жизненного цикла человека по средствам извлечения данных POI (Point-of-Interest) из LBSN (Location Based Social Network). На основе алгоритмов классификации и прогнозирования строится рекомендательная система, которая в свою очередь предугадывает поведение пользователя при смене локации в соответствии с его собственными привычками и традициями. На основе алгоритмов кластеризации строится модель рекомендательной системы, включающая временные, хронологические, географические и личные особенности пользователя.

Ключевые слова: POI (Point-of-Interest), LBSN (Location Based Social Network), рекомендательные системы, кластеризация, жизненный цикл человека, численный вероятностный анализ.

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

1 Аспирант кафедры Искусственного интеллекта.