CC BY
12
УДК 621
О МОДЕЛИРОВАНИИ ТЕЧЕНИЯ ВОЗДУШНЫХ МАСС, ОСЛОЖНЕННЫХ ФИЗИКО-ХИМИЧЕСКИМИ ПРОЦЕССАМИ, ПРИ БОЛЬШОМ ЧИСЛЕ РЕАГИРУЮЩИХ ВЕЩЕСТВ С ИСПОЛЬЗОВАНИЕМ ГРАФИЧЕСКИХ УСКОРИТЕЛЕЙ
А.В. ЕВСЕЕВ, асп.
Предлагается постановка задачи моделирования течения воздушных масс в районе над лесным пожаром с множеством реагирующих веществ. Используются уравнения Навье-Стокса и динамики концентрации веществ, решаемые методом расщепления, а также уравнения химической кинетики, вычисляемые методом Гира. Предлагаются параллельные схемы методов в системе CUDA для ускорения вычислений.
Ключевые слова: численное моделирование, вычислительная гидродинамика, уравнение Навье-Стокса, уравнение Пуассона, метод прогонки, CUDA, метод Гира.
ON SUMILATION OF AIR MASS FLOWS DURING PHYSICAL AND CHEMICAL PROCESSES WITH A LARGE NUMBER OF REACTING SUBSTANCES USING GRAPHICS ACCELERATOR
А.У. EVSEEV, Post Graduate Student
The article describes the flow simulation of the air mass over the forest fire with the large number of reactants. The mathematical model is based on the Navier-Stokes equation, dynamic concentration of agents and chemical kinetics equations. The first and the second equations are solved using the projection method. The latter group of equations is computed by using the Gear method. CUDA Parallel system is used to speed up computations.
Key words: numerical simulation, computing hydrodynamics, Navier-Stokes equation, Poisson equation, sweep method, CUDA technology, Gear method.
Исследование физико-химических течений при большом числе реагирующих веществ представляет значительный интерес во многих приложениях, например, при сжигании топлив, моделировании процессов, происходящих при пожаре, распространении загрязнений, в газодинамических и химических лазерах и др.
Моделирование аэрогидродинамических процессов само по себе является нетривиальной задачей. Кроме того, эти процессы осложнены множественными реакциями большого количества веществ.
Таким образом, адекватное моделирование таких реагирующих сред возможно лишь при использовании многопроцессорной вычислительной техники, которая позволяет существенно сократить затраты машинного времени.
Одним из основополагающих уравнений гидродинамики является уравнение Навье-Стокса, которое традиционно записывается в системе «скорость - давление» [1, 3]:
— = -(• V)u + иАй - — VP + F,
<dt { ’ p ' (1)
divu = 0,
где й - вектор скорости; P - давление среды; F -ускорение макрочастиц среды, вызванное обменом количеством движения с соседними макрочастицами; и - кинематическая вязкость.
Для нахождения поля давления выводится уравнение Пуассона:
др = -р{д О + (0 V) О +
+2
3UX )2 + (dUy ^
дх J l дУ J
dUx dUy + 2 dU
dUz л2
dz ,
(2)
ду дх дz дх
д0у д0z -^
+2—у_о°^ -идо - дыр
дz ду
где О - дивергенция среды:
- д0 д0у д0
О = 6Ы0 = 0 + —^ +0 . (3)
дх ду дz
Таким образом, процесс решения гидродинамической задачи сводится к нахождению поля давления в будущий момент времени из известных
полей скорости Рт+1 ^ от, / = (х,у,z), а из него поля скоростей: от+1 ^ Рт+1, \ = (х,у,z).
Для моделирования поведения течения воздушных масс в районе над лесным пожаром применяют модель турбулентности, основанную на методе осреднения по Рейнольдсу (РАМЭ). Идея этого метода заключается в замене случайно изменяющихся характеристик потока (скорость, давление, плотность) суммами осредненных и пульса-
+
ционных составляющих [7]. Весьма существенным является то обстоятельство, что напряжения Рейнольдса являются случайными величинами, поэтому в расчетах используют статистические данные об их величине (модели турбулентности), которые получают путем анализа результатов эксперимента. К настоящему времени разработано значительное количество моделей турбулентности различной сложности, позволяющих оценить (смоделировать) величину турбулентных напряжений в различных условиях [4, 5].
Граничные условия для твердых и неподвижных граничных поверхностей состоят в обращении в ноль компонентов скорости и нормальной составляющей градиента давления. В том случае, если имеются входные отверстия, то для давления используются граничные условия первого рода, и оно задаётся непосредственно. Для компонентов скорости принято задавать параболический профиль для проекции скорости на ось, параллельную отверстию, и обращать в ноль проекции, перпендикулярные ему. На вытяжных отверстиях часто используются граничные условия второго рода, при котором в ноль обращается производная скорости по нормали к границе [1, 3].
Чтобы корректно описывать движение реагирующей среды в модель, добавляются уравнения диффузии реагирующих веществ и уравнения химической кинетики. Первые уравнения позволяют определить концентрации вещества в зависимости от координаты и времени, а уравнения химической кинетики позволяют описать изменение концентрации веществ в результате проходящих химических реакций (прямая кинетическая задача) [1].
Уравнения диффузии имеют следующий вид:
80, 2 80, 2 8 0, dCі
—J- + Y иі—J- = DC У —J + —J-
dt дхі Jf~1 8xf dt
J = 1,N,
(4)
ществ; член
dC,
~dT
N N ____
У LjP і - j =1-K -
(5)
і=1
і=1
где Р/ - символ /-го вещества; Цу - стехиометрический коэффициент при /-м веществе в левой части
7-й реакции; Яу - стехиометрический коэффициент при /-м веществе в правой части у-й реакции. Матрицы L = Ц] и R = [Яу] имеют размерность (МхК).
Система дифференциальных уравнений химической кинетики полностью определяется набором химических реакций. Это система обыкновенных дифференциальных уравнений первого порядка следующего вида:
K
N
K
N
-0-У RAn te-Y. LnAindm
u 1 j=1 m=1 j=1 m=1
і = Ш
■mj
m
(6)
где Ау - константа скорости у-й реакции. При рассмотрении равновесных и квазиравновесных процессов ограничимся для оценки скоростей реакций уравнением Аррениуса
Ai = QjTNj>
‘-j і RT
(7)
где Су - концентрация у-го вещества; йс. - коэффициент диффузии у-го вещества; N - число ве-
определяет изменение концен-
трации у-го вещества в результате химических реакций.
Граничные условия на твердых стенках для концентраций Су устанавливаются с учетов инертности этих поверхностей обращением в ноль производной по нормали к границе. На открытых границах используются мягкие граничные условия. Случай, когда границы не являются инертными, не рассматривался.
При моделировании химических реакций предположим, что имеется N веществ, участвующих в К реакциях, тогда система уравнений будет иметь вид [1, 6]
где Qj - предэкспонент; Nj - показатель степени в температурном множителе; T - температура проведения реакции (в градусах Кельвина); Ej - энергия активации (в ккал/моль); R = 1,987-10-3 ккал-град-1-моль-1.
Таким образом, объем вычислений сильно зависит от числа участвующих веществ и происходящих реакций. Уравнения Навье-Стокса (1) удобно решать с помощью метода расщепления по направлениям, который разбивает процесс вычисления на ряд этапов, решаемых методом одномерной прогонки. Также с помощью этого метода вычисляются уравнения диффузии веществ (4), однако присутствующий в них член изменения концентраций в результате химических реакций необходимо вычислять перед каждым шагом по времени с использованием текущего значения т. Система уравнений химической кинетики (6) достаточно часто является жёсткой, поэтому шаг интегрирования должен согласовываться с характерным временем наиболее быстрого процесса. В этом случае необходимо использовать жестко устойчивые численные методы, такие как метод Гира [6, 8].
В сфере научных вычислений все большую популярность набирают вычисления с использованием графических ускорителей. Современные видеоакселераторы представляют собой массивнопараллельные процессоры с общей памятью. В отличие от центрального процессора с несколькими ядрами, один графический процессор содержит несколько сотен ядер, которые могут проводить вычисления параллельно.
К настоящему времени существуют несколько технологий и программных библиотек для работы с GPU: NVIDIA CUDA, AMD/ATI APP (ранее ATI Stream), OpenCL, Brook++, DirectCompute, Microsoft C++ AMP. В приведенном списке первые две технологии являются архитектурозависимыми и используются только в графических ускорителях конкретного производителя.
Однако наиболее развитой на данный момент технологией является система CUDA, предложенная компанией Nvidia в 2007 году [10]. Она позволяет достаточно просто использовать графические
ускорители как массивно-параллельные процессоры (SIMD) с помощью простой и понятной модели выполнения и использования мощностей видеоакселератора. В данной технологии вычисления производятся множеством блоков, состоящих из некоторого числа обрабатывающих потоков.
Уравнение Пуассона является очень трудоемким для вычисления. Поэтому необходимо выбрать метод его интегрирования, который позволяет решать его с очень высокой точностью и эффективно распараллеливается. Одним из таких методов является метод Ричардсона. Л. Ричардсон предложил рассматривать стационарное уравнение Пуассона как эволюционное, решая его до того времени t, пока решение перестанет меняться в пределах интересующей точности (метод установления). Идея Л. Ричардсона состоит в том, чтобы на каждом шаге менять итерационный параметр т [2], основанный на полиномах Чебышева.
При построении параллельной схемы метода в системе CUDA важным является правильное и оптимальное разбиение задачи на блоки, которое естественным образом получается с использованием геометрического параллелизма. Уравнение Пуассона разбивается на ряд подзадач, для каждой из которых поставлены свои граничные условия.
Отличительной особенностью CUDA является то, что в ней нет возможности осуществить глобальную синхронизацию вычислений между блоками иным способом, кроме как закончить вычисления во всех блоках и произвести обмен данными на компьютере (хосте). В связи с этим каждый этап вычислительного алгоритма может использовать собственное, оптимальное для него разделение задачи между блоками. Однако эта особенность не позволяет полностью перенести итерационные алгоритмы на GPU, поскольку каждый блок после проведения итерации на своем участке данных и вычисления погрешности должен соотнести их с другими блоками, а это не представляется возможным. Таким образом, для итерационных алгоритмов естественным решением является вычисление одной итерации метода на устройстве, а проверка сходимости и принятие решения о продолжении расчета принимает хост.
Проверка сходимости метода заключается в выполнении редакционного суммирования по массиву сходимости. В результате получается количество элементов, не достигших сходимости [3].
Существует два основных подхода к организации вычислительного алгоритма метода прогонки:
1) применение специализированных алгоритмов параллельной прогонки, основанных на методах циклической редукции [9, 10];
2) использование обычного последовательного алгоритма прогонки каждым потоком [9].
В параллельной циклической редукции (PCR) решение находится построением линейной комбинации уравнений за один проход, при этом на каждом шаге система делится на две подсистемы, пока не останется N/2 систем с двумя неизвестными, и находится решение каждой системы (рис. 1).
4 рабочих потока
Ю 0) G) О 0) © © (Л
(О (О S3 В) (0 <й «]»—Д—
К) (£1 К) в) (»1 EHD—
GO (Р) vv G*) G*) (Р) нахождение всех
неизвестных
1одг (8) = 3 шага Рис. 1. Параллельная циклическая редукция
Алгоритм параллельной циклической редукции накладывает ограничения на максимально возможный размер одной системы, который ограничивается возможностями аппаратуры, поскольку одна система решается одним блоком, а максимально возможное число потоков в блоке, и соответственно, размер системы, равно 512 [10].
Второй подход лишен данного недостатка, но в нём появляется проблема организации массивов входных данных. Если удается использовать возможности связанной загрузки данных, то сокращается общее время ожидания загрузки и возрастает эффективность. При организации массивов данных (a, b, c, f, H) можно располагать данные последовательно, т.е. элементы массива каждой системы идут друг за другом либо перемежаясь между собой. Первый способ подходит для алгоритмов, основанных на циклической редукции, поскольку в них данные каждым блоком считываются только для одной системы. Второй способ эффективен при использовании обычного алгоритма прогонки, так как в данном случае один блок решает несколько систем одновременно, и каждый поток считывает один и тот же элемент данных только у разных систем. Таким образом, достигается связанная загрузка данных для каждого элемента (рис. 2) [9].
система 1 система 2
аО а1 а2 аЗ аО al а2 а) аО аО аО аО а1 |а11 а1 a1
Потв«1 Потса J ftorocl
а) б)
Рис. 2. Расположение данных: а - последовательное; б - перемежающееся
При вычислении метода Гира с использованием библиотеки CVODE [12] были протестированы три доступные реализации векторных операций: последовательная (serial), параллельная с использованием MPI и CUDA [3]. Результаты сравнения производительности векторных операций представлены на рис. 3. Анализ графиков показывает, что CUDA-реализация векторных операций начинает выигрывать в производительности при увеличении размерности векторов, характеризующих количество веществ и реакций между ними. Таким образом, это позволит строить химические модели многокомпонентных смесей со сложными реакциями между ними, и их число может измеряться тысячами и сотнями тысяч.
Рис. 3. Сравнение векторных операций на различных платформах
Заключение
Моделирование движения среды, осложнённой множеством химических реакций большого числа веществ, изначально предполагает наличие больших расчетных сеток для более детального изучения. Кроме того, число веществ и реакций между ними может исчисляться десятками тысяч. В этом случае использование графических ускорителей становится оправданным, поскольку позволяет их полностью загрузить вычислениями.
Предложенная методика позволяет заметно сократить накладные расходы на рассмотрение одного варианта течения и при поиске наилучшего инженерного решения, например, в области лазерных технологий, существенно увеличить число исследуемых вариантов. Также в некоторых случаях становится возможным моделирование процесса в режиме реального времени с возможностью интерактивного взаимодействия.
1. Численные методы и параллельные вычисления для задач механики жидкости, газа и плазмы / Э.Ф. Балаев, Н.В. Нуж-дин, В.В. Пекунов и др.; Иван. гос. энерг. ун-т. - Иваново, 2003.
2. Евсеев А.В. Методы решения уравнения Пуассона / ИГТА. - Иваново, 2009.
3. Евсеев А.В. О моделировании физико-химических течений при большом числе реагирующих веществ // Физикохимическая кинетика в газовой динамике. - 2011. -Т. 12 / http:// www.chemphys.edu.ru/pdf/2011-05-05-001.pdf
4. Перминов В.А. Численное решение задачи о возникновении верхового лесного пожара в трехмерной постановке // Вестник Томского государственного университета. - 2009. -Т.6. - № 1. - С. 41-48.
5. Perminov V. Numerical Solution of Reynolds equations for Forest Fire Spread // Lecture Notes in Computer Science. -2002. - V. 2329. - P. 823-832.
6. Полак Л.С., Гольденберг М.Я., Левицкий А.А. Вычислительные методы в химической кинетике. - М.: Наука, 1984.
7. Метод моделирования отсоединенных вихрей для расчета отрывных турбулентных течений: предпосылки, основная идея и примеры применения / М.Х. Стрелец, А.К. Травин, М.Л. Шур, Ф.Р. Спаларт // Научно технические ведомости. -2004. - № 2.
8. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциальноалгебраические задачи: пер. с англ. - М.: Мир, 1999.
9. Efficient Tridiagonal Solvers for ADI methods and Fluid
Simulation / Nvidia [Электронный ресурс] -
http://www.nvidia.com/content/GTC-2010/pdfs/ 2015_GTC2010.pdf
10. Zhang Y., Cohen J., Owens J.D. Fast tridiagonal solvers on the GPU // Principles and Practice of Parallel Programming. - 2010. - P. 127-136.
11. NVIDIA GPU Computing Developer / Nvidia [Электронный ресурс] - http://developer.nvidia.com/object
/gpucomputing.html
12. SUNDIALS (SUite of Nonlinear and Dlfferen-
tial/ALgebraic equation Solvers) / Lawrence Livermore National Laboratory [Электронный ресурс] -
https://computation.llnl.gov/casc/sundials/main.html
Евсеев Александр Владимирович,
ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», аспирант кафедры высокопроизводительных вычислительных систем, телефон (4932) 26-98-29.
