CC BY
11
УДК 533.6.011 Вестник СПбГУ. Математика. Механика. Астрономия. 2020. Т. 7 (65). Вып. 3 МБС 76Р05, 82В40, 80А32, 82Б05
Оптимизация моделирования колебательной кинетики углекислого газа в полном поуровневом приближении*
В. И. Гориховский, Е. А. Нагнибеда
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Гориховский В. И., Нагнибеда Е. А. Оптимизация моделирования колебательной кинетики углекислого газа в полном поуровневом приближении // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2020. Т. 7(65). Вып. 3. С. 527-538. https://doi.org/10.21638/spbu01.2020.315
Численное моделирование неравновесной кинетики углекислого газа в полном поуровневом приближении представляет собой трудоемкую вычислительную задачу, которая включает решение огромной жесткой системы дифференциальных уравнений и требует оптимизированных методов для ее решения. В настоящем исследовании мы предлагаем оптимизацию для расширенного метода Гира (ЕВБР). Использование адаптивного шага по времени вместо фиксированного сокращает количество шагов в алгоритме во много тысяч раз, хотя и увеличивает время выполнения каждого шага. Использование параллельных вычислений для нахождения релаксационных членов в уравнениях кинетики позволяет еще больше сократить время расчетов. Численные эксперименты по моделированию колебательной релаксации пространственно-однородного углекислого газа проводились с помощью оптимизированных вычислительных схем разных порядков. Рекомендован оптимальный алгоритм расчетов: параллельная ЕВБР-схема четвертого порядка с адаптивным временным шагом. Этот метод требует меньше вычислительных затрат времени и памяти, а также обладает высокой стабильностью.
Ключевые слова: колебательная кинетика, углекислый газ, поуровневое приближение, численные методы, параллельные алгоритмы, оптимизация численных расчетов.
1. Введение. Моделирование неравновесной кинетики углекислого газа является одной из наиболее важных и вычислительно сложных задач в современной физической механике. Исследование неравновесных процессов в потоках углекислого газа и их влияния на параметры течений необходимо для решения задач космической газодинамики, включая вход космических аппаратов в атмосферы Марса и Венеры [1, 2], а также задач химии плазмы [3] и лазерных технологий [4]. Неравновесная физико-химическая кинетика в потоках смесей, содержащих углекислый газ, изучалась в течение последних десятилетий с использованием строгого поуровневого описания колебательной и химической релаксации [2, 5, 6], а также многотемпературных и однотемпературных подходов [7, 8].
Поуровневый подход основан на решении уравнений для заселенности колебательных уровней трех типов колебаний молекул углекислого газа и других молекул смеси совместно с уравнениями газовой динамики. Это наиболее строгий подход, но он чрезвычайно затратен в вычислительном отношении, поскольку требует ре-
* Работа выполнена при финансовой поддержке РФФИ (гранты №18-01-00493, 19-31-90059). (¡5 Санкт-Петербургский государственный университет, 2020
шения жесткой системы, содержащей более 6000 уравнений для заселенности колебательных уровней трех типов колебаний молекулы углекислого газа. Кроме того, необходимо вычислить сотни тысяч коэффициентов скорости переходов энергии для нахождения релаксационных членов в уравнениях кинетики на каждом шаге вычислений. Существующие решатели жестких систем обыкновенных дифференциальных уравнений либо не позволяют работать с системами такого размера, либо имеют недостаточную точность. Для решения такой системы в настоящей работе предложена оптимизированная схема расчета. Предлагаемый подход основан на методе Гира (backward differential formula — BDF) [9].
Численное решение жестких систем предполагает применение методов с высокой точностью и стабильностью. Однако неявные численные методы требуют больших затрат времени и памяти. Расширенный метод Гира использует предсказания о значениях в последующие моменты времени для перехода к явной схеме вычислений [10], которая позволяет численно решать систему уравнений кинетики без распределенного хранения данных и кластерных вычислений. Применение улучшенных предикторов для предсказания может повысить точность решения без потери стабильности [11]. Несмотря на то, что эта схема позволяет решать жесткие системы на обычных процессорах, подобные расчеты требуют огромного количества шагов, поскольку алгоритм предполагает использование фиксированного временного шага. В схемах BDF может применяться адаптивный временной шаг, основанный на экстраполяции Ричардсона локальной ошибки дискретизации [12]. Это позволяет уменьшить количество выполняемых шагов с миллионов до тысяч. В данной статье анализируется реализация адаптивных шагов для расширенного метода Гира.
Высокая сложность численного моделирования кинетики углекислого газа связана не только с жесткостью системы дифференциальных уравнений, но и с ее размером. Для оптимизации вычислений в настоящей работе используется эффективная параллельная структура данных [13] и алгоритмы параллельного решения разреженных систем линейных алгебраических уравнений [14]. В работе показана эффективность полученной вычислительной схемы и проведены численные эксперименты.
2. Основные уравнения. Сначала рассмотрим уравнения, описывающие колебательную релаксацию в пространственно-однородном углекислом газе. Учитываются следующие процессы, происходящие при столкновениях молекул: VT-обмены поступательной и колебательной энергией трех колебательных мод и VV-обмены колебательной энергией между модами молекулы углекислого газа, сталкивающейся с инертным партнером M:
VTi : CO2(ii, ¿2,is) + M ^ CO2(n ± 1, «2, is) + M,
(1)
VT2 : CO2(ii, «2, is) + м ^ CO2(ii, «2 ± 1, is) + M,
(2)
VTs : CO2(ii, «2, is) + M ^ CO2(ii, «2, is ± 1) + M,
(3)
VVi-2 : CO2(ii, «2, is) + M ^ CO2(ii ± 1, «2 T 2, is) + M, VV2-s : CO2(ii, «2, is) + M ^ CO2(ii, «2 ± 3, is T 1) + M,
(4)
(5)
VVi-2-s : CO2(ii, i2, is) + M ^ CO2(ii ± 1, «2 ± 1, is T 1) + M.
(6)
Здесь i i, ¿2, ¿3 — колебательные квантовые числа, соответствующие симметричной, деформационной и антисимметричной модам CO2. Считается, что партнер по столкновению (M) не меняет своего состояния в результате обменов энергией. Колебательные спектры мод рассчитывались на основе модели ангармонического осциллятора [15]. Следует отметить, что в настоящей работе в расчетах учитываются все возбужденные уровни с энергией, меньшей энергии диссоциации. Ранее в существующих работах в расчетах кинетики углекислого газа рассматривались «обрезанные» спектры трех мод молекулы CO2 [2, 16]. В настоящем исследовании химические реакции не учитываются.
Система уравнений, описывающая изменение заселенности ni1ji2ji3 (t) колебательных уровней (ii,¿2,¿3) и температуры газа T(t) со временем, имеет вид [5, 8]:
dn^(t) =щ;%,м im=0,...,lm, т =1,2,3, (7)
pU = const, (8)
где RV1ibi2 ¿3 — релаксационные члены, p — плотность, U — полная удельная энергия:
pU = p(Etr + Erot + Evibr). (9)
Полная удельная энергия может быть выражена через Etr, Erot, Evibr — поступательную, вращательную и колебательную энергии, определяемые следующим образом:
3
pEtr =-пкТ, pErot = пкТ, pEVibr = e»i,»2,»3n»i,»2,i3> (Ю)
¿1 ,i2,i3
к — постоянная Больцмана, п — числовая плотность, — колебательная энер-
гия молекулы СО2 на уровнях ¿1, ¿2,¿3.
Релаксационные члены характеризуют изменение заселенности колебательных уровней в результате УТ- и УУ-обменов в молекулах СО2:
д^ . = д^ . + д^ . + д^ . + д^-.2 + дУУТ3 + д^-.2-3. (11)
«1,«2 >®3 «1,«2>®3 «1,«2>г3 «1,»2,г3 «1,«2>г3 «1 ,.2 ,г3 «1 ,«2 >®3 У '
Релаксационные члены для каждого обмена энергией могут быть записаны в виде
ДУТ1 — П- т ■■ кМ + П , • • кМ -П- • • (кМ + кМ ) (12)
Д«1,«2 ,«3 — п«1-1,«2,«3 к«1 — 1—«1 + п«1 + 1,«2,®3 "^ + 1—«1 п«1,«2,«3 (к«1 —«1 -1 + к«1 —«1 +1/, (12)
Дг1,«2 ,«3 — П«1,«2 —1,«3 кМ — 1—¿2 + п«1 ,«2 + 1,«3 к«2 + 1—«2 — П«1,«2,г3 (кМ2—¿2 — 1 + кМ2—¿2 + 1), (13) ДУТ3 — п. . . 1кМ + п , , >М _п. . . (кМ + кМ ) (14)
д«1,«2 ,«3 — ' 1 ,«2 ,«3 —1кг3 — 1—«3 + п«1,«2,«3 +1 к.3 +1—п«1,«2,«3 (к«3—«3 — 1 + к«3—«3 + 1), (14)
T-»V Vl-2 __7.M I 7.M
Ri1,i2,i3 = nii-1,i2 + 2,i3ki1 -1,i2 + 2^i1,i2 + nii + 1,i2-2,i3 ki1 + 1,i2-2^i1
— nil,i2,i3 (ki1 ,i2 ^¿1 —1,i2 + 2 + ki1,i2^i1 + 1,i2-2
), (15)
rVV2-3 = kM I kM —
Ri1,i2,i3 = ni1.i2—3,i3 + 1ki2—3,i3 + 1^i2,i3 + ni1,i2+3,i3 — 1ki2+3,i3 — 1^i2,i3
— ni1,i2,i3 (ki2,i3^i2— 3,i3 + 1 + k^2^i3^i2+3,i3 —1), (16)
RVVl-2-3 _ n kM +
Ril,i2,i3 _ n»1 -1,i2-1,i3 + 1kii-1,i2-1,i3 + 1^il,i2,i3 + + nil + 1,i2 + 1,i3-1ki1 + 1,i2 + 1,i3-1^i1,-
21,22,23 ( ii ,22 ,23^21 -1,22-1,23 + 1 + 21,22,23^21 + 1,22 + 1,23-1 ), (17)
где k2122 23^2/ 2/ 2/ — зависящие от температуры коэффициенты скорости переходов энергии с уровня (i 1, , ¿з) на уровень (¿1, ¿2, ¿3)- Для расчета коэффициентов скорости мы используем формулы SSH-теории, обобщенные для молекул CO2 [17] -
3. Вычислительная схема. Численный метод решения жесткой системы уравнений (7), (8) должен быть точным в обширной области устойчивости. Запишем систему в виде задачи Коши:
^ = R(y,t), y(t0) = yo, (18)
где y — вектор неизвестных, состоящий из заселенности уровней и температуры газа, yo — начальное условие.
Чтобы решить эту систему, мы используем явные методы, основанные на методе Гира (BDF), введенном в [9]. Метод BDF очень точен, но требует решения большой системы линейных алгебраических уравнений- Это приводит к высокому потреблению памяти при расчетах. Следовательно, необходимы явные методы с сопоставимой точностью и стабильностью. Эти методы используют предсказания значений в последующие моменты времени для решения сложных задач и известны как расширенный метод Гира (extended backward differential formula — EBDF) [10]. Оптимизация EBDF с различными предикторами анализировалась в [11].
Несмотря на высокую точность, расширенные явные методы трудно применять из-за необходимости проведения большого числа шагов. Для сокращения времени вычислений могут быть использованы адаптивные шаги по времени [12]. В работе применяется метод адаптивных шагов для схемы EBDF. Для дальнейшего сокращения времени расчетов используются параллельные вычисления для нахождения релаксационных членов и решения рассматриваемых систем линейных алгебраических уравнений для заселенности колебательных уровней CO2.
3.1. Расширенным метод Гира. n-й шаг метода Гира порядка k выражается следующим образом:
к
аУп+j _ AtR(yn+k), (19)
j=0
где At — шаг по времени.
Для повышения стабильности BDF-схема была расширена до EBDF путем введения предсказания правых частей уравнений в последующих точках:
к
У^дУп+j _ AtAkд(Уп+к) + Д£вк+1Й(Уп+к+1), (20)
j=o
где коэффициенты в взяты из [11].
Предполагая известными yn, yn+1,..., yn+k-1, один шаг EBDF-схемы можно разбить на подшаги:
1) вычисление первого предиктора уп+к как решения системы линейных уравнений
к
У^д = Д^Дп+к; (21)
¿=0
2) вычисление второго предиктора уп+к+1 как решения новой системы линейных уравнений
к
У^дyn+j+1 _ AtRn+k+1; (22)
j=o
3) оценка предикторов релаксационных членов по решению уравнения (22)
Rn+k+1 _ R(tn+k+1 ,yn+k+1); (23)
4) решение системы (20) с предсказаниями значений релаксационных членов (23) для получения решения yn+k EBDF-схемы.
Такая схема имеет порядок k + 1 и локальную ошибку дискретизации (local truncation error — LTE) на (n + к)-м шаге [11]:
LTEn+k = Aifc+2 (l - ^zl) д-§У{к+1) + ^y(fc+2)) (tn) + o(Atk+% (24)
где C — локальная ошибка дискретизации (19) и C2 — локальная ошибка дискретизации (20).
Несмотря на высокую точность метода, использование постоянного шага по времени приводит к большому времени расчета. Явный алгоритм прост в реализации и распараллеливании, но требует дополнительного анализа для каждого набора входных данных при нахождении оптимального шага по времени. Для обеспечения стабильности метода временной шаг необходимо делать небольшим. Система (7), (8) содержит более шести тысяч уравнений. Из-за этого применение постоянного временного шага потребует миллионов шагов алгоритма и огромного времени вычисления. Поэтому для оптимизации расчетов необходимо использовать переменный временной шаг. Выбор адаптивного шага по времени позволяет достичь желаемой точности при относительно низких вычислительных затратах.
3.2. Адаптивный временной шаг. Адаптивная стратегия временного шага основана на оценках локальных ошибок дискретизации (24) [12]:
LTEn+k _ Aik+V(yk,tk) + o(Atk+3). (25)
На основании этой оценки можно построить экстраполяцию Ричардсона. Эта экстраполяция позволяет выбрать наибольший временной шаг с ошибкой, не превосходящей LTEn+k. Подходящая норма и допуск TOL указываются пользователем. Результат должен соответствовать критерию «error per step»:
||LTE„+k||< TOL. (26)
Исходя из этого критерия, если At уменьшается до величины 7At, то значение LTEn+k уменьшается до Yk+2LTEk. Таким образом, наибольшее значение 7, для которого локальная ошибка дискретизации удовлетворяет критерию «error per step», определяется выражением
1
Л ( TOL
1 = (рщ) ■ (27)
Если оценка нормы LTE на шаге n+k меньше допуска, временной шаг можно увеличить до Atn+k+1 _ 7Atn+k. Глобальная ошибка этих схем меньше TOL• Ns, где Ns — общее число шагов. Алгоритм выбора шага по времени может быть сформулирован следующим образом:
1) нахождение решения y1,..., ym и временного шага Atm;
2) проведение шага EBDF-схемы;
3) вычисление оценки LTEm;
4) вычисление 7 на основе экстраполяции Ричардсона;
5) вычисление «оптимального» шага по времени для следующего шага EBDF-схемы.
Схема EBDF с таким алгоритмом выбора шага по времени (adaptive timestep extended backward differential formula — AT-EBDF) позволяет решать жесткие системы дифференциальных уравнений с заданной точностью. Однако такой подход требует много времени и памяти. Поэтому для оптимизации расчетов целесообразно использовать параллельные вычисления.
3.3. Параллельная версия AT-EBDF. Для оптимизации численного моделирования кинетики углекислого газа в работе используются параллельные структуры данных и распределенные алгоритмы. Самыми дорогостоящими по времени и памяти в AT-EBDF-схеме являются этапы, связанные: 1) с вычислением релаксационных членов и 2) с вычислением предикторов и неизвестных y, получаемых из решения соответствующих линейных систем уравнений.
Нахождение релаксационных членов основано на расчете коэффициентов скорости переходов энергии. Эффективная параллельная структура данных для хранения и расчета коэффициентов скорости была предложена в нашей предыдущей работе [13]. Эта структура данных использует кусочно-линейное хранение данных и единовременное выделение потоков, что позволяет вычислять полные наборы коэффициентов скорости переходов энергии, используя формулы SSH-теории, без излишних затрат времени на обращение к данным и инициализацию параллельных вычислений.
Для ускорения алгоритма решения линейной системы можно использовать факт разреженности матриц в BDF-схеме. Метод подматриц, основанный на LU-разложении, рассматривается в [14]. Алгоритм решения линейных систем можно сформулировать следующим образом:
1) построение подматриц;
2) идентификация нулевых подматриц;
Input: initial conditions, intitial time step,
Parallel evaluate right hand parts: R(Tk_r
Parallel calculate predictors nj, tJ as a solution of the k step BDF-4
Calculate time step
1
Parallel evaluate right hand parts: R(Tj)
Parallel calculate solution n , Tk
no ^—• 6heck ~~ —^^^ ■еСТ"""" equilibrium
—.^.^condition^^—-----
I yes done
Рис. 1. Основной алгоритм параллельного расширенного метода Гира с адаптивным временным шагом.
3) выполнение операций с подматрицами;
4) сборка результирующей матрицы.
На основе этих оптимизаций алгоритм AT-EBDF был преобразован для параллельных вычислений. В общем виде основной алгоритм Parallel AT-EBDF-схемы представлен на рис. 1. Все линейные системы решаются с использованием метода параллельного обращения разреженной матрицы на основе метода подматриц. Релаксационные члены оцениваются с помощью оптимальной структуры данных, основанной на формулах SSH-теории.
4. Численный эксперимент. Для проверки эффективности вычислительной схемы были проведены численные эксперименты. Все расчеты проводились на процессоре Intel® CoreTM i7-4770k. При проведении экспериментов использовалась система с активированной технологией Hyper-Threading. Из доступных 4 ядер и 8 потоков при параллельных расчетах были задействованы 3 ядра и 6 потоков. Схемы были реализованы на языке C+—+ с использованием потоков POSIX.
В эксперименте сравнивались различные EBDF-схемы:
• EBDF-3, EBDF-4, EBDF-5 — схемы с постоянным временным шагом и порядками точности 3, 4, 5 соответственно;
• AT-EBDF-3, AT-EBDF-4, AT-EBDF-5 — схемы с адаптивным временным шагом и порядками точности 3, 4, 5;
• Parallel AT-EBDF-3, AT-EBDF-4, AT-EBDF-5 — параллельные схемы с адаптивным временным шагом и порядками точности 3, 4, 5.
1Е-9 1Е-8 1Е-7 1Е-6 1Е-5 1Е-4 1Е-3 0,01 0,1 1
t, sec
Рис. 2. Температура газа.
Рис. 3. Заселенности уровней энергии.
Во всех схемах использовалась параллельная структура данных для расчета коэффициентов скорости переходов энергии. Для анализа эффективности алгоритмов была рассмотрена задача моделирования пространственно-однородной релаксации углекислого газа. Система уравнений (7), (8) решалась совместно с уравнениями для релаксационных членов (11)-(17). Начальные данные для моделирования: температура газа To = 1000 K; давление Po = 10000 Па; начальные заселенности колебательных уровней определяются распределением Больцмана с колебательной температурой Туo = 5000 K. Начальный допуск рассчитывается на основе локальной ошибки дискретизации схем EBDF с временным шагом, равным 10-6, 10-8, 10-10с.
Результаты моделирования близки для всех рассмотренных вычислительных схем. Далее представлены результаты, полученные по схеме Parallel AT-EBDF-4. Рис. 2 показывает изменение температуры во время всего процесса релаксации для схемы с точностью, основанной на TOL = 10-6. На рис. 3 показано отношение заселенности определенных уровней к общей числовой плотности для уровней (0, 0, 0), (2, 2, 2), (5, 5, 5), (7, 7, 7), (10,10,10). Число молекул с низкоэнергетическими состояниями растет, а популяции с высокими энергиями со временем уменьшаются. Для уровня (2, 2, 2) можно видеть немонотонное изменение заселенности.
Таблица 1 показывает среднее время, затрачиваемое на вычисление одного шага алгоритма для разных схем. Время расчета одного шага близко для каждого типа временного шага. При использовании постоянного временного шага вычисление одного шага проводится быстрее, но требуется гораздо больше таких шагов, как будет показано далее. Адаптивные шаги являются в несколько раз более дорогостоящими в связи с использованием двух оценок для релаксационных членов, решения трех систем линейных алгебраических уравнений и экстраполяции временного шага. Используя распараллеливание, можно значительно сократить время вычислений для адаптивного временного шага. Также применение этой схемы для вычислений на GPU является вполне естественным для повышения скорости счета.
Таблица 1. Среднее время вычисления одного шага в рассматриваемых схемах (в мс)
EBDF-3 EBDF-4 EBDF-5
Схема с постоянным шагом 330 340 340
Схема с адаптивным шагом 2240 2250 2250
Параллельная схема с адаптивным шагом 565 570 570
Таблица 2. Количество шагов для схем с постоянным и адаптивным временным шагом
Точность EBDF-3 AT-EBDF-3 EBDF-4 AT-EBDF-4 EBDF-5 AT-EBDF-5
КГЬ ~ 1.5 • 10ь 1476 ~ 3.6 • 10ь 1311 ~ 7.1 • 104 1198
ю-у ~ 5.5 • 10й 1556 ~ 4.4 • 10' 1401 ~ 2.1 • 10ь 1362
10"^ ~ 7.2 • 10iu 1721 ~ 2.4 • 10й 1492 ~ 3.7- 10» 1502
Таблица 2 показывает количество шагов схем с постоянным и адаптивным временным шагом. При использовании постоянного временного шага необходимое количество шагов значительно растет при увеличении требуемой точности. Для схем с адаптивным временным шагом необходимо значительно меньшее количество шагов для достижения той же точности.
Основываясь на результатах численных экспериментов, можно сделать вывод, что различные схемы AT-EBDF дают аналогичные результаты с точки зрения времени. Расчет по схеме AT-EBDF-3 выполняется на 15 % дольше, чем другие. Оставшиеся две схемы дают близкую эффективность, но схема AT-EBDF-5 использует на 20 % больше оперативной памяти без дополнительного выигрыша по времени для хранения систем линейных алгебраических уравнений, коэффициентов скорости реакций и предикторов релаксационных членов. Схемы более высокого порядка могут привести к потере устойчивости методов [11].
5. Заключение. Численное моделирование неравновесной кинетики углекислого газа является трудоемкой вычислительной задачей, которая включает решение жесткой системы, состоящей из более чем 6000 дифференциальных уравнений. Для решения этой задачи требуются оптимизированные явные методы. В работе предложено несколько методов оптимизации для расширенного метода Гира. Первая оптимизация — это переход от постоянного временного шага к адаптивному временному шагу. Это уменьшает количество шагов алгоритма во много тысяч раз, при этом сложность выполнения одного шага увеличивается в 6 раз. Следующее усовершенствование расчетной схемы использует параллельные вычисления для еще большего сокращения времени вычислений. Использование параллельной структуры данных для нахождения и хранения коэффициентов скорости обмена энергии, эффективный расчет релаксационных членов и использование распределенных вычислений при решении линейных систем сокращает затрачиваемое время до четырех раз.
Численные эксперименты по моделированию колебательной релаксации в пространственно-однородном углекислом газе проводились для оптимизированных вычислительных схем разных порядков. Был выбран наиболее оптимальный по времени и памяти алгоритм вычислений: параллельная схема, основанная на расширенном методе Гира четвертого порядка с адаптивным временным шагом. Численный эксперимент показал, что этот метод обладает высокой эффективностью по времени и памяти, а также хорошей стабильностью. Более того, расширенный метод Гира может быть естественным образом адаптирован для вычислений на графических процессорах, поскольку его основная часть — решение систем линейных алгебраических уравнений.
Литература
1. Park C., Howe J. T., Howe R.L. Review of Chemical Kinetic Problems of Future
NASA Missions, II: Mars Entries // J. Thermophys. Heat Transfer. 1994. Vol.8, no. 1. P. 9-23. https://doi.Org/10.2514/3.496
2. Kustova E. V., Nagnibeda E.A., Armenise I. Vibrational-Chemical Kinetics in Mars Entry Problems // The Open Plasma Physics Journal. 2014. Vol.7, N Suppl 1: M5. P. 76-87. https://doi.org/10.2174/1876534301407010076
3. Silva T., Grofulovic M., Klarenaar B.L.M., Morillo-Candas A.S., Guaitella O., Engeln R., Pintassilgo C.D., Guerra V. Kinetic study of low-temperature CO2 plasmas under non-equilibrium conditions. I. Relaxation of vibrational energy // Plasma Sources Sci. Technol. 2018. Vol. 27, N 1. https: //doi.org/10.1088/1361-6595/aadb60
4. Невдах В. В., Ганджали М. Колебательная кинетика активных сред непрерывных CO2-лазеров // Журнал прикладной спектроскопии. 2005. Т. 72, №1. С. 72-79.
5. Kustova E. V., Nagnibeda E. A. State-to-State Theory of Vibrational Kinetics and Dissociation in Three-Atomic Gases // Rarefied Gas Dynamics. 2001. Vol. 585. P. 620-627. (AIP Conference Proceedings.) https://doi.org/10.1063/L1407618
6. Sahai A., Lopez B.E., Johnston C. O., Panesi M. Novel approach for CO2 state-to-state modeling and application to multi-dimensional entry flows // AIAA SciTech Forum — 55th AIAA Aerospace Sciences Meeting. 2017. Vol. AIAA 2017-0213. P. 20. https://doi.org/10.2514/6.2017-0213
7. Kustova E. V., Nagnibeda E. A. On a correct description of a multi-temperature dissociating CO2 flow // Chem. Phys. 2006. Vol.321. P. 293-310.
8. Kustova E. V., Nagnibeda E. A. Kinetic model for multi-temperature flows of reacting carbon dioxide mixture // Chem. Phys. 2012. Vol.398. P. 111-117. https://doi.org/10.1016 /j.chemphys.2005.08.026
9. Gear C. W. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice Hall PTR, 1971.
10. Cash J. R. On the Integration of Stiff Systems of O.D.E.S Using Extended Backward Differentiation Formulae // Numer. Math. 1980. Vol.34. P. 235-246.
11. Alberdi E., Anza J.J. A Predictor Modification to the EBDF Method for Stiff Systems // Journal of Computational Mathematics. 2011. Vol.29. P. 199-214.
12. Shampine L. F. Error Estimation and Control for ODEs //J. Sci. Comput. 2016. Vol. 25. P. 3-16. https://doi.org/10.1007/s10915-004-4629-3
13. Гориховский В. И., Нагнибеда Е. А. Коэффициенты скорости переходов колебательной энергии при столкновениях молекул углекислого газа: оптимизация вычислений // Вестник С.-Петерб. ун-та. Математика. Механика. Астрономия. 2019. Т. 6(64). Вып. 4. С. 659-670. https://doi.org/10.21638/11701/spbu01.2019.411
14. Lass M., Mohr S., Wiebeler H., Kiihne T.D., Plessl C. A Massively Parallel Algorithm for the Approximate Calculation of Inverse P-th Roots of Large Sparse Matrices // PASC'18: Proceedings of the Platform for Advanced Scientific Computing Conference. 2018. No. 7. P. 1-11. https://doi.org/10.1145/3218176.3218231
15. Herzberg G. Infrared and Raman Spectra of Polyatomic Molecules. D. Van Nostrand Company, Inc., 1951.
16. Armenise I., Kustova E. V. Mechanisms of Coupled Vibrational Relaxation and Dissociation in Carbon Dioxide // J. Phys. Chem. 2018. Vol. 122. P. 5107-5120. https://doi.org/10.1021/acs.jpca.8b03266
17. Шварц Р. Н., Славский З. И., Герцфельд К. Ф. Расчет времени колебательной релаксации в газах. В кн.: Газодинамика и теплообмен при наличии химических реакций. М.: Наука, 1962. С. 399-420.
Статья поступила в редакцию 6 февраля 2020 г.;
после доработки 15 марта 2020 г.; рекомендована в печать 19 марта 2020 г.
Контактная информация:
Гориховский Вячеслав Игоревич — аспирант; v.gorikhovskii@spbu.ru
Нагнибеда Екатерина Алексеевна — д-р физ.-мат. наук, проф.; e_nagnibeda@mail.ru
Optimization of CO2 vibrational kinetics modeling in the full state-to-state approach*
V. I. Gorikhovsky, E. A. Nagnibeda
St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation
For citation: Gorikhovsky V. I., Nagnibeda E. A. Optimization of CO2 vibrational kinetics modeling in the full state-to-state approach. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2020, vol. 7(65), issue 3, pp. 527-538. https://doi.org/10.21638/spbu01.2020.315 (In Russian)
Numerical modeling of nonequilibrium state-to-state carbon dioxide kinetics is a challenging time-consuming computational task that involves solving a huge system of stiff differential equations and requires optimized methods to solve it. In the present study, we propose and investigate optimizations for the extended backward differential formula (EBDF) scheme. Using adaptive timesteps instead of fixed ones reduces the number of steps in the algorithm many thousands of times, although with an increase in step complexity. The use of parallel computations to calculate relaxation terms allows one to further reduce the computation time. Numerical experiments on the modeling of spatially homogeneous carbon dioxide vibrational relaxation were performed for optimized computational schemes of different orders. Based on them, the most optimal algorithm of calculations was recommended: a parallel EBDF-scheme of fourth-order with an adaptive timestep. This method takes less computational time and memory costs and has the high stability.
Keywords: vibrational kinetics, carbon dioxide, parallel algorithms, state-to-state approach, optimization of numerical calculations.
References
1. Park C., Howe J.T., Howe R. L., "Review of Chemical Kinetic Problems of Future NASA Missions, II: Mars Entries", J. Thermophys. Heat Transfer 8(1), 9-23 (1994). https://doi.org/10.2514/3.496
2. Kustova E. V., Nagnibeda E. A., Armenise I., "Vibrational-Chemical Kinetics in Mars Entry Problems", The Open Plasma Physics Journal 7, N Suppl 1: M5, 76-87 (2014). https://doi.org/10.2174/1876534301407010076
3. Silva T., Grofulovic M., Klarenaar B.L. M., Morillo-Candas A. S., Guaitella O., Engeln R., Pintassilgo C.D., Guerra V., "Kinetic study of low-temperature CO2 plasmas under non-equilibrium conditions. I. Relaxation of vibrational energy", Plasma Sources Sci. Technol. 27(1) (2018). https: //doi.org/10.1088/1361-6595/aadb60
4. Nevdakh V., Ganjali M., "Vibrational kinetics of the active media of CW CO2 lasers", Journal of Applied Spectroscopy 72, 75-83 (2005). https://doi.org/10.1007/s10812-005-0034-4
5. Kustova E. V., Nagnibeda E. A., "State-to-State Theory of Vibrational Kinetics and Dissociation in Three-Atomic Gases", Rarefied Gas Dynamics 585, 620-627 (AIP Conference Proceeding, 2001). https://doi.org/10.1063/1.1407618
6. Sahai A., Lopez B.E., Johnston C.O., Panesi M., "Novel approach for CO2 state-to-state modeling and application to multi-dimensional entry flows", AIAA SciTech Forum — 55th AIAA Aerospace Sciences Meeting AIAA 201 7-0213, 20 (2017). https://doi.org/10.2514/6.2017-0213
7. Kustova E. V., Nagnibeda E. A., "On a correct description of a multi-temperature dissociating CO2 flow", Chem. Phys. 321, 293-310 (2006).
8. Kustova E. V., Nagnibeda E. A., "Kinetic model for multi-temperature flows of reacting carbon dioxide mixture", Chem. Phys. 398, 111-117 (2012). https://doi.org/10.1016/jxhemphys.2005.08.026
9. Gear C. W., Numerical Initial Value Problems in Ordinary Differential Equations (Prentice Hall PTR, 1971).
10. Cash J.R., "On the Integration of Stiff Systems of O.D.E.S Using Extended Backward Differentiation Formulae", Numer. Math. 34, 235-246 (1980).
11. Alberdi E., Anza J. J., "A Predictor Modification to the EBDF Method for Stiff Systems", Journal of Computational Mathematics 29, 199-214 (2011).
*The work is supported by Russian Foundation for Basic Research (grants no. 18-01-00493, 19-3190059).
12. Shampine L. F., "Error Estimation and Control for ODEs", J. Sci. Comput. 25, 3-16 (2016). https: //doi.org/10.1007/s10915-004-4629-3
13. Gorikhovskii V. I., Nagnibeda E. A., "Energy exchange rate coefficients in modeling carbon dioxide kinetics: calculation optimization", Vestnik St. Petersb. Univ. Math. 52, 428-435 (2019). https://doi.org/10.1134/S1063454119040046
14. Lass M., Mohr S., Wiebeler H., Kühne T.D., Plessl C., "A Massively Parallel Algorithm for the Approximate Calculation of Inverse P-th Roots of Large Sparse Matrices", PASC'18: Proceedings of the Platform for Advanced Scientific Computing Conference (7), 1-11. https://doi.org/10.1145/3218176.3218231
15. Herzberg G., Infrared and Raman Spectra of Polyatomic Molecules (D. Van Nostrand Company, Inc., 1951).
16. Armenise I., Kustova E. V., "Mechanisms of Coupled Vibrational Relaxation and Dissociation in Carbon Dioxide", J. Phys. Chem. 122, 5107-5120 (2018). https://doi.org/10.1021/acs.jpca.8b03266
17. Shwartz R. N., Slavsky Z. I., Herzfeld K.F., "Calculation of Vibrational Relaxation Times in Gases", in: Thermodynamic Properties of Individual Substances, 399-420 (Nauka Publ., Moscow, 1962). (In Russian)
Received: February 6, 2020 Revised: March 15, 2020 Accepted: March 19, 2020
Authors' information:
Viacheslav I. Gorikhovskii — v.gorikhovskii@spbu.ru Ekaterina A. Nagnibeda — e_nagnibeda@mail.ru
ХРОНИКА
12 февраля 2020 г. на заседании секции теоретической механики им. проф. Н. Н. Поляхова в Санкт-Петербургском Доме ученых РАН выступили: 1) доктор физ.-мат. наук, профессор А. А. Тихонов (Санкт-Петербургский государственный университет) с докладом на тему «Дом ученых: начало пути. Документы и факты»; 2) доктор физ.-мат. наук, профессор М. П. Юшков и кандидат физ.-мат. наук, доцент Г. В. Па-вилайнен (Санкт-Петербургский государственный университет) с докладом на тему «Очерк истории секции теоретической механики».
Краткое содержание докладов:
Заседание секции посвящено 100-летию Дома ученых им. А. М. Горького РАН. Оба доклада приурочены к знаменательному юбилею первого в России Дома ученых. Анализируются предпосылки к возникновению Дома ученых в стенах дворца Великого князя Владимира Александровича. Рассматривается личный вклад выдающихся представителей отечественной науки и культуры, принимавших участие в формировании Дома ученых — от идеи до реализации. Деятельность Дома ученых в период становления иллюстрируется с использованием многочисленных архивных документов. Цитируются малоизвестные публицистические произведения А. М. Горького, не вошедшие в полное собрание его сочинений. Анализируется роль и значение Санкт-Петербургского Дома ученых в жизни отечественной науки и культуры. Выясняются причины успехов в становлении Дома ученых в 1920-1921 гг. Дается очерк истории секции теоретической механики им. проф. Н. Н. По-ляхова.
