Некоторые результаты для различных методов моделирования несжимаемой гидродинамики свободной поверхностью на графических процессорах Текст научной статьи по специальности «Физика»

УДК 532

НЕКОТОРЫЕ РЕЗУЛЬТАТЫ ДЛЯ РАЗЛИЧНЫХ МЕТОДОВ МОДЕЛИРОВАНИЯ НЕСЖИМАЕМОЙ ГИДРОДИНАМИКИ СВОБОДНОЙ ПОВЕРХНОСТЬЮ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ *

© 2010 И. В. Гугушвили1, Н. М. Евстигнеев2

1 аспирант всероссийского НИИ гидротехники и мелиорации e-mail: gugushviliagmail. com

ст. науч. сотрудник Института системного анализа РАН

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

Ключевые слова: уравнение Навье-Стокса, Сеточный метод Больцмана,

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

Введение

Для решения начально-краевых задач вычислительной гидродинамики со свободной поверхностью одними из важнейших факторов алгоритма являются точность отслеживания положения фазы раздела сред (свободной поверхности), вычислительная эффективность и скорость расчета. Кроме того, к методам решения предъявляются особые требования, связанные с выполнением законов сохранения. Численная схема должна отображать фазу раздела четко, не смазывая ее при длительном интегрировании: алгоритм должен учитывать возможность образования невыпуклых областей свободной поверхности, и на свободной поверхности для двух жидкостей в общем случае должно выполняться условие:

dF

(P1 - P2 + FS ' k) ' Пг = (Tlik - T2,k ) ' Пк + 7FL ,

dxi

где 1, 2 - номера жидкостей; P - давление, n - вектор внешней нормали к фазе раздела; Fs - поверхностное натяжение (Н/м); x 3, i,k ={1,2,3}.

В данной работе выполнен анализ трех различных подходов, которые в настоящее время применяются для решения задач трехмерной гидродинамики со свободной поверхностью. Анализ проводится как с точки зрения точности получаемых результатов, так и с точки зрения вычислительной эффективности. В связи с развитием дешевых и экономически эффективных вычислений для графических процессорных устройств (ГПУ) предполагается выполнить сравнение кодов для реализации на ГПУ

*

Работа поддержана Российским фондом фундаментальных исследований (гранты 08-07-00074а и 09-07-00078а).

фирмы NVIDIA на языке CUDA, являющемся расширением языка программирования C.

Для всех нижеперечисленных начально-краевых задач предполагается, что представлена система жидкость-газ, в которой газ не оказывает влияние на жидкость и атмосферное давление равно нулю. В связи с этим условие на свободной поверхности упрощается. Турбулентные напряжения, возникающие при соответствующих режимах течения, определяются по модели динамики больших вихрей [Evstigneev 2008], если не оговорено отдельно.

1. Общая начально-краевая задача

Определение 1. Пусть задана стационарная произвольная ограниченная область ОєШ3 с границами ГОєШ2. Найти вектор-функцию скорости V : Qx [0, T] ^ Ш3, скалярную функцию давления P : Q x [0, T] ^ Ш и положения свободной поверхности F : Qx [0,T] ^ Ш такие, чтобы для любой замкнутой области внутри Qx [0,T] они удовлетворяли бы уравнению движения вязкой несжимаемой жидкости, при условии постановки на ГО корректных краевых условий в Ш2 и постановки в Q корректных начальных условий в Ш3 для уравнений гидродинамики. Под корректными начальными условиями будем понимать такие, что V ■ V(xt ,0) = 0, P(xt ,0) є C2, F(xt ,0) є C2. Под корректными граничными условиями будем понимать такие, что V(хГ1, t) = 0,

V- V(хГ2, t) = 0, ---( Г3Г4, ) = 0, —( Г1,Г2, ) = 0, P(хГ3, t) = Const, P(xГ4, t) = 0,

^Г3, Г4 ^ Г1, Г 2

дF ( x Г1Г 3, t)

-------:---= 0, F(хГ3,t) = F(xГ3,0); здесь Гі - набор всех возможных граничных

^ Г1,Г 3

условий i={1 - твердая стенка, 2 - входное условие, 3 - выходное условие, 4 -свободная поверхность}. Для корректной начально-краевой задачи необходимо наличие хотя бы двух граничных условий.

2. Метод на основе решения уравнений Навье-Стокса

Будем рассматривать систему уравнений Навье-Стокса, записанную в интегральном виде, и начально-краевую задачу (Опр.1) для них при условии наличия фазы раздела - свободной поверхности:

— | q • ёЖ + | /. • п • ёэ - | • п - Рх - Рг - ^ = 0, (1)

д о до до

где: q = [°; м; У; - вектор переменных количества движения, ^ = {и''у, м/} в

1і • П = [(©) (м0 + п • Р); (у© + п • Р) (^0 + п • Р)]г декартовой системе координат; і ' 4 х л\ у р\ г л -

^ п = {пх, пу, п,}

невязкие потоки для произвольной подобласти внутри О; ^ х у г, - вектор

нормали к границам замкнутой подобласти;

А • п = [°; (пх^хх + пуТху + пТхг + тхм + ™хг ) (пхТух + п уТ уу + пТу, + тух + т уг )

(пхТгх + пуТгу + пТгг + ™ гх гу ^ - вектор

вязких потоков для произвольной подобласти; © = пхм + пуу + п^ ;{і, у} = {1,2,3} -

1 Ж дГ

декартовы направления; тт = — (------------1-----) - тензор вязких напряжений; g -

дх ■ дХі

] і

единичный вектор направления действия гравитации; Я = иЬ / V, Рг = и / - числа

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

осреднения для турбулентного течения модели динамики больших вихрей; р -

функция, отвечающая за поверхностное натяжение.

Для моделирования свободной поверхности для уравнений Навье-Стокса существует несколько подходов [Белоцерковский 1994]. Здесь будем рассматривать подход, который стал стандартом для задач со свободной поверхностью, а именно задание цветовой функции [Евстигнеев 2010а]. Сила поверхностного натяжения определяется как

р = акУР, (2)

где а - коэффициент смачивания; к - кривизна свободной поверхности; F -

цветовая функция, определяемая уравнением адвекционного переноса в О:

Л 7—'

— = -У(Т • Р) + РУ-V, (3)

07

УР

при этом кривизна свободной поверхности к = ^ . (4)

Поставив начально-краевую задачу для (1)-(4), необходимо разбить О на дискретные симплексы с целью записи (1)-(4) на этих полученных симплексах. Здесь будем рассматривать наиболее сложную для моделирования свободной поверхности сетку - неструктурированную сетку тетраэдров и призм. Метод решения изложен в [Евстигнеев 2010а].

3. Метод на основе гидродинамического приближения уравнений Больцмана

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

Пусть задана вероятностная функция распределения одной частицы (PPDF) f (г, е, 1), такая, что выражение [ /(г, е, 7) • ё3г • ёЗе ] представляет собой количество частиц в

момент времени t внутри фазового объема [ ё3 г • ё3 е ] в окрестности г и е, где г -координаты в физическом пространстве, е - скорость частицы. Тогда эволюция PPDF будет описываться уравнением вида:

I+e-V,+a-v e if (r,e,,)=(f

где a - внешняя сила, действующая на частицу,

coll , (5)

- член уравнения,

coll

отвечающий за соударения частиц.

Уравнение (5) переводится к форме, введенной Больцманом [Huang 1988]. Если

система находится в околоравновесном положении, то интеграл

I дt

может быть

coll

значительно упрощен. В соответствии с теорией интервальных соударений [Там же] можно сказать, что за временной интервал St часть частиц, определяемая как St/т = 1/г* (Т - время релаксации), в заданном малом объеме подвергаются соударениям и при этом их PPDF f переходят в равновесное состояние f(0), определяемое распределением Максвелла. Тогда интеграл соударения можно переписать в приближении Бхатнагара-Гросса-Крука (BGK) и для выражения силового

взаимодействия вводится предположение, что ^е1 ~ ^е1 (что в большинстве случаев верно) и тогда (1) перепишем так:

1 + е .V Г = - 1 Ї(0) + Д • (е - и)

дг г т ЯТ . (6)

Связь между микроскопической характеристикой частицы (PPDF) и макроскопическими характеристиками среды в случае гидродинамического приближения выполняется путем интегрирования f по фазовому объему, например:

Р = {[/¥е ; ри = |[/ • е]ёе ; рє = 0.5|[(е - и)2/. (7)

Рассматривая процесс эволюции (5) при малом числе Кнудсена

є = ХІЬ < 0.1 (8)

и малых значениях скорости при постоянной температуре (малое число Маха), можно применить разложение Чапмена-Энскога [Евстигнеев 2010], получив

(е - и)2

p^ exp

f(0) =

2RT

' e ■ U (e ■ U)2 U2 ^

1 +--------+

П /2 - ■ 2 + °(^3) (9)

(2лЯТ) П°/2 ^ ЯТ 2(ЯТ)2 2ЯТ

и применив (7), можно получить следующее гидродинамическое приближение для нулевых и первых моментов разложения для PPDF по числу Кнудсена:

(pU) = 0, дt V '

д(pU) , . , . _ д

(10)

+ (U ■ V)pU + V(pRT) = — ^-rV■ U

дt дх.

дх 3

i V i J

где 77 = spRT - аналог динамической вязкости. Очевидно, что при определенном выборе RT уравнения (10) переходят в уравнения Навье-Стокса для слабосжимаемой жидкости с уравнением состояния вида P = pRT, что при названных предположениях

о постоянной температуре приведет к модельным уравнениям искусственной сжимаемости, аналогичным модели Черина [Белоцерковский 1994]. При этом дополнительный член в вязкостном тензоре исчезает при V • (U) = 0, то есть при полной несжимаемости. Условие на свободной поверхности для макроскопических переменных (7) задается аналогично (2)-(4).

3.2. Метод решения

Одним из эффективных методов решения (4) является применение дискретных клеточных автоматов. Такой метод в литературе получил название ячеечного метода Больцмана [Huang 1988]. Не будем детально останавливаться на данном методе. Вся область покрывается системой симплексов, состоящих из кубов. На ней решается (4) методом клеточных автоматов, при этом для аппроксимации сложной границы вводится метод погружной границы [Евстигнеев 2010]. Рассмотрим особенности решения уравнения свободной поверхности.

В отличие от методов решения (4) в разделе 2.2, здесь автором предлагается применение самого факта переноса массы клеточным автоматом по заданным направлениям. Пусть задана начально-краевая задача (Опр. 1). Движение свободной поверхности определяется перемещением массы жидкости. Разобьем Q на область, где полностью задана жидкостью, область, где полностью задан газ и фазу раздела. Тогда будем учитывать перенос массы только для области фазы раздела, поскольку в области, занятой жидкостью, перенос массы симметричный, что следует из (4). Следовательно, расчет переноса массы и вычисления (2)-(3) необходимо проводить только в области фазы раздела.

Поскольку РРББ фактически выражают перенос массы, то тогда ее перенос может быть описан в рамках только РРББ для области фазы раздела как

Ат,(х; г + М) = (х + е,; г) - (х; г), (11)

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

Ат, (х; г + Аг) = (/. (х + е{; г) - fг (х; г)) • 0,5(е( х + е{; г) + е( х; г)), (12)

где е(х; г) = т (х; г)/ р (х; г) - доля заполнения жидкостью ячейки.

Подсчитав всю массу

т(х; г + Аг) = т(х; г) + £ Ат. (х; г + Аг), (13)

2

определяем тип ячейки после переноса. Если в данной ячейке е(х; г + Аг) > 1, то ячейка фазы раздела становится полной, а дополнительная масса Ат( х; г + Аг) = т( х; г + Аг) - р( х; г + Аг) перераспределяется по соседним ячейкам, имеющим свойства фазы раздела или пустой ячейки. Если же в данной ячейке £■(х; г + Аг) < 0, то тогда ячейка помечается как заполненная газом, а дефицит массы Ат(х; г + Аг) = т(х; г + Аг) < 0 перераспределяется по соседним ячейкам фазы раздела или полным ячейкам.

Учет граничных условий на свободной поверхности (2)-(3) переписывается в формате РРББ, таким образом, на ячейке, помеченной как фаза раздела, функция распределения со стороны газа определяется так:

(х; I + А!) = !, т(рл ,и) + /™(рА; и) - /, (х; I). (14).

Тут предполагается, что плотность газа рА = 1 и кинематическая вязкость газа много больше кинематической вязкости жидкости. А скорость в газе переносится со

ди

стороны жидкости, то есть —

дп

= 0 . Перенос со стороны газа будет выполняться,

п • е > 0; п = — 2

если

е([ х +1; у; г]; г) - е([ х -1; у; г ]; г)

а([х;у +1; г];г)-е([х;у -1; г];г) . (15)

_е([х; у; г +1]; г) - е([х; у; г -1]; г)_

Таким образом, выполняется условие поддержания тонкого слоя фазы раздела и консервативность уравнения переноса (4).

4. Метод на основе гидродинамики сглаженных частиц

Будем решать начально-краевую задачу (Опр. 1) следующим образом.

Пусть в О есть частица ‘1’, с координатами х;, постоянной массой ш; и некоторым свойством А;. Тогда, если известны значения этого свойства для частиц ‘^ в некоторой окрестности И, то А; определяется так [Ке1а§ег 2006]:

А,

А( х), = £

т— • Ж (х. - х.; И) Р ]

(16)

Здесь р - плотность частицы; W - некоторая сглаживающая функция. Сумма ведется по всем частицам ‘]’, для которых ||х. - х}. || < 2И .

С учетом симметрии взаимодействия, сглаживающая функция W выбирается так, чтобы

4

|Ж(х, - х,;к)йх = 1; Ж(х, - х,;Л) = 0, при ||х, - х,|| > 2Л. (17)

х - х,

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

ШУ УР У2У Р • § ,1С.

— +----= <и----+-----^. 08)

ш р р р

При этом закон сохранения массы выполняется строго и автоматически (простое

суммирование по количеству частиц в О), что является значительным плюсом для

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

представленных выше.

Очевидно, что в связи с конечностью окрестности И получаемое уравнение

сохранения массы имеет сжимаемость. При этом плотность жидкости

Р(х), = £т, •Ж(Х - х,;Л)’ (19)

1

а давление Р = С 2(р -р0), где С - скорость распространения возмущений в

жидкости, р0 - плотность жидкости в равновесном состоянии. Тогда сжимаемость можно определить так:

1 = (А -Ро) = (А -Ро) 1 = 1 (20)

Р к Р - Р0 Р0 Р С 2р’ К }

и для выбранной жидкости, например для воды, приняв модуль объемной

упругости к, можно найти требуемое значение С = — для удовлетворения

\Р0

физического критерия сжимаемости, в отличие от метода 3, в котором сжимаемость является исключительно численной особенностью алгоритма разложения.

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

Р + Р

/.г"" = -£т,-^-•УЖ(х, -х/;Л), (21) ,

, 2а

и для силы трения

У + У/Г = м£т,----------1 • У2Ж(х, -х,;Л). (22)

, Р,

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

После вычисления цветовой функции

^ т

р, = £^-• ж(х, -х,;Л) (22)

1 А

все остальные величины, отвечающие за взаимодействия на свободной поверхности, находятся аналогично (2)-(3).

Очевидно, что производительность метода пропорциональна 0(К2), где N -количество частиц, в связи с тем что необходимо находить окрестность сглаживания. Но, поскольку И конечна, можно применять алгоритм секционного поиска, основанный на разбиении О на вспомогательные ячейки со стороной порядка И, несущие информационную нагрузку. Такие ячейки быстро определяют местонахождение

частицы в Q и значительно снижают вычислительные затраты, до O(M- N), где M -количество частиц, для которых ||хг. - Xj || < 2h.

5. Особенности реализации алгоритмов на графическом процессоре

Все вышеперечисленные алгоритмы реализованы на платформе графических процессоров NVIDIA, и именно на ней будет выполняться сопоставление результатов. Методы решения задач гидродинамики со свободной поверхностью, описанные в разделах 2 и 3, реализованы автором ранее [Евстигнеев 2009, 2010а]. Метод решения задач на основе гидродинамики сглаженных частиц разрабатывается автором в настоящее время. Для всех задач использовался графический процессор NVIDIA TESLA C1060, установленный в лаборатории 11-3 ИСА РАН.

Для всех задач исполняемый код переписывается с учетом архитектуры графического процессора. Выделяются области, для которых используется быстрая разделяемая память устройства. Так, с точки зрения алгоритмов наиболее приспособленным к распараллеливанию под графический процессор является алгоритм сеточного метода Больцмана, поскольку он обладает локальностью и не требует значительного обмена данных в рамках различных потоков. Самым сложным является алгоритм решения уравнений Навье-Стокса, поскольку он включает в себя отдельно решение уравнения Пуассона [Евстигнеев 2009], решение нелинейного уравнения конвективного переноса, а также нелокальный обмен данных.

6. Сопоставление решения некоторых начально-краевых задач

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

6.1. Обрушение столба жидкости. Классическая задача, для которой имеется произвольный разрыв (одна часть области расчета задается как занятая жидкостью, другая - газом) и рассматривается его эволюция во времени. Сравнение проводилось по трем кодам и с литературными данными [Evstigneev 2008; Zhou, Kat, Butchner 1999]. Расчетная область представляла из себя прямоугольный параллелепипед, со сторонами x, y, z = 800; 800; 600. Вектор гравитации задан по направлению оси Z. Для метода Навье-Стокса и для сеточного метода Больцмана количество ячеек составило 2,5 млн., для метода частиц - 700 тыс. частиц. Сравнение показало, что наибольшее совпадение с экспериментом имеет метод гидродинамики сглаженных частиц (отличие профиля свободной поверхности составляет 2,5-3%). Метод, основанный на решении уравнения Навье-Стокса, показал отличие от профиля в 5%, а сеточный метод Больцмана имеет локально отличие от экспериментальных данных [Zhou, Kat, Butchner 1999] до 14% в связи с эффектами сжимаемости. Так, при ударе обрушенной массы жидкости о противоположную стенку в методе для уравнений Больцмана наблюдается сжатие жидкости и, таким образом, потеря кинетической энергии.

Необходимо отметить, что время интегрирования одного шага по времени для каждого метода значительно отличается. Самый бастрый - сеточный метод Больцмана

- 0,094 сек/шаг. Далее идет метод сглаженных частиц - 0,367 сек/шаг, и далее метод на основе уравнений Навье-Стокса - 1,2 сек/шаг. Таким образом, сеточный метод Больцмана быстрее метода сглаженных частиц в 3,9 раза и быстрее метода решения полных уравнений Навье-Стокса в 12,7 раза.

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

совпало. Сравнение амплитуды затухающих колебаний свободной поверхности для различных начальных условий показало, что наибольшее отличие имеет сеточный метод Больцмана. Отличие от метода решения полных уравнений Навье-Стокса составило 1,2% в начальный момент и 9% к моменту затухания волны. Метод сглаженных частиц имеет отличие от метода Навье-Стокса до 5% к моменту затухания волны.

6.3. Образование опрокидывающейся волны. Задача аналогична 6.1, только область расчета выбирается так, что часть ее заполнена жидкостью во много раз меньше характеристической глубины течения, с тем чтобы на ней образовывались волны на мелкой воде. Результат моделирования показал, что наиболее близкими значениями к физически наблюдаемой картине обладает метод сглаженных частиц (см. рис.). Время распространения и время глубины отличаются незначительно, но частицы позволяют показать полностью образование сопутствующих брызг и неустойчивость фронта волны в момент опрокидывания. Сеточный метод Больцмана показал худшие результаты, при этом отличие момента обрушения волны от других методов составило 12,5%. Картина скорости расчета аналогична задаче 6.1.

Сопоставление результатов продольного разреза уровней фазы раздела для задачи обрушения столба жидкости с образованием опрокидывающейся волны, шаг 0,25 сек: сверху - метод решения полных уравнений Навье-Стокса, снизу - метод сглаженных частиц

Приведенный анализ показывает, что для приложений, в которых важно точное отслеживание свободной поверхности, наиболее экономным с точки зрения

вычислений на графических процессорах является сеточный метод Больцмана, ввиду присущей клеточным автоматам локальности. При этом результаты, полученные этим методом, страдают от чрезмерной сжимаемости среды, что наблюдалось при решении тестовых задач. Однозначно можно сказать, что для моделирования явлений свободной поверхности на больших масштабах (в<0,0001) целесообразно применять лагранжевы методы (метод сглаженных частиц), которые значительно выигрывают в производительности по сравнению с методами на основе полных уравнений Навье-Стокса и при этом получают близкие, а иногда и более детальные результаты.

т

г

Заключение

Библиографический список

Белоцерковский О. М. Численное моделирование в механике сплошных сред. М.: Физматлит, 1994.

Евстигнеев Н. М. Интегрирование уравнения Пуассона с использованием графического процессора технологии CUDA // Вычислительные методы и программирование. 2010. Т. 10. С. 268-274.

Евстигнеев Н. М. Применение графического процессора для ускорения численного сеточного метода Больцмана с энтропийной стабилизацией // Труды ИСА РАН. 2010 (принята в печать).

Евстигнеев Н. М. Численный метод решения уравнений Навье-Стокса на неструктурированных сетках с применением Лагранжево - Эйлерового метода // Научно-технические ведомости СПбГПУ. 2010. 1 (93). С. 163-170.

Evstigneev N. Integration of 3D incompressible free surface Navier-Stokes equations on unstructured tetrahedral grid using distributed computation on TCP/IP networks // Proc. of the VII International conf. “Advances in Fluid Mechanics”. Oxford, 2008. 15-20 may. P.194-208,

Huang. K. Statistical Mechanics, 2nd ed. WELLY, ISBN: 978-0-471-81518-1, 1988

http://developer.download.nvidia.com/compute/cuda/2_3/toolkit/docs/NVIDIA_CUD

A_Programming_Guide_2.3.pdf

Kelager Micky. Lagrangian Fluid Dynamics Using Smoothed Particle Hydrodynamics. Department of Computer Science. University of Copenhagen, 2006.

Zhou Z. Q., Kat J. O. D. and Buchner B., A nonlinear 3-D approach to simulate green water dynamics on deck // Nantes, 1999, 7th Intl. Conf. Num. Ship Hydrodynamics.