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

УДК 004.94

МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ СЛАБОСЖИМАЕМОЙ ЖИДКОСТИ МЕТОДОМ СГЛАЖЕННЫХ ЧАСТИЦ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ

© 2016 Е.И. Николаев

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

Статья поступила в редакцию 14.10.2016

В работе предложен алгоритм для моделирования течения идеальной жидкости при наличии в области решения твердых тел с произвольной геометрией. Алгоритм основан на применении метода гидродинамики сглаженных частиц и адаптирован для параллельных вычислительных систем на базе графических процессоров. Разработанный алгоритм реализован в виде программного комплекса на языке С++. Модификации данного алгоритма могут быть использованы для моделирования процессов газодинамики.

Ключевые слова: моделирование, численные методы, гидродинамика сглаженных частиц, графический процессор

На современном этапе развития численных методов широкое распространение получили подходы, основанные на сеточных методах решения систем дифференциальных уравнений в частных производных. Распространение и повышение доступности высокопроизводительных вычислителей на основе графических ускорителей (graphical processing unit, GPU) создает перспективы для эффективного применения методов бессеточного численного моделирования. Наиболее эффективным, с точки зрения особенностей реализации на GPU, является метод гидродинамики сглаженных частиц (smoothed particles hydrodynamics, SPH). Применение метода SPH позволяет существенно расширить круг решаемых гидродинамических задач, а также уменьшить время получения решения для проблем, решаемых сеточными методами. Метод SPH позволяет моделировать следующие явления, математическая постановка которых представлена системами дифференциальных уравнений (СДУ) в частных производных: вход в жидкость объектов сложной геометрии; обтекание потоком жидкости объектов сложной геометрии; обтекание потоком жидкости твердых или пластичных объектов, поверхностные явления в жидкостях в областях с произвольной геометрией; моделирование разрывных течений жидкостей с различными физическими характеристиками [1]. Метод позволяет комбинировать условия моделирования. Отсутствие необходимости генерации сеток, как в методе конечных разностей (МКР) или конечных элементов (МКЭ), позволяет более эффективно решать задачи с деформируемыми телами в жидкостях, подвижными границами и свободными поверхностями.

Николаев Евгений Иванович, кандидат технических наук, доцент. E-mail: [email protected]

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

Применение бессеточных методов сопряжено с высокими вычислительными затратами, но, одновременно, метод SPH характеризуется возможностью распараллеливания и отсутствием необходимости генерации сеток, что существенно снижает вычислительные затраты по сравнению с МКЭ и МКР [2]. Применение GPU для реализации основных стадий процесса моделирования позволяет сократить время получения решения, а для ограниченного круга задач - получать решение и визуализировать моделируемое течение в режиме реального времени. Реализация SPH на GPU сопряжена с необходимостью решения следующих проблем: постановка задачи моделирования, применение допущений; дискретизация исходных уравнений модели; проектирование алгоритма решения задачи с учетом высокой степени параллелизма решаемой задачи; адаптация алгоритма к возможностям вычислительной системы на базе GPU; визуализация результатов с учетом особенностей GPU-вычислителя.

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

1. Постановка задачи. Применяется модель слабосжимаемой жидкости, движение которой описывается уравнениями Навье-Стокса [3].

= -- Vp + V. (vVv) + Fx v

dt p \dt

Здесь v=(vx, Vy, Vz) - скорость жидкости в некоторой точке; р - плотность жидкости; p - давление; v - кинематическая вязкость; Fext - внешние силы, приведенные к единице массы (Fext=g, g - ускорение свободного падения). Для слабосжимае-мой жидкости давление p в точке с плотностью р получим по формуле p=kpo ((р / ро)у - 1)/у, где k -коэффициент жесткости, коэффициент у=7, ро -плотность воды.

На рис. 1 представлена схема расчетной области. Замкнутый резервуар О с границей ЭФ=Г имеет линейные размеры: высота - H, длинна - L. Препятствие fio с границей дйо=Го имеет линейные размеры: высота - h, длинна - l.

Рис. 1. Схема расчетной области

К уравнениям (1) добавляются начальные граничные условия:

^иЦ, = 0 ; РГиЦ, = 0 (2)

2. Дискретизация исходной задачи. Построение дискретного аналога уравнений (1), (2) при использовании метода БРИ основано на ла-гранжевом представлении сплошной среды множеством частиц Ту = (г(х), г(у), г(г)) (у = ТЖ, Ы-количество частиц), каждая из которых обладает своей скоростью у-„ массой ш/, плотностью р/, конечным объемом Ау. При этом дискретизация исходной функции /(г) заключается в двухэтап-ном преобразовании: аппроксимация ядром (3) и аппроксимация частицами (4).

f (r) = j f (r' )W (:

r - r'

(3)

Здесь Ж(х — х', к) - функция ядра. Интегрирование производится по области решения О. Аппроксимация ядром представляет функцию /(г) через множество значений /(т' ) , где Т, т' ей;

к - длинна сглаживания. В качестве ядра применяется функция, удовлетворяющая условиям:

(т — т', =1 Ит Ж(г — г', к) = 3(т — г')

й , к^0 ,

где 3 - дельта-функция Дирака.

(ДТ »=Е Р-Ку) ЖУ.Жу= Ж(Т — Ту, к)

] =1 Р] ; 7 7 (4)

Здесь {/(гг )} - приближенное значение функции /(г), соответствующее частице г, г = 1, N, выраженное через характеристики ш/, р/ оставшихся частиц г/. Значение производной функции аппроксимируется с использованием (5).

N

m,

(V-f (r,)) = f (r,)-VWj

j=1 pj

r — r dW-

V Ж = r- j ,j

' j

dr.

r dW-

j j

r dr■■

j j

При моделировании гидродинамических процессов бессеточными методами важное значение имеет функция ядра ж(Т — Т', к). Выбор конкретного ядра, исследование влияния Ж(т — т' , к) на точность, корректность и время получения решения, являются предметом отдельного исследования. В рамках данной работы применяется ядро (6).

W (R) = k ■

kern(R), R е [0;2], 0, R е (—да;0) U (2;+да).

(6)

Здесь параметр к определяется из соотношения Ж(т — т' , к) = 0 при |г — г'| > кк. Функция

кегп(Я) = 2 — 9Я2 +19Я3 -— Я4

3 8 24 З2 [3]. Здесь Я = г / к,

где г = |г — г'|. В двумерном случае к = 15/7як2, в

трехмерном - к = 315/208як3. При практической реализации метода БРИ аппроксимацию (4) производят не по всем N частицам г/, а используют конечное подмножество частиц

йшр = {гу | |т — г] < кк } области решения Ф (рис. 2).

Рис. 2. Область поддержки Цщр для частицы г, Для каждой частицы справедливо равенст-

во (7).

N

a=Z mW

j=1 (7)

Уравнения (1) приводятся к дискретной форме с учетом (3) - (7).

r

j

r

Q

N

avL = у dt

vi j=1

(

m

j

Pi P2

VW

J У

Vr

Vi Pj Vj Vr,

N

m,

8W„-

(8)

Здесь Vj = Vi - vj - дискретная аппроксимация относительной скорости пары частиц. Каждое уравнение (8) записано для частицы г, но вычисление осуществляется по параметрам частиц г-„ входящих в область поддержки Ц^.

3. Проектирование итерационного алгоритма решения задачи методом SPH. Обобщенный алгоритм решения дискретного аналога (7) исходной задачи представлен на рис. 3. Каждая стадия вычислительного эксперимента может быть реализована различными способами. На стадии инициализации необходимо задать геометрические размеры подобластей области моделирования; задать физические характеристики моделируемой жидкости, начальные условия вычислительного эксперимента. Для оптимизации данного этапа разработан алгоритм импорта геометрических объектов из SD-редактора Blender - это позволило сократить время подготовки геометрических объектов, присутствующих в области решения. Геометрия тел задается множеством точек, соединенных ребрами; поверхности представлены множеством треугольных плоских элементов. Через графические примитивы передаются и начальные условия вычислительного эксперимента. Например, цвет точки, заданный в редакторе Blender, интерпретируется программой моделирования как начальная скорость частицы в данной точке.

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

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

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

Инициализация

Генерация частиц

И

Поиск соседей

Вычисление ядра

Препроцессинг

Вычисление сил

Расчет вязкости

Т

Обновление параметров частиц i-

Визуализация

Рис. 3. Схема алгоритма решения задачи методом SPH

Для каждой частицы г, вычисляется функция хэширования (9).

/ьы, (г) =г(^(х)а(у) + г(^(х) + г(х)

В

(7)

r 1» _ r 1» f(.a) _ 'i 'min ,

Ja)

od dX),

d = ^(х, d(у), d(2)) - вектор, содержащий размерность сетки по каждому измерению; а - перечисление координат x, y, г; гтщ = О^П,-минимально возможные координаты по каждой оси.

После вычисления функции хэширования необходимо осуществить сортировку массива частиц по значению /м, при этом частицы, имеющие близкие хэш-функции будут располагаться в массиве рядом друг с другом. Сортировка массива частиц на каждой итерации алгоритма приводит к затратам времени, но, учитывая необходимость поиска соседних частиц для каждой частицы, механизм хэширования снижает трудоемкость поиска соседей с 0(№) до 0(с„М), с„, - среднее количество соседних частиц. На рис. 4 представлена схема однородной сетки хэш-функции для задачи, представленной на рис. 1. Показаны схематично фиктивные частицы, покрывающие поверхность дамбы в два слоя. Хэш функция не является уникальной для каждой частицы, поэтому несколько частиц могут попасть в одну ячейку.

Рис. 4. Схема однородной сетки для вычисления хэш-функций частиц

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

> » 1

К

! { и

А

Рис. 5. Применение квадро-дерева для поиска соседних частиц

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

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

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

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

- силы, действующие на каждую частицу;

- значения вязкости частиц;

- скорости частиц;

- обновленные координаты частиц;

- интегральные характеристики (момент, энергия).

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

При запуске программы в качестве параметра командной строки передается количество итераций или длительность вычислительного эксперимента. Эти данные будут переданы на GPU и передача между GPU и CPU будет прервана до момента завершения вычислительного процесса. После копирования результатов моделирования на CPU производится визуализация данных с использованием пакетов языка Python. На рис. 6 показана стадия вычислительного эксперимента, схема которого приведена на рис. 1, в момент времени i=0,01 с. На рис. 7 приводится развитие процесса обрушения столба жидкости (i=0,05 с). Развитие обрушения в момент времени i=0,09 с показано на рис. 8. На рис. 9 показана

заключительная стадия вычислительного эксперимента (в момент времени ¿=0,13 с), после которой жидкость переходит в состояние покоя.

Рис. 6. Визуализация столба жидкости при t=0,01 с

Рис. 7. Визуализация столба жидкости при t=0,05 с

Рис. 8. Визуализация столба жидкости при t=0,09 с

Рис. 9. Визуализация столба жидкости при t=0,13 с

В представленном вычислительном эксперименте использовалось 105 частиц для представления жидкости. Из-за простой геометрической формы дамбы и ее неподвижности виртуальные частицы не применялись. Полное время, затраченное на вычислительный эксперимент, составило 10 с при использовании видеокарты GeForce 560 Ti и 5,7 с при проведении моделирования на видеокарте GeForce GTX Titan.

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

СПИСОК ЛИТЕРАТУРЫ:

1. Liu, G. Spatial Smoothed Particle Hydrodynamics: A Meshfree Particle Method / G.R. Liu, M.B. Liu. - World Scientific Publishing Co. 2003. 473 p.

2. Дроздова, В.И. Адаптация конечно-разностных и конечно-элементных методов решения задачи левитации к реляционной архитектуре программных комплексов / В.И. Дроздова, Г.В. Шагро-ва, Е.И. Николаев // Наукоемкие технологии. 2012. №7. С. 59-65.

3. Суравкин, А.Ю. Реализация метода SPH на CUDA для моделирования несжимаемых жидкостей // Наука и образование: электронное научно-техническое издание 2012. №7 [Электронный ресурс]. URL:http://technomag.edu.rU/doc/423582.html.

SIMULATION THE FLOW OF WEAKLY COMPRESSIBLE LIQUID BY MEANS OF SMOOTHED PARTICLE METHOD HYDRODYNAMICS ON GRAPHICAL

PROCESSING UNITS

© 2016 E.I. Nikolaev

North-Caucasus Federal University, Stavropol

An algorithm for the simulation of ideal fluid in the presence of solutions in the field of solids with arbitrary geometry is offered. The algorithm is based on the method of smoothed particle hydrodynamics and adapted for parallel computing systems based on graphics processors. The algorithm is implemented as a software package in the language C ++. Modification of this algorithm can be used to simulate gas dynamics processes.

Key words: modeling, numerical methods, smoothed particle hydrodynamics, graphic processor

Evgeniy Nikolaev, Candidate of Technical Sciences, Associate Professor. E-mail: [email protected]