Математическое моделирование инициирования детонации в канале с профилированным торцом с использованием параллельных вычислений

Лопато Александр Игоревич

УДК 519.63

DOI: 10.18698/2309-3684-2023-4-1526

Математическое моделирование инициирования детонации в канале с профилированным торцом с использованием параллельных вычислений

ИАП РАН, Москва, 123056, Россия

Работа посвящена численным исследованиям инициирования детонации в газовой смеси в прямоугольном канале с профилированным торцом. Инициирование детонации происходит в результате взаимодействия ударных волн, которые образуются при отражении падающей ударной волны относительно малой интенсивности от торца канала. Математическая модель представляет собой систему уравнений газовой динамики, дополненную кинетикой Аррениуса для модельной водородно-кислородной смеси с табличными параметрами кинетики, отвечающими рабочему диапазону давлений и температур смеси. Численные расчеты проводятся с использованием метода конечных объемов. Построение расчетных сеток, состоящих из треугольных ячеек, осуществляется с использованием свободно распространяемого пакета SALOME. Численный алгоритм распараллелен методом декомпозиции расчетной области с использованием библиотеки METIS. Обмен сеточными функциями между вычислительными ядрами осуществляется с использованием функций библиотеки MPI. Рассмотрен вопрос ускорения реализованного в коде параллельного алгоритма по сравнению со случаем линейной зависимости ускорения от числа вычислительных ядер. Проведен ряд расчетов с использованием различного числа треугольных ячеек и сравнение картин инициирования детонации. Для всех проведенных расчетов время инициирования детонации примерно одинаковое. Основная разница в картинах детонации связана с разрешением структур, связанных с течением газа и физико-химическими реакциями в смеси.

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

Введение. Математическое моделирование газодинамических процессов с химическими реакциями, ударными волнами (УВ), детонационными волнами (ДВ) представляет интерес для фундаментальных и прикладных исследований. Физико-химические параметры течения с такими процессами могут меняться в широких диапазонах. Например, в [1] наблюдаются изменения давлений и температур на несколько порядков. Течения с волнами детонации также характеризуются сложной нелинейной природой, исследование которой сопряжено с определенными сложностями при проведении экспериментов, включая измерение характеристик течения, величины которых могут варьироваться в широких пределах.

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

включая геометрию областей, химические свойства смеси, масштабы изменения физических параметров газа. Другой проблемой, где процесс распространения ДВ играет важную роль, является разработка импульсного детонационного двигателя, в котором распространение ДВ происходит в кольцевом канале [2]. Распространение ДВ в таком двигателе сопряжено с выделением энергии, используемой в качестве полезной работы. Известно, что КПД двигателя, основанного на детонационном горении может быть на несколько десятков процентов выше, чем в случае дефлаграционного горения.

В работе [1] рассматривается процесс инициирования ДВ в результате отражения падающей УВ от профилированного торца канала. Процесс инициирования является достаточно сложным из-за ряда факторов, включая формирование новых УВ, контактных поверхностей (КП), волн разрежения (ВР) и других структур течения. Использование профилированных торцов в каналах является одним из возможных решений задачи об оптимизации процесса инициирования ДВ. Так, использование профилированных торцов позволяет уменьшить время инициирования ДВ по сравнению с использованием каналов с плоским торцом. Дело в том, что выбор геометрии торца влияет на число волн, образующихся при отражении падающей УВ, а также на интенсивность взаимодействия волн. Инициатором ДВ в таких случаях могут стать так называемые "горячие точки" - области повышенных значений давления и температуры смеси (см, например, [3]). В работе [1] рассматривается использование профилированных торцов в форме участков эллиптических кривых. Получен ряд характеристик процесса инициирования детонации, описывающих процесс в целом. Так, показано, что использование геометрии с двумя разнесенными эллиптическими участками обеспечивает меньшее время инициирования смеси по сравнению с пятью соприкасающимися участками, а также плоским торцом. В данной работе проводится численное исследование инициирования ДВ в канале с геометрией торца из [1].

Численные исследования проблем инициирования и распространения детонации в сложных областях с криволинейными границами, как правило, проводятся с использованием неструктурированных расчетных сеток. Как отмечено в [4], автоматическое построение блочно-структурированных сеток в произвольных областях с достаточной аппроксимацией границ является достаточно сложной задачей. С другой стороны, за последние десятилетия разработаны и реализованы алгоритмы построения неструктурированных сеток для областей произвольной

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

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

Постановка задачи. Рассматривается плоский канал с торцом, состоящим из трех вертикальных участков и двух эллиптических. Канал заполнен модельной водородно-кислородной смесью. Давление в канале 0.04 атм, температура — 298 К. В канал подается плоская ударная волна с числом Маха 2.5, которая распространяется в сторону торца. Учитывая, что в данной задаче имеется плоскость симметрии, будем рассматривать половину геометрии канала, содержащую один эллиптический участок, и считать ее расчетной областью задачи. Уменьшение области задачи в два раза дает возможность сократить время расчетов.

У [мм] 30

28 30

[[мм]

0

Рис. 1. Схема расчетной области. Схема содержит размеры расчетной области, постановку граничных условий, положение падающей ударной волны с числом Маха 2.5 в начальный момент времени

Геометрия расчетной области приведена на рис.1. Граничные условия устанавливаются следующим образом. Через левую стенку подается ударная волна с числом Маха 2.5. На верхней границе — условие симметрии. На нижней и правой стенках — условие непротекания. Геометрия задачи отвечает геометрии экспериментальной установки из [1].

Рассматривается процесс инициирования газовой смеси, реализующийся в результате взаимодействия падающей ударной волны с профилированным торцом канала. При построении картин течений в численных экспериментах внимание уделяется распределениям температур и градиентов плотности. Распределения градиентов плотности отвечают визуализации Шлирен-методом, которые применяется в натурных экспериментах.

Математическая модель. Математическая модель основана на двумерной системе газовой динамики, записанной в декартовой системе координат и дополненной одностадийной моделью кинетики горения Аррениуса. Система определяющих уравнений приведена ниже:

Ш д¥ дС 0 — + — + — = 8, дХ дх ду

и =

р ри рл " 0 "

ри р + ри2 рuv 0

рл , Р = рuv , С = р + рл2 , 8 = 0

е ( р + е ) и (р + е ) V 0

Р2 _ р2и р2л р(

£ =

/ . + ZQ, р = -ЯТ, ( = -Лр2 ехр (-— |. р(у- 1) ^ V ^ I ЯТ)

Здесь, Х — время, р — плотность газовой смеси, и и V — компоненты скорости частиц газа, р — давление, Т — температура, е — полная энергия на единицу объема, 2 — массовая доля реагента, ( — скорость химической реакции, £ — удельная внутренняя энергия газа, Q — тепловой эффект химической реакции, и — молярная масса смеси, Я — универсальная газовая постоянная, Л — предэкспоненциальный множитель. Рассматривается модель совершенного газа с постоянным значением показателя адиабаты у .

В работе используются следующие значения физико-химических параметров смеси:

МДж

у = 1.23, и = 12-, Q = 7.37-

моль

кг

E = 76.2 ^, A = 9.16.10"-^

моль кг•с

Значения параметров взяты из базы данных [5]. Значения, взятые из [5] отвечают давлению смеси около 0.2 атм, что соответствует давлению за ударной волной, отраженной от стенки, в постановке данной задачи с числом Маха 2.5 и давлением 0.04 атм в покоящейся смеси.

Вычислительный алгоритм. Основной элемент вычислительной технологии в решении данной задачи заключается в использовании полностью неструктурированных сеток с треугольными ячейками. При построении таких сеток используется триангуляция Делоне. Процесс построения треугольных расчетных сеток осуществляется в свободно распространяемом пакете SALOME [6].

Рис. 2. Пример декомпозиции расчетной области на 4 подобласти. Разными цветами обозначены разные подобласти. Серым цветом обозначена подобласть с номером т. Жирной линией ограничены гало-ячейки подобласти т

Вычислительный алгоритм основан на методе расщепления физических процессов Стренга. Для расчета течения за один шаг по времени сначала уравнения газовой динамики интегрируются по времени без учета химической кинетики, а затем учитывается вклад химической кинетики без учета конвекции. На этапе учета конвекции дискретизация уравнений осуществляется с использованием метода конечных объемов. Для учета потоков через грани ячеек используется поток ЛИБМ [7]. Для повышения порядка аппроксимации по координате применяется реконструкция сеточных функций с ограничителем штшоё.

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

Вычислительный алгоритм распараллелен с использованием библиотек METIS, MPI. Расчетная область разбивается на подобласти с использованием свободно распространяемой библиотеки METIS. Число подобластей выбирается равным числу вычислительных ядер, используемых в расчетах.

Будем, для определенности, считать, что подобласть с номером m, отвечает вычислительному ядру с номером id, равным m. Для определения газодинамических параметров течения в ячейках подобласти m в новый момент времени используются, как массивы значений в этих ячейках в текущий момент времени, так и массивы значений в гало-ячейках [8] подобласти m. Под гало-ячейками подобласти m понимаются ячейки из соседних подобластей, актуальные значения газодинамических параметров в которых передаются вычислительному ядру m для расчета потоков через грани ячеек соответствующей подобласти.

Пример разбиения области на 4 подобласти приведен на рис.2. Серым цветом отмечены внутренние ячейки подобласти m. Выделенной линией ограничены гало-ячейки подобласти m. Набор гало-ячеек определяется расположением границ между подобластями, а также вычислительным алгоритмом расчета реконструкций сеточных функций на этапе учета конвекции.

Учитывая, что интегрирование по времени осуществляется с использованием двухстадийного метода Рунге-Кутты, обмен значениями в гало-ячейках между вычислительными ядрами осуществляется два раза за один шаг по времени. Обмен значениями в гало-ячейках осуществляется с использованием функций библиотеки MPI. Блокирующие функции MPI_Send и MPI_Recv являются эффективными в случае, когда число соседних подобластей для любой области не выше двух.

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

Параллельные расчеты выполнены с использованием суперкомпьютеров МВС-10П и МВС-10П ОП (межведомственный суперкомпьютерный центр РАН) [9].

При выполнении расчетов распределения параметров течения в определенные моменты записываются в файлы CGNS формата. Название файла отвечает моменту времени и номеру вычислительного ядра, который записывает распределения параметров в подобласти.

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

число вычислительных ядер

Рис. 3. График зависимости ускорения реализованного в коде параллельного вычислительного алгоритма от числа вычислительных ядер. Пунктирная линия отвечает теоретическому ускорению, которое на практике не достигается из-за ряда факторов. Точками отмечены измеренные в ходе расчетов значения ускорений.

Сплошная линия состоит из отрезков прямых, соединяющих точки

Результаты моделирования. Для определения характеристик реализации параллельного алгоритма расчета газодинамического течения в коде была проведена серия расчетов задачи из данной работы с использованием разного числа вычислительных ядер. Расчеты проводились с использованием одной из неструктурированных сеток, с числом ячеек около 6 миллионов.

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

ускорениям реализации алгоритма. Пунктирная линия отвечает предельно возможной с теоретической точки зрения зависимости ускорения от числа ядер (максимальному ускорению).

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

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

У [мм]

60

50

40

30

20

10

Неустойчивость Кельвина - * Гкльмгольца

У [ мм]

60

50

40

30

20

10

30 х [мм]

а)

б)

Рис. 4. Пространственные распределения градиента плотности (верхняя половина рисунка) и температуры (нижняя половина рисунка) в момент времени 41 мкс, полученные в расчетах на сетке, состоящей из: a — 28 млн ячеек; б — 1.2 млн ячеек

]

I

0

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

Наличие в канале вертикальных стенок усиливает газодинамические параметры среды в области эллиптических участков. Усиление параметров среды оказывает влияние на картину воспламенения смеси и время воспламенения смеси. Воспламенение модельной водородно-кислородной смеси происходит в момент времени около 20 мкс. Далее, происходит распространение волн горения в канале.

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

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

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

Была также проведена серия расчетов на сетках с другим числом ячеек. Проведены расчеты на сетках, состоящих из: 1.2 млн, 6 млн, 18 млн ячеек. Пространственные распределения газодинамических параметров в момент времени 41 мкс для сетки, состоящей из 1.2 млн ячеек, приведены на рис.4б. Картины инициирования детонации на разных расчетных сетках соотносятся друг с другом. В частности, момент инициирования детонации одинаковый во всех расчетах и равен приблизительно 41 мкс. Однако, детальное разрешение таких элементов и структур газодинамического течения, как неустойчивость Кельвина-Гельмгольца, контактных поверхностей, областей горения смеси наблюдается на самых детальных расчетных сетках.

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

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

Декомпозиция области проводится с использованием функций библиотеки METIS. Обмен данными между вычислительными ядрами осуществляется с применением библиотеки MPI. Рассмотрен вопрос ускорения реализации параллельного вычислительного алгоритма.

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

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

При проведении исследований использовались суперкомпьютеры МСЦ РАН МВС-10П и МВС-10П ОП.

Работа выполнена в рамках госзадания ИАП РАН.

ЛИТЕРАТУРА

[1] Vasil'ev A.A. Cellular structures of a multifront detonation wave and initiation (review). Combust. Explos. Shock Waves, 2005, vol. 5, iss. 11, pp. 1-20.

[2] Ivanov V.S., Frolov S.M., Zangiev A.E., Zvegintsev V.I., Shamshin I.O. Hydrogen fueled detonation ramjet:Conceptual design and test fires at Mach 1.5 and 2.0. Aerospace Science and Technology, 2021, no. 109, pp. 1-12.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

[3] Lopato A.I. Numerical simulation of shock-to-detonation transitions using one-stage and detailed chemical kinetics mechanism. In: Favorskaya, M.N., Favorskaya, A.V., Petrov, I.B., Lakhmi, C.J. (eds.) Smart Modeling for Engineering Systems, 2021, pp. 79-88.

[4] Hirsch C. Numerical computation of internal and external flows: the fundamentals of computational fluid dynamics. Amsterdam, Elsevier, 2007, 680 p.

[5] Schultz E., Shepherd J. Validation of detailed reaction mechanisms for detonation simulation. Caltech Explosion Dynamics Lab, 2000, no. FM99-5.

[6] Salome. The open space integration platform for numerical simulation. https://www.salome-platform.org. Last accessed 17 June 2023.

[7] Liou M.S., Steffen C.J.Jr. A new flux splitting scheme. Journal of Computational Physics, 1993, no. 107, pp. 23-39.

[8] Besnard J., Malony A., Shende S., Perache M., Carribault P., Jaeger J. An MPI halo-cell implementation for zero-copy abstraction. EuroMPI 15: Proceedings of the 22nd European MPI Users' Group Meeting, 2015, pp. 1-9.

[9] Savin G.I., Shabanov B.M., Telegin P.N., Baranov A.V. Joint Supercomputer Center of the Russian Academy of Sciences: Present and Future. Lobachevskii Journal of Mathematics, 2019, no. 40, pp. 1853-1862.

[10] Visit. https://https://visit-dav.github.io/visit-website. Last accessed 17 June 2023.

[11] Utkin P.S., Lopato A.I., Vasil'ev A.A. Mechanisms of detonation initiation in multi-focusing systems. Shock Waves, 2020, vol. 30, iss. 4, pp. 741-753.

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

Ссылку на эту статью просим оформлять следующим образом:

Лопато А.И. Математическое моделирование инициирования детонации в канале с профилированным торцом с использованием параллельных вычислений.

Математическое моделирование и численные методы, 2023, № 4, с. 15-26.

Лопато Александр Игоревич — канд. физ.-мат. наук, с.н.с. ИАП РАН. e-mail: lopato2008@mail.ru.

Mathematical Modeling of Detonation Initiation in the Channel with the Profiled End Using Parallel Computations

The work is devoted to numerical studies of detonation initiation in a gas mixture in a rectangular channel with a profiled end. Detonation is initiated as a result of the interaction of shock waves which are formed as a result of the reflection of an incident shock wave of relatively low intensity from the end of the channel. The mathematical model includes the system of gas dynamics equations supplemented by Arrhenius kinetics for a model hydrogen-oxygen mixture with tabular kinetic parameters corresponding to the operating range ofpressures and temperatures of the mixture. Numerical calculations are carried out using the finite volume method. The construction of computational grids consisting of triangular cells is carried out using the free software SALOME. The numerical algorithm is parallelized by the computational domain decomposition method using the METIS library. The exchange of grid functions between computational cores is carried out using the functions of the MPI library. The problem of acceleration of the parallel algorithm realized in the code is considered in comparison with the case of the linear dependence of the number of computational cores. A number of calculations were carried out using a different number of triangular cells and a comparison of patterns of detonation initiation was carried out. The performed calculations show that the detonation initiation time is approximately the same in all computations. The main difference in detonation patterns is associated with gas flow and Physical and chemical reactions in the mixture.

Keywords: mathematical modeling, detonation initiation, gas mixture, unstructured grids, parallel calculations

REFERENCES