Гидродинамика сглаженных частиц

Схематическое изображение свертки SPH

Обтекание цилиндра со свободной поверхностью, смоделированной с помощью SPH. Видеть для подобных симуляций.

Гидродинамика сглаженных частиц ( SPH ) — это вычислительный метод, используемый для моделирования механики сплошных сред, таких как механика твердого тела и потоков жидкости . Его разработали Джингольд и Монаган. и Люси в 1977 году, первоначально для астрофизических проблем. Его использовали во многих областях исследований, включая астрофизику , баллистику , вулканологию и океанографию . Это бессеточный лагранжев метод (где координаты движутся вместе с жидкостью), и разрешение метода можно легко настроить по отношению к таким переменным, как плотность .

Метод

Преимущества

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

• SPH можно легко распространить на широкий спектр областей и гибридизировать с некоторыми другими моделями, как обсуждается в разделе «Физическое моделирование» .
• Как обсуждалось в разделе о слабосжимаемом SPH , этот метод обладает отличными свойствами сохранения.
• Вычислительная стоимость моделирования SPH на количество частиц значительно меньше, чем стоимость моделирования на основе сетки на количество ячеек, когда интересующая метрика связана с плотностью жидкости (например, функция плотности вероятности флуктуаций плотности). Это так, потому что в SPH резолюция ставится там, где находится вопрос.

Ограничения

• Установка граничных условий в SPH, таких как входы и выходы. и стены сложнее, чем с использованием сеточных методов. Фактически было заявлено, что «обработка граничных условий, безусловно, является одним из самых сложных технических моментов метода SPH». Эта проблема отчасти связана с тем, что в SPH частицы вблизи границы изменяются со временем. Тем не менее, граничные условия для SPH доступны.
• Вычислительная стоимость моделирования SPH на количество частиц значительно превышает стоимость моделирования на основе сетки на количество ячеек, когда интересующая метрика не связана (напрямую) с плотностью (например, спектр кинетической энергии). Поэтому, игнорируя проблемы параллельного ускорения , моделирование потоков постоянной плотности (например, внешней аэродинамики ) более эффективно с использованием сеточных методов, чем с помощью SPH.

Примеры

Гидродинамика

Рис. SPH-моделирование океанских волн с использованием FLUIDS v.1 (Hoetzlein)

Гидродинамика сглаженных частиц все чаще используется для моделирования движения жидкости . Это связано с рядом преимуществ по сравнению с традиционными методами, основанными на сетке. Во-первых, SPH гарантирует сохранение массы без дополнительных вычислений, поскольку частицы сами по себе представляют массу. Во-вторых, SPH вычисляет давление на основе взвешенных вкладов соседних частиц, а не путем решения линейных систем уравнений. Наконец, в отличие от методов на основе сетки, которые должны отслеживать границы жидкости, SPH создает свободную поверхность для двухфазных взаимодействующих жидкостей напрямую, поскольку частицы представляют собой более плотную жидкость (обычно воду), а пустое пространство представляет собой более легкую жидкость (обычно воздух). По этим причинам можно моделировать движение жидкости с помощью SPH в реальном времени. Однако как методы на основе сетки, так и методы SPH по-прежнему требуют создания визуализируемой геометрии свободной поверхности с использованием техники полигонизации, такой как меташары и марширующие кубы , точечные пятна или визуализация «ковра». Для газовой динамики более уместно использовать саму функцию ядра для визуализации плотности столба газа (например, как это сделано в пакете визуализации SPLASH).

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

Недавние работы в области SPH для моделирования жидкости увеличили производительность, точность и области применения:

• Б. Соленталер, 2009 г., разрабатывает прогнозно-корректирующий SPH (PCISPH), позволяющий лучше ограничить несжимаемость.
• М. Имсен и др., 2010 г., представляют обработку границ и адаптивное изменение времени для PCISPH для точного взаимодействия с твердыми телами.
• К. Бодин и др., 2011, заменяют стандартное уравнение давления состояния ограничением плотности и применяют вариационный интегратор времени.
• Р. Хетцляйн, 2012 г., разрабатывает эффективный SPH на базе графического процессора для больших сцен в Fluids v.3.
• Н. Акинчи и др., 2012 г., представляют универсальную обработку границ и метод двустороннего SPH-жесткого соединения, который полностью основан на гидродинамических силах; подход применим к различным типам решателей SPH
• M. Macklin et al., 2013 моделирует несжимаемые потоки в рамках позиционной динамики для больших временных шагов.
• Н. Акинчи и др., 2013, представляют универсальную технику поверхностного натяжения и двусторонней адгезии жидкость-твердое тело, которая позволяет моделировать множество интересных физических эффектов, наблюдаемых в реальности.
• Дж. Кайл и Э. Террелл, 2013 г., применяют SPH для полнопленочной смазки.
• А. Махдави и Н. Талеббейдохти, 2015 г., предлагают гибридный алгоритм для реализации твердого граничного условия и моделирования течения через плотину с острым гребнем.
• С. Таваккол и др., 2016 г., разработали curvSPH, который делает горизонтальный и вертикальный размер частиц независимыми и создает равномерное распределение массы вдоль изогнутых границ.
• В. Косторц и А. Эсмаил-Якас, 2020, предлагают общий, эффективный и простой метод оценки коэффициентов нормализации вблизи кусочно-плоских границ.
• Колагросси и др., 2019, изучают обтекание цилиндра вблизи свободной поверхности и сравнивают его с другими методами.

Астрофизика

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

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

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

Твердая механика

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

Еще одним важным преимуществом бессеточных методов в целом и SPH в частности является то, что проблемы зависимости от сетки естественным образом избегаются, учитывая бессеточную природу метода. В частности, выравнивание сетки связано с проблемами, связанными с трещинами, и его избегают в SPH благодаря изотропной поддержке функций ядра. Однако классические составы SPH страдают нестабильностью растяжения. и отсутствие последовательности. За последние годы для повышения точности решения SPH были введены различные поправки, что привело к RKPM , предложенному Liu et al. Рэндлс и Либерски и Джонсон и Бейссель пытались решить проблему непротиворечивости в своих исследованиях ударных явлений.

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

Недавние улучшения в понимании сходимости и стабильности SPH позволили найти более широкое применение в механике твердого тела. Другие примеры применения и развития метода включают:

• Моделирование обработки металлов давлением.
• Метод SPAM (прикладная механика сглаженных частиц) на основе SPH для ударного разрушения твердых тел Уильяма Г. Гувера .
• Модифицированный SPH (SPH/MLSPH) для лечения переломов и фрагментаций.
• Тейлор-SPH (TSPH) для распространения ударных волн в твердых телах.
• Обобщенная координата SPH (GSPH) распределяет частицы неоднородно в декартовой системе координат и упорядочивает их посредством отображения в обобщенной системе координат, в которой частицы выровнены на одинаковом расстоянии.

Числовые инструменты

Интерполяции

Метод гидродинамики сглаженных частиц (SPH) работает путем разделения жидкости на набор дискретных движущихся элементов. $i,j$, называемые частицами. Их лагранжева природа позволяет установить их положение $\mathbf {r} _{i}$ путем интегрирования их скорости $\mathbf {v} _{i}$ как:

${\frac {\mathrm {d} {\boldsymbol {r}}_{i}}{\mathrm {d} t}}={\boldsymbol {v}}_{i}.$

Эти частицы взаимодействуют через ядерную функцию с характерным радиусом, известным как «длина сглаживания», обычно представляемую в уравнениях как $h$. Это означает, что физическая величина любой частицы может быть получена путем суммирования соответствующих свойств всех частиц, находящихся в пределах диапазона ядра, причем последнее используется в качестве весовой функции. $W$. Это можно понять в два этапа. Сначала произвольное поле $A$ записывается как свертка с $W$:

$A({\boldsymbol {r}})=\int A\left({\boldsymbol {r^{\prime }}}\right)W(|{\boldsymbol {r}}-{\boldsymbol {r^{\prime }}}|,h)\,\mathrm {d} V\!\left({\boldsymbol {r'}}\right).$

Ошибка в приведенном выше приближении составляет порядок $h^{2}$. Во-вторых, интеграл аппроксимируется суммированием Римана по частицам:

$A({\boldsymbol {r}})=\sum _{j}V_{j}A_{j}W(|{\boldsymbol {r}}-{\boldsymbol {r}}_{j}|,h),$

где суммирование по $j$ включает все частицы в моделировании. $V_{j}$ это объем частицы $j$, $A_{j}$ это значение количества $A$ для частицы $j$ и ${\boldsymbol {r}}$ обозначает позицию. Например, плотность $\rho _{i}$ частицы $i$ может быть выражено как:

$\rho _{i}=\rho ({\boldsymbol {r}}_{i})=\sum _{j}m_{j}W_{ij},$

где $m_{j}=\rho _{j}V_{j}$ обозначает массу частицы и $\rho _{j}$ плотность частиц, в то время как $W_{ij}=W_{ji}$ это краткое обозначение $W(|{\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{j}|,h)$. Ошибка приближения интеграла дискретной суммой зависит от $h$, от размера частиц (т.е. $V_{j}^{1/d}$, $d$ будучи измерением пространства), а также от расположения частиц в пространстве. Последний эффект пока малоизучен.

Обычно используемые функции ядра включают функцию Гаусса , сплайн пятой степени и функцию Вендланда. $C^{2}$ ядро. Последние два ядра имеют компактный носитель (в отличие от гауссовского, где вклад небольшой на любом конечном расстоянии) с поддержкой, пропорциональной $h$. Преимущество этого подхода заключается в экономии вычислительных усилий за счет исключения относительно незначительного вклада удаленных частиц.

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

Дискретизация основных уравнений

Для частиц постоянной массы, дифференцируя интерполированную плотность $\rho _{i}$ относительно временной доходности

${\frac {d\rho _{i}}{dt}}=\sum _{j}m_{j}\left({\boldsymbol {v}}_{i}-{\boldsymbol {v}}_{j}\right)\cdot \nabla W_{ij},$

где $\nabla W_{ij}=-\nabla W_{ji}$ это градиент $W_{ij}$ относительно ${\boldsymbol {r}}_{i}$. Сравнивая это уравнение с уравнением неразрывности в лагранжевом описании (с использованием материальных производных ),

${\frac {d\rho }{dt}}=-\rho \nabla \cdot {\boldsymbol {v}},$

очевидно, что его правая часть является приближением $-\rho \nabla \cdot \mathbf {v}$; следовательно, дискретный оператор дивергенции определяется следующим образом:

$\operatorname {D} _{i}\left\{{\boldsymbol {v}}_{j}\right\}=-{\frac {1}{\rho _{i}}}\sum _{j}m_{j}\left({\boldsymbol {v}}_{i}-{\boldsymbol {v}}_{j}\right)\cdot \nabla W_{ij}.$

Этот оператор дает SPH-аппроксимацию $\nabla \cdot \mathbf {v}$ на частице $i$ для данного набора частиц с заданными массами $m_{j}$, позиции $\left\{\mathbf {r} _{j}\right\}$ и скорости $\left\{\mathbf {v} _{j}\right\}$.

Другим важным уравнением для сжимаемой невязкой жидкости является уравнение Эйлера для баланса импульса:

${\frac {d{\boldsymbol {v}}}{dt}}=-{\frac {1}{\rho }}\nabla p+{\boldsymbol {g}}$

Как и в случае с непрерывностью, задача состоит в том, чтобы определить дискретный оператор градиента, чтобы записать

${\frac {d{\boldsymbol {v}}_{i}}{dt}}=-{\frac {1}{\rho }}\operatorname {\mathbf {G} } _{i}\left\{p_{j}\right\}+{\boldsymbol {g}}$

Один из вариантов

$\operatorname {\mathbf {G} } _{i}\left\{p_{j}\right\}=\rho _{i}\sum _{j}m_{j}\left({\frac {p_{i}}{\rho _{i}^{2}}}+{\frac {p_{j}}{\rho _{j}^{2}}}\right)\nabla W_{ij},$

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

$\sum _{i}V_{i}{\boldsymbol {v}}_{i}\cdot \operatorname {\mathbf {G} } _{i}\left\{p_{j}\right\}=-\sum _{i}V_{i}p_{i}\operatorname {D} _{i}\left\{{\boldsymbol {v}}_{j}\right\},$

это дискретная версия континуальной идентичности

$\int {\boldsymbol {v}}\cdot \operatorname {grad} p=-\int p\operatorname {div} \cdot {\boldsymbol {v}}.$

Это свойство приводит к хорошим консервационным свойствам.

Обратите также внимание, что этот выбор приводит к симметричному оператору дивергенции и антисимметричному градиенту. Хотя существует несколько способов дискретизации градиента давления в уравнениях Эйлера, приведенная выше антисимметричная форма является наиболее признанной. Он поддерживает строгое сохранение линейного и углового момента. Это означает, что сила, действующая на частицу $i$ по частице $j$ равна той, которая действует на частицу $j$ по частице $i$ включая смену знака эффективного направления, благодаря свойству антисимметрии $\nabla W_{ij}=-\nabla W_{ji}$.

Тем не менее, были предложены другие операторы, которые могут работать лучше численно или физически. Например, одним из недостатков этих операторов является то, что, хотя расхождение $\operatorname {D}$ является согласованным нулевого порядка (т.е. дает ноль при применении к постоянному векторному полю), можно видеть, что градиент $\operatorname {\mathbf {G} }$ нет. Чтобы обойти эту проблему, было предложено несколько методов, ведущих к перенормированным операторам (см., например, ).

Вариационный принцип

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

${\mathcal {L}}=\sum _{j}m_{j}\left({\tfrac {1}{2}}{\boldsymbol {v}}_{j}^{2}-e_{j}+{\boldsymbol {g}}\cdot {\boldsymbol {r}}_{j}\right)$,

где $e_{j}$ частицы – удельная внутренняя энергия . Уравнение Эйлера-Лагранжа вариационной механики для каждой частицы гласит:

${\frac {\mathrm {d} }{\mathrm {d} t}}{\frac {\partial {\mathcal {L}}}{\partial {\boldsymbol {v}}_{i}}}={\frac {\partial {\mathcal {L}}}{\partial {\boldsymbol {r}}_{i}}}.$

Применительно к приведенному выше лагранжиану он дает следующее уравнение импульса:

$m_{i}{\frac {\mathrm {d} {\boldsymbol {v}}_{i}}{\mathrm {d} t}}=-\sum _{j}m_{j}{\frac {\partial e_{j}}{\partial {\boldsymbol {r}}_{i}}}+m_{i}{\boldsymbol {g}}=-\sum _{j}m_{j}{\frac {\partial e_{j}}{\partial \rho _{j}}}{\frac {\partial \rho _{j}}{\partial {\boldsymbol {r}}_{i}}}+m_{i}{\boldsymbol {g}}$

где использовалось правило цепочки, поскольку $e_{j}$ зависит от $\rho _{j}$, а последнее – от положения частиц. Использование термодинамического свойства $\mathrm {d} e=\left(p/\rho ^{2}\right)\mathrm {d} \rho$ мы можем написать

$m_{i}{\frac {\mathrm {d} {\boldsymbol {v}}_{i}}{\mathrm {d} t}}=-\sum _{j}m_{j}{\frac {p_{j}}{\rho _{j}^{2}}}{\frac {\partial \rho _{j}}{\partial {\boldsymbol {r}}_{i}}}+m_{i}{\boldsymbol {g}},$

Подключение интерполяции плотности SPH и явное дифференцирование ${\tfrac {\partial \rho _{j}}{\partial {\boldsymbol {r}}_{i}}}$ приводит к

${\frac {\mathrm {d} {\boldsymbol {v}}_{i}}{\mathrm {d} t}}=-\sum _{j}m_{j}\left({\frac {p_{i}}{\rho _{i}^{2}}}+{\frac {p_{j}}{\rho _{j}^{2}}}\right)\nabla W_{ij}+{\boldsymbol {g}},$

которое представляет собой уже упомянутое уравнение импульса SPH, где мы признаем $\operatorname {\mathbf {G} }$ оператор. Это объясняет, почему сохраняется линейный момент, а также позволяет сохранить угловой момент и энергию.

Интеграция времени

На основе работ, проделанных в 80-х и 90-х годах по численному интегрированию точечных частиц в больших ускорителях, были разработаны подходящие интеграторы времени с точными свойствами сохранения в долгосрочной перспективе; их называют симплектическими интеграторами . Наиболее популярной в литературе по SPH является схема чехарды , при которой для каждой частицы считывается $i$:

${\begin{aligned}{\boldsymbol {v}}_{i}^{n+1/2}&={\boldsymbol {v}}_{i}^{n}+{\boldsymbol {a}}_{i}^{n}{\frac {\Delta t}{2}},\\{\boldsymbol {r}}_{i}^{n+1}&={\boldsymbol {r}}_{i}^{n}+{\boldsymbol {v}}_{i}^{i+1/2}\Delta t,\\{\boldsymbol {v}}_{i}^{n+1}&={\boldsymbol {v}}_{i}^{n+1/2}+{\boldsymbol {a}}_{i}^{i+1}{\frac {\Delta t}{2}},\end{aligned}}$

где $\Delta t$ — это шаг по времени, верхние индексы обозначают итерации по времени, а ${\boldsymbol {a}}_{i}$ – ускорение частицы, определяемое правой частью уравнения количества движения.

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

Интегрирование плотности не изучалось подробно ( см. Ниже подробнее ).

Симплектические схемы консервативны, но явны, поэтому для их численной устойчивости требуются условия устойчивости, аналогичные условию Куранта-Фридрихса-Леви (см. ниже ).

Граничные методы

Поддержка SPH Convolution раскололась вблизи границы

В случае, если свертка SPH должна выполняться близко к границе, т.е. ближе, чем s · h , тогда интегральная поддержка усекается. Действительно, когда на свертку влияет граница, свертка должна быть разделена на два интеграла:

$A({\boldsymbol {r}})=\int _{\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)W(|{\boldsymbol {r}}-{\boldsymbol {r^{\prime }}}|,h)d{\boldsymbol {r^{\prime }}}+\int _{B({\boldsymbol {r}})-\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)W(|{\boldsymbol {r}}-{\boldsymbol {r^{\prime }}}|,h)d{\boldsymbol {r^{\prime }}},$

где B( r ) — компактный опорный шар с центром в r и радиусом s · h , а Ω( r ) обозначает часть компактного носителя внутри расчетной области, Ω ∩ B( r ) . Следовательно, наложение граничных условий в SPH полностью основано на аппроксимации второго интеграла в правой части. То же самое, конечно, можно применить и к вычислению дифференциальных операторов:

$\nabla A({\boldsymbol {r}})=\int _{\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)\nabla W({\boldsymbol {r}}-{\boldsymbol {r^{\prime }}},h)d{\boldsymbol {r^{\prime }}}+\int _{B({\boldsymbol {r}})-\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)\nabla W({\boldsymbol {r}}-{\boldsymbol {r^{\prime }}},h)d{\boldsymbol {r^{\prime }}}.$

В прошлом было предложено несколько методов моделирования границ в SPH.

Комплексное пренебрежение

Модель свободной поверхности SPH посредством интегрального пренебрежения

Самая простая граничная модель пренебрегает интегралом:

$\int _{B({\boldsymbol {r}})-\Omega ({\boldsymbol {r}})}A\left({\boldsymbol {r^{\prime }}}\right)\nabla W({\boldsymbol {r}}-{\boldsymbol {r^{\prime }}},h)d{\boldsymbol {r^{\prime }}}\simeq {\boldsymbol {0}},$

такие, что учитываются только объемные взаимодействия,

$\nabla A_{i}=\sum _{j\in \Omega _{i}}V_{j}A_{j}\nabla W_{ij}.$

Это популярный подход, когда свободная поверхность рассматривается в однофазном моделировании.

Основным преимуществом этого граничного условия является его очевидная простота. Однако при применении этого метода границ необходимо учитывать несколько проблем согласованности. На самом деле это серьезное ограничение на его потенциальное применение.

Расширение жидкости

Техника расширения границ жидкости SPH

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

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

Основным преимуществом этой методологии является простота при условии, что вклад границ вычисляется как часть объемных взаимодействий. Кроме того, эта методология была глубоко проанализирована в литературе.

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

Граничный интеграл

Интегральная модель границы SPH

Новейшей граничной методикой является методология граничного интеграла. В этой методологии интеграл пустого объема заменяется поверхностным интегралом и перенормировкой:

$\nabla A_{i}={\frac {1}{\gamma _{i}}}\left(\sum _{j\in \Omega _{i}}V_{j}A_{j}\nabla W_{ij}+\sum _{j\in \partial \Omega _{i}}S_{j}A_{j}{\boldsymbol {n}}_{j}W_{ij}\right),$

$\gamma _{i}=\sum _{j\in \Omega _{i}}V_{j}W_{ij},$

где n _{j —} нормаль общего j -го граничного элемента. Поверхностный член также можно решить, используя полуаналитическое выражение.

Моделирование физики

Гидродинамика

Слабо сжимаемый подход

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

Слабо сжимаемый SPH в гидродинамике основан на дискретизации уравнений Навье – Стокса или уравнений Эйлера для сжимаемых жидкостей. Чтобы закрыть систему, соответствующее уравнение состояния , связывающее давление используется $p$ и плотность $\rho$. В общем случае так называемое уравнение Коула (иногда ошибочно называемое « уравнением Тейта ») используется в SPH. Он читает

$p={\frac {\rho _{0}c^{2}}{\gamma }}\left(\left({\frac {\rho }{\rho _{0}}}\right)^{\gamma }-1\right)+p_{0},$

где $\rho _{0}$ - эталонная плотность и $c$ скорость звука . Для воды, $\gamma =7$ обычно используется. Фоновое давление $p_{0}$ добавляется, чтобы избежать отрицательных значений давления.

Реальные почти несжимаемые жидкости, такие как вода, характеризуются очень высокими скоростями звука порядка $10^{3}\mathrm {m/s}$. Следовательно, информация о давлении распространяется быстрее по сравнению с фактическим объемным потоком, что приводит к очень малым числам Маха. $M$. Уравнение количества движения приводит к следующему соотношению:

${\frac {\delta \rho }{\rho _{0}}}\approx {\frac {|{\boldsymbol {v}}|^{2}}{c^{2}}}=M^{2}$

где $\rho$ - изменение плотности и $v$ вектор скорости. На практике значение c меньше фактического принимается, чтобы избежать слишком малых временных шагов в схеме интегрирования по времени. Обычно числовая скорость звука принимается такой, что допускается изменение плотности менее 1%. Это так называемое предположение о слабой сжимаемости. Это соответствует числу Маха меньше 0,1, что означает:

$c=10v_{\text{max}}$

где максимальная скорость $v_{\text{max}}$ необходимо оценить, например, с помощью закона Торричелли или обоснованного предположения. Поскольку происходят лишь небольшие изменения плотности, можно принять линейное уравнение состояния:

$p=c^{2}\left(\rho -\rho _{0}\right)$

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

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

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

Метод фильтра плотности

В схемах первой группы фильтр применяется непосредственно к полю плотности для удаления паразитного числового шума. Наиболее часто используемые фильтры — это MLS (по методу скользящих наименьших квадратов) и фильтр Шепарда. который может применяться на каждом временном шаге или на каждом n временных шагах. Чем чаще используется процедура фильтрации, тем более регулярные поля плотности и давления получаются. С другой стороны, это приводит к увеличению вычислительных затрат. При длительном моделировании использование процедуры фильтрации может привести к нарушению гидростатической составляющей давления и к несоответствию между глобальным объемом жидкости и полем плотности. Кроме того, это не обеспечивает соблюдение динамического граничного условия для свободной поверхности.

Диффузионная техника терминов

Другой способ сгладить поле плотности и давления — добавить диффузионный член внутрь уравнения неразрывности (группа 2):

${\displaystyle {\frac {d\rho _{i}}{dt}}=\sum _{j}m_{j}\left({\boldsymbol {v}}_{i}-{\boldsymbol {v}}_{j}\right)\cdot \nabla W_{ij}+{\mathcal {D}}_{i}(\rho ),}$

Первые схемы, в которых применялся такой подход, были описаны в Ferrari. и в Молтени где диффузионный член моделировался как лапласиан поля плотности. Подобный подход также использовался в делах Фатехи и Манзари. .

Моделирование SPH: распределение давления потока при прорыве плотины с использованием стандартной формулы SPH.

Моделирование SPH: распределение давления потока при прорыве плотины с использованием стандартной формулы δ-SPH.

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

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

Техника решения Римана

Моделирование SPH: распределение давления потока при прорыве плотины с использованием решателя Римана с малодиссипативным ограничителем.

Для метода SPH, основанного на решателях Римана, межчастичная задача Римана строится вдоль единичного вектора $\mathbf {e} _{ij}=-\mathbf {r} _{ij}/r_{ij}$ частица указывающей формы $i$ частица $j$. В этой задаче Римана начальные левые и правые состояния находятся на частицах $i$ и $j$ , соответственно. $L$ и $R$ государства являются

$${\begin{cases}(\rho _{L},U_{L},P_{L})=(\rho _{i},\mathbf {v} _{i}\cdot \mathbf {e} _{ij},P_{i})\\(\rho _{R},U_{R},P_{R})=(\rho _{j},\mathbf {v} _{j}\cdot \mathbf {e} _{ij},P_{j}).\end{cases}}$$

Решение задачи Римана приводит к появлению трех волн, исходящих из разрыва. Две волны, которые могут быть ударными или волнами разрежения, распространяющимися с наименьшей или наибольшей волновой скоростью. Средняя волна всегда представляет собой контактный разрыв и разделяет два промежуточных состояния, обозначаемых $(\rho _{L}^{\ast },U_{L}^{\ast },P_{L}^{\ast })$ и $(\rho _{R}^{\ast },U_{R}^{\ast },P_{R}^{\ast })$. Предполагая, что промежуточное состояние удовлетворяет $U_{L}^{\ast }=U_{R}^{\ast }=U^{\ast }$ и $P_{L}^{\ast }=P_{R}^{\ast }=P^{\ast }$, линеаризованный решатель Римана для гладких потоков или только с умеренно сильными толчками можно записать как

$${\begin{cases}U^{\ast }={\overline {U}}+{\frac {1}{2}}{\frac {(P_{L}-P_{R})}{{\bar {\rho }}c_{0}}}\\P^{\ast }={\overline {P}}+{\frac {1}{2}}{\bar {\rho }}c_{0}{(U_{L}-U_{R})},\end{cases}}$$

где ${\overline {U}}=(U_{L}+U_{R})/2$ и ${\overline {P}}=(P_{L}+P_{R})/2$ являются средними значениями между частицами. При решении задачи Римана, т.е. $U^{\ast }$ и $P^{\ast }$, дискретизация метода SPH равна

$${\frac {d\rho _{i}}{dt}}=2\rho _{i}\sum _{j}{\frac {m_{j}}{\rho _{j}}}(\mathbf {v} _{i}-\mathbf {v} ^{\ast })\cdot \nabla _{i}W_{ij},$$

$${\frac {d\mathbf {v} _{i}}{dt}}=-2\sum _{j}m_{j}\left({\frac {P^{\ast }}{\rho _{i}\rho _{j}}}\right)\nabla _{i}W_{ij}.$$

где $\mathbf {v} ^{\ast }=U^{\ast }\mathbf {e} _{ij}+({\overline {\mathbf {v} }}_{ij}-{\overline {U}}\mathbf {e} _{ij})$. Это указывает на то, что средние скорость и давление между частицами просто заменяются решением задачи Римана. Сравнивая оба значения, можно увидеть, что промежуточная скорость и давление из средних значений между частицами представляют собой неявную диссипацию, т.е. регуляризацию плотности и численную вязкость соответственно.

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

$$P^{\ast }={\overline {P}}+{\frac {1}{2}}\beta {\overline {\rho }}{(U_{L}-U_{R})},$$

где ограничитель определяется как

$$\beta =\min {\big (}\eta \max(U_{L}-U_{R},0),{\overline {c}}{\big )}.$$

Обратите внимание, что $\beta$ обеспечивает отсутствие диссипации при нахождении жидкости под действием волны расширения, т.е. $U_{L}<U_{R}$, и что параметр $\eta$, используется для модуляции диссипации, когда жидкость находится под действием волны сжатия, т.е. $U_{L}\geq U_{R}$. Численные эксперименты выявили $\eta =3$ в целом эффективен. Также обратите внимание, что диссипация, вносимая промежуточной скоростью, не ограничена.

Моделирование вязкости

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

• Искусственная вязкость

Представлено Монаганом и Джингольдом искусственная вязкость использовалась для борьбы с потоками жидкости с высоким числом Маха . Он читает

$\Pi _{ij}={\begin{cases}{\dfrac {-\alpha {\bar {c}}_{ij}\phi _{ij}+\beta \phi _{ij}^{2}}{{\bar {\rho }}_{ij}}}&\quad {\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}<0\\0&\quad {\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}\geq 0\end{cases}}$

Здесь, $\alpha$ контролирует объемную вязкость, в то время как $\beta$ действует аналогично искусственной вязкости Неймана-Рихтмейра. $\phi _{ij}$ определяется

$\phi _{ij}={\frac {h{\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}}{\Vert {\boldsymbol {r}}_{ij}\Vert ^{2}+\eta _{h}^{2}}},$

где η _h — небольшая часть h (например, 0,01 h ), чтобы предотвратить возможную числовую бесконечность на близких расстояниях.

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

$\Pi _{ij}=\alpha hc{\frac {{\boldsymbol {v}}_{ij}\cdot {\boldsymbol {r}}_{ij}}{\Vert {\boldsymbol {r}}_{ij}\Vert ^{2}+\eta _{h}^{2}}}.$

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

$\alpha hc=2(n+2){\frac {\mu }{\rho }}$

подставляется в приведенное выше уравнение, где $n$ — количество пространственных измерений модели. Этот подход вводит объемную вязкость $\zeta ={\frac {5}{3}}\mu$.

• Моррис

Для низких чисел Рейнольдса модель вязкости Морриса было предложено.

$[\nu \Delta {\boldsymbol {v}}]_{ij}=2\nu {\frac {m_{j}}{\rho _{j}}}\,{\frac {{\boldsymbol {r}}_{ij}\cdot \nabla w_{h,ij}}{\Vert {\boldsymbol {r}}_{ij}\Vert ^{2}+\eta _{h}^{2}}}\,{\boldsymbol {v}}_{ij}.$

• ЛоШао

Дополнительная физика

• Поверхностное натяжение
• Теплопередача
• Турбулентность

Астрофизика

Часто в астрофизике помимо чистой гидродинамики хотят смоделировать самогравитацию. Основанный на частицах характер SPH делает его идеальным в сочетании с гравитационным решателем на основе частиц, например, с кодом гравитации дерева . сетка частиц или частица-частица частица-сетка .

Механика твердого тела и взаимодействие жидкости со структурой (FSI)

Полная лагранжева формулировка для механики твердого тела

Чтобы дискретизировать основные уравнения динамики твердого тела, используется корректирующая матрица $\mathbb {B} ^{0}$ впервые представлено воспроизведение вращения твердого тела как

$\mathbb {B} _{a}^{0}=\left(\sum _{b}V_{b}^{0}\left(\mathbf {r} _{b}^{0}-\mathbf {r} _{a}^{0}\right)\otimes \nabla _{a}^{0}W_{ab}\right)^{-1}$ ( 1 )

где

$$\nabla _{a}^{0}W_{a}={\frac {\partial W\left(|\mathbf {r} _{ab}^{0}|,h\right)}{\partial |\mathbf {r} _{ab}^{0}|}}\mathbf {e} _{ab}^{0}$$

обозначает градиент функции ядра, оцененной в исходной эталонной конфигурации. Обратите внимание, что индексы $a$ и $b$ используются для обозначения твердых частиц и длины сглаживания $h$ идентично таковому при дискретизации уравнений жидкости.

Используя исходную конфигурацию в качестве эталона, плотность твердого тела напрямую оценивается как

$\rho _{a}=\rho _{a}^{0}{\frac {1}{J}}$ ( 2 )

где $J=\det(\mathbb {F} )$ – определитель Якобиана тензора деформации $\mathbb {F}$.

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

$m_{a}{\frac {{\text{d}}\mathbf {v} }{{\text{d}}t}}=2\sum _{b}V_{a}V_{b}{\tilde {\mathbb {P} }}_{ab}\nabla _{a}^{0}W_{ab}+\mathbf {g} +\mathbf {f} _{a}^{F:p}+\mathbf {f} _{a}^{F:v}$ ( 3 )

где межчастичное среднее первое напряжение Пиолы-Кирхгофа ${\tilde {\mathbb {P} }}$ определяется как

${\tilde {\mathbb {P} }}_{ab}={\frac {1}{2}}\left(\mathbb {P} _{a}\mathbb {B} _{a}^{0}+\mathbb {P} _{b}\mathbb {B} _{b}^{0}\right).$ ( 4 )

Также $\mathbf {f} _{a}^{F:p}$ и $\mathbf {f} _{a}^{F:v}$ соответствуют давлению жидкости и силам вязкости, действующим на твердую частицу $a$, соответственно.

Муфта гидроконструкции

При соединении жидкость-структура окружающая твердая структура ведет себя как движущаяся граница для жидкости, а на границе раздела жидкость-структура накладывается граничное условие прилипания. Силы взаимодействия $\mathbf {f} _{i}^{S:p}$ и $\mathbf {f} _{i}^{S:v}$ действие на частицу жидкости $i$, из-за присутствия соседней твердой частицы $a$, можно получить как

$\mathbf {f} _{i}^{S:p}=-2\sum _{a}V_{i}V_{a}{\frac {p_{i}\rho _{a}^{d}+p_{a}^{d}\rho _{i}}{\rho _{i}+\rho _{a}^{d}}}\nabla _{i}W(\mathbf {r} _{ia},h)$ ( 5 )

и

$\mathbf {f} _{i}^{S:v}=2\sum _{a}\eta V_{i}V_{a}{\frac {\mathbf {v} _{i}-\mathbf {v} _{a}^{d}}{r_{ia}}}{\frac {\partial W(\mathbf {r} _{ia},h)}{\partial {{r}_{ia}}}}.$ ( 6 )

Здесь мнимое давление $p_{a}^{d}$ и скорость $\mathbf {v} _{a}^{d}$ определяются

${\begin{cases}p_{a}^{d}=p_{i}+\rho _{i}\max(0,(\mathbf {g} -{\frac {{\text{d}}\mathbf {v} _{a}}{{\text{d}}t}})\cdot \mathbf {n} ^{S})(\mathbf {r} _{ia}\cdot \mathbf {n} ^{S})\\\mathbf {v} _{a}^{d}=2\mathbf {v} _{i}-\mathbf {v} _{a}.\end{cases}}$ ( 7 )

где $\mathbf {n} ^{S}$ обозначает направление нормали к поверхности твердой конструкции, и мнимая плотность частиц $\rho _{a}^{d}$ рассчитывается через уравнение состояния.

Соответственно, силы взаимодействия $\mathbf {f} _{a}^{F:p}$ и $\mathbf {f} _{a}^{F:v}$ действие на твердую частицу $a$ даны

$\mathbf {f} _{a}^{F:p}=-2\sum _{i}V_{a}V_{i}{\frac {p_{a}^{d}\rho _{i}+p_{i}\rho _{a}^{d}}{\rho _{i}+\rho _{a}^{d}}}\nabla _{a}W(\mathbf {r} _{ai},h)$ ( 8 )

и

$\mathbf {f} _{a}^{F:v}=2\sum _{i}\eta V_{a}V_{i}{\frac {\mathbf {v} _{a}^{d}-\mathbf {v} _{i}}{r_{ia}}}{\frac {\partial W(\mathbf {r} _{ia},h)}{\partial {{r}_{ai}}}}.$ ( 9 )

Антисимметричность производной ядерной функции обеспечит сохранение импульса для каждой пары взаимодействующих частиц. $i$ и $a$.

Другие

Метод дискретных элементов , используемый для моделирования сыпучих материалов , относится к SPH.

Дальнейшее чтение

• Гувер, РГ (2006 г.); Прикладная механика гладких частиц: современное состояние, World Scientific.
• Стеллингверф, РФ; Вингейт, Калифорния; «Моделирование воздействия с помощью SPH» , Мемуары Итальянского астрономического общества, Том 65, стр. 1117 (1994).
• Амада, Т.; Имура, М.; Ясумуро, Ю.; Манабе, Ю.; и Чихара, К. (2004); «Моделирование жидкости на основе частиц на графическом процессоре», в Трудах семинара ACM по универсальным вычислениям на графических процессорах (август 2004 г., Лос-Анджелес, Калифорния).
• Дебрун, М.; и Кани, М.-П. (1996). «Сглаженные частицы: новая парадигма анимации сильно деформируемых тел» в Proceedings of Eurographics Workshop on Computer Animation and Simulation (август 1996 г., Пуатье, Франция).
• Хегеман, К.; Карр, Северная Каролина; и Миллер, GSP; «Моделирование жидкости на основе частиц на графическом процессоре», в материалах Международной конференции по вычислительной науке (Ридинг, Великобритания, май 2006 г.), Конспекты лекций по информатике, версия 3994/2006 (Springer-Verlag).
• Келагер, М. (2006) Лагранжева гидродинамика с использованием гидродинамики сглаженных частиц , магистерская диссертация, Univ. Копенгаген.
• Колб, А.; и Кунц Н. (2005); «Динамическое взаимодействие частиц для моделирования жидкости на основе графического процессора», в материалах 18-го симпозиума по методам моделирования (2005), стр. 722–727.
• Лю, гр.; и Лю, МБ; Гидродинамика сглаженных частиц: метод бессеточных частиц , Сингапур: World Scientific (2003).
• Монаган, Джозеф Дж. (1992). «Гидродинамика сглаженных частиц», Ежегодный обзор астрономии и астрофизики (1992). 30: 543–74.
• Мюллер, М.; Чарыпар, Д.; и Гросс, М.; «Моделирование жидкости на основе частиц для интерактивных приложений», Брин, Д.; и Лин, М. (ред.), Труды симпозиума Eurographics/SIGGRAPH по компьютерной анимации (2003).
• Вестерлунд, М.; Моделирование и визуализация вязкой жидкости с использованием гидродинамики сглаженных частиц , магистерская диссертация, Университет Умео, Швеция.
• Виоло, Д.; Механика жидкости и метод SPH , Oxford University Press (2012).

Внешние ссылки

• Первое крупное моделирование звездообразования с использованием SPH.
• SPHERIC (Международное сообщество исследований и инженерии SPH)
• ITVO — это веб-сайт Итальянской теоретической виртуальной обсерватории, созданный для запроса базы данных архива численного моделирования.
• Галерея изображений SPHC отображает широкий спектр тестовых случаев, экспериментальных проверок и коммерческих применений кода SPH SPHC.
• Вывод модели SPH на основе уравнений Навье-Стокса.

Программное обеспечение

• Algodoo — это среда 2D-моделирования для образования с использованием SPH.
• AQUAgpusph — это бесплатная (GPLv3) SPH для исследователей, созданная исследователями и для исследователей.
• Dive Solutions — это коммерческое веб-программное обеспечение SPH для проектирования CFD.
• DualSPhysics — это код SPH с открытым исходным кодом, основанный на SPHysics и использующий вычисления на графическом процессоре. Компоненты с открытым исходным кодом доступны по лицензии LGPL.
• FLUIDS v.1 — это простая реализация 3D SPH в реальном времени с открытым исходным кодом (Zlib) на C++ для жидкостей для ЦП и ГП.
• Fluidix  — это API моделирования частиц на базе графического процессора, доступный в OneZero Software.
• GADGET [1] — это свободно доступный ( GPL ) код для космологического моделирования N-тел/SPH.
• Симулятор GPUSPH SPH с вязкостью (GPLv3)
• Pasimodo — это программный пакет для методов моделирования на основе частиц, например SPH.
• LAMMPS — это массовый параллельный код классической молекулярной динамики с открытым исходным кодом, который может выполнять SPH-моделирование.
• Уровень абстракции физики — это система абстракции с открытым исходным кодом, которая поддерживает физические движки реального времени с поддержкой SPH.
• PreonLab — это коммерческое инженерное программное обеспечение, разработанное FIFTY2 Technology, реализующее неявный метод SPH.
• Punto — это бесплатно доступный инструмент визуализации для моделирования частиц.
• pysph Платформа с открытым исходным кодом для гидродинамики сглаженных частиц в Python (новая лицензия BSD)
• Py-SPHViewer Инструмент визуализации Python с открытым исходным кодом для моделирования гидродинамики сглаженных частиц.
• RealFlow Commercial SPH решатель для киноиндустрии.
• RheoCube — это коммерческий SaaS- продукт от Lorenz Research, предназначенный для изучения и прогнозирования реологии и стабильности сложных жидкостей.
• SimPARTIX — это коммерческий пакет моделирования для моделирования SPH и метода дискретных элементов (DEM) от Fraunhofer IWM.
• SPH-поток
• СФЕРА
• SPHinXsys — это мультифизическая библиотека SPH с открытым исходным кодом и различными разрешениями. Он предоставляет API-интерфейсы C++ для точного физического моделирования и предназначен для моделирования связанных промышленных динамических систем, включая динамику жидкости, твердого тела, многих тел и многое другое.
• SPHysics — это реализация SPH с открытым исходным кодом на Фортране.
• SPLASH — это инструмент визуализации с открытым исходным кодом (GPL) для моделирования SPH.
• SYMPLER : бесплатный симулятор символических частиц от Фрайбургского университета.
• Nauticle — это универсальный вычислительный инструмент для численных методов, основанных на частицах.
• NDYNAMICS — это коммерческое программное обеспечение для моделирования жидкости, основанное на неявном SPH, разработанное CENTROID LAB, которое в настоящее время используется для внутреннего/внешнего заводнения, ядерной/химической инженерии.