CC BY
5УДК 531.745:533.27:519.6
ВЫЧИСЛИТЕЛЬНОЕ МОДЕЛИРОВАНИЕ ЗАДАЧ ГОРЕНИЯ ГРЕМУЧИХ ГАЗОВЫХ
СМЕСЕЙ
В. Б. Бетелин1,2, Н. Н. Смирнов1,2, В. Ф. Никитин1,2, М. Н. Смирнова2,3, Л. И. Стамов1,2,
В. В. Тюренкова1,2*
1 Научно-исследовательский институт системных исследований 2 Московский государственный университет имени М. В. Ломоносова 3 Санкт-Петербургский государственный политехнический университет
* tyurenkova.v.v@yandex.ru
Целью настоящей работы является исследование особенностей вычислительного моделирования на многопроцессорных ЭВМ универсальной и гибридной архитектуры задач горения и детонации газовых смесей, распространяющихся в нестационарном режиме. Проводится исследование влияния разрешения сетки на результат моделирования. Также проводится сравнение производительности графической платы и универсальных процессоров, подключаемых параллельно, при решении задач инициирования детонации.
Ключевые слова: горение, детонация, химическая кинетика, газовая динамика, параллельные вычисления, графические процессоры.
NUMERICAL SIMULATIONS FOR SOLVING PROBLEMS OF COMBUSTION IN GASES
V. B. Betelin12, N. N. Smirnov12, V. F. Nikitin12, M. N. Smirnova23, L. I. Stamov12,
V. V. Tyurenkova12*
1 System Analysis Research Institute, Russian Academy of Sciences
2 Moscow State University 3 St. Petersburg State Polytechnic University
* tyurenkova.v.v@yandex.ru
The paper studies the features of numerical simulation and parallel computing as applied to combustion and detonation problems in gases for unsteady modes. The effect of grid resolution on the simulation results has been studied. A comparison of the GPU and CPUs performace operated in paralle, for solving detonation initiation problems is provided.
Keywords: combustion, detonation, chemical kinetics, gas dynamics, parallel computing, GPU.
Введение
Горение — один из основополагающих процессов в деятельности человека, который был освоен и используется в течении тысячелетий. Научные же исследования процессов горения начались относительно недавно — в середине 19-го века. Наиболее интригующим в науке о горении является тот факт, что при одинаковом химизме реакций и энерговыделении процесс может протекать в режимах принципиально отличающихся друг от друга, что необходимо учитывать при управлении техническими или природными процессами, включающими горение как одну из стадий.
В 1881 году французские физики Малляр и Ле Шателье [1] и, независимо от них, Бертело и Вьей [2], проводя опыты по распространению пламени, обнаружили, что в обычных условиях пламя в трубе, заполненной горючей газовой смесью, распространяется с небольшой скоростью порядка нескольких метров в секунду. Но при некоторых обстоятельствах медленный процесс горения внезапно переходит в очень быстрый процесс, распространяющийся со скоростью 2 000-3 000 м/с. Этот быстрый процесс горения был назван «фальшивым горением», или детонацией от фр. detonner — фальшивить, звучать не в тон. Чтобы отличать от этого сверхзвукового режима распространения, дозвуковые режимы горения назвали «дефлаграция». Факт наличия двух скоростей горения, различающихся на три порядка, требовал теоретического объяснения, которое впервые было дано профессором Московского университета В. А. Михельсоном в 1893 г. в его работе «О нормальной скорости воспламенения
гремучих газовых смесей», опубликованной в Учёных записках Московского Университета [3]. Опираясь на появившиеся к тому времени работы Рэнкина [4] и Гюгонио [5], В. А. Михельсон объяснил, что механизмом распространения горения при детонации является не теплопередача, а «нагревание до точки воспламенения..., благодаря чрезвычайно быстрому адиабатному сжатию» в ударных волнах. Таким образом, была начата разработка классической теории детонации, получившей дальнейшее развитие в работах Чепмена (1899) [6] и Жуге (1905) [7]. Это явилось важным этапом в развитии в нашей стране школы наук о горении, заложенной ещё М. В. Ломоносовым в XVIII в., и занявшей лидирующее положение в мире во второй половине ХХ столетия. Несмотря на бесспорный приоритет В. А. Михельсона в разработке теоретической модели детонации, детонационная волна, распространяющаяся в самоподдерживающемся режиме, получила название детонации Чепмена—Жуге, которое и используется по сей день.
В настоящее время предсказательное моделирование сложных трёхмерных процессов распространения детонационных волн возможно только с применением современной мощной вычислительной техники. При решении задач газовой динамики с химическими превращениями определяющая система уравнений может интегрироваться методами разделения по физическим процессам. В этом случае анализ расходования процессорного времени на интегрирование частей системы, связанных с конвективной газовой динамикой, процессами переноса (теплопроводностью и диффузией) и химической кинетикой, убедительно показывает, что на интегрирование системы кинетических уравнений, в зависимости от выбранной модели, тратится от 60 до 90 % времени расчёта. Важно отметить, что при этом расчёт стационарных режимов распространения фронта детонации или горения занимает существенно меньше времени, чем нестационарный процесс перехода горения в детонацию. Поэтому для повышения эффективности расчёта задач газодинамики горения необходимо развивать как новые алгоритмические, так и аппаратные подходы, ориентированные на решение данного класса задач. В предлагаемой работе приведены результаты модельной трёхмерной задачи газодинамики горения и детонации. Проводится изучение ускорения исполнения расчётных программ при использовании различного числа параллельных исполнителей на различных системах и с различными алгоритмами распараллеливания.
Постановка модельной трёхмерной задачи газодинамики горения и детонации
Рассматривается замкнутая область в виде прямоугольного параллелепипеда размерами Lx х Ly х Lz, заполненная однородно перемешанной смесью водорода, кислорода и азота с заданными соотношениями молярных плотностей [H2] : [O2] : [А2] при заданной начальной температуре Тщ и давлении Рщ. Стенки считаем адиабатическими и некаталитическими; адсорбция и эмиссия газов на стенках отсутствуют. Газы считаем совершенными (но не обязательно калорически совершенными).
Начало координат поместим в угол рассматриваемой области, оси направим параллельно ее сторонам. В начальный момент времени t = 0 в шаре радиусом rign с центром в точке Xign = {Xign,yign,Zign} начинается выделение энергии и массы из постороннего по отношению к рассматриваемой системе источника. Выделение массы и энергии происходит в течение достаточно короткого времени tign. Подобный процесс может моделировать взрыв детонатора из твёрдого ВВ (без конкретизации процесса в самом детонаторе). Параметры зажигания выбираются таким образом, чтобы привести к прямому инициированию детонации в газовой смеси для заданных значений начальной температуры и давления.
Математическая модель
Уравнения. Для расчёта многокомпонентной газовой динамики с химическими превращениями, без учёта явлений переноса и с внешними источниками массы и энергии используем следующие уравнения:
я
(1)
дри д , ч
"5Г + ^ (pkj =ш k + Mи;
д1 дх1 дрщ д , ч др
"¿г + + = 0; (2)
1(Е+4) + щ{(Е+р+р,2)и) = 0+^ (3)
В системе (1)-(3) индекс £ пробегает значения \,...,Ыс (перечень компонент), индексы I,/ — значения 1,2,3 (перечень координат); по повторяющимся индексам ведётся суммирование. Всего в системе
дифференциальных уравнений (1)-(3) N +4 уравнения.
Алгебраические соотношения. Дифференциальные уравнения (1)-(3) дополняются тремя алгебраическими соотношениями, а также алгебраическими выражениями для химических источников и для внешних источников массы и энергии. Первые три алгебраических соотношения выглядят как:
Мс к=1
Мс
Р =
Рк
Мс
к=1
к=1
Е = ^Ц^к {Т )_1) •
(4)
Выражениями (4) определяются, соответственно, плотность газовой смеси р, давление р и внутренняя энергия Е единицы объёма смеси. Обозначено: Яо — универсальная газовая постоянная, Щ — молярные массы компонент, а функции температуры Нк(Т) — безразмерная энтальпия компонент, включающая в себя энтальпию формирования компонент при заданной температуре Тг^ («химическую энергию»). Подробно концепция безразмерных термодинамических данных описана в [8]. В алгебраических соотношениях здесь и в дальнейшем часто удобно пользоваться не парциальной плотностью компонент рк, а молярной плотностью (иногда, особенно в химической литературе, называемой «концентрацией») Хк:
Хк =р к/Щ. (5)
Источник массы Шк связан с химическими взаимодействиями присутствующих в системе газов. Источник массы Мк — выход компонент из внешнего источника. Источник в уравнении энергии — это суммарный приход энергии вместе с внешней массой Мк; источник 0 в уравнении энергии — дополнительный приход энергии из внешнего источника. Источник тепловой энергии от химических взаимодействий внутри системы в настоящей модели не участвует, поскольку энергия Е включает в себя и химическую энергию компонент. Внешние источники не связаны непосредственно с химическими взаимодействиями внутри системы; это внешний подвод энергии и массы, служащий для зажигания или для возбуждения движения в газовой смеси. В данной постановке задачи они явно зависят от времени и локализации.
Химические источники Шк для большинства систем выглядят сложно; их можно представить зависящими от температуры Т и набора молярных плотностей X = {Хк}; сумма этих источников вследствие закона сохранения массы в химических взаимодействиях равна нулю:
Шк = Щ Ш к(Т ,Х),
Мс
ТсШ к = 0.
к=1
(6)
Имеются и более строгие законы химических взаимодействий (закон сохранения массы каждого элемента), которые учитываются в кинетическом механизме и которые могут учитываться в численной реализации модели для упрощения расчётов и повышения точности. Общий вид химических источников довольно сложен и включает члены, нелинейные по всем аргументам; типичным выражением является
Шк = Vг,кшг,
шг = Мг (X)
кг,г(Мг,Т) Д Х^' - кК,г(Мг,Т) ДХ^'
(7)
где ьоГ — скорость (интенсивность) реакции г, 1>г,к — алгебраический стехиометрический коэффициент компонента к в реакции г, Мг — коэффициент влияния третьих (не изменяемых) компонент в реакции г, который равен 1 при отсутствии такого влияния, кг,г — коэффициент скорости прямой реакции, обычно зависящий лишь от температуры, но для некоторых («выпадающих») реакций также от Мг , кв,г — коэффициент скорости обратной реакции, аг,к — степени компонент в прямой реакции (обычно, хотя не всегда, ненулевые степени только у входящих компонент), г,к — степени компонент в обратной реакции.
Зависимости Нк (Т) обычно реализуются как полиномы с различными коэффициентами для различных компонент. Расчёт этих компонент на основе табличных данных о теплоёмкости при различной температуре, и данных об энтальпии формирования — предмет термодинамического моделирования газовой смеси. Мы брали данные из [9]; формат этих данных подробно описан в руководствах [7] и [10].
Преобразование уравнений массы компонент. Для практических расчётов уравнение (1) может использоваться как непосредственно, так и быть преобразовано к следующим уравнениям:
% + ¿И) = М = Е м (8)
1 ъ
дУи_шдУл ш ъ + Мъ Здесь обозначено: Уъ — массовые доли компонент:
т + идХъ = = 5ъ- (9)
V Ръ Ръ ХъЩ (10)
Уъ = _ = ^- = ^ V ш ■ (10)
Преобразование уравнения (1) в систему (8)-(9) приводит к тому, что в системе становится больше уравнений, чем алгебраически независимых переменных. Возможны следующие решения этой ситуации (не считая использования изначальной подсистемы уравнений (1)):
1. Замена одного из уравнений (9) (по выбору, с номером т), алгебраическим соотношением, следующим из определения (10):
Ут =1-13 Уъ ■ (11)
ъ=т
2. Расчёт на каждом временном шаге переопределённой системы так, как если бы переменные Уъ были независимы. После каждого шага проводится коррекция, обеспечивающая выполнение условий (10). Пусть Уъ — значения массовых долей до коррекции, Уъ — после. Сохраняющая пропорции компонент и гарантирующая от смены знака коррекция выглядит как
У (12)
3. Расчёт на каждом временном шаге переопределённой системы так, как будто переменные Уъ являются независимыми, но без какой-либо коррекции, по крайней мере, на каждом шаге. Вместо этого величина дисбаланса
£у =
1
ъ
(13)
выводится как одна из мер погрешности расчётной схемы. В случае превышения этой величиной заданного значения могут приниматься решения по коррекции перед дальнейшим расчётом или по пределам пригодности численной схемы.
Источники массы и энергии. Указанная нами постановка задачи приводит к тому, что источники массы Мъ и энергии Шм, Ш могут быть представлены как:
щ = АМъ [ 1> \х ~ хтп\< ^и,0< * < иёп; (14)
Г 1, X - хк
\ 0, е^е.
\ 1, |х Xig
\ 0, е^е.
Шм = £ Мъеъ (Гч) = ЯоГ^Мъ (4(Т ) - 1) ■ (16)
ш ДШ Г 1,
Iх xign\ ^ rign,0 * tign; (15)
tignQign 1 0,
ъ
Здесь обозначено: АМъ — общая масса компонента ъ из внешнего источника, АШ — общая дополнительная энергия из внешнего источника. Энергия, притекающая вместе с массой (16) определяется из условия, что до смешения с газовыми компонентами системы температура газа равна Т . Объём, в котором выделяется энергия — это объём шара
Qign = 3 (17)
при условии, что центр зажигания Х(ёп лежит внутри рабочей области на расстоянии не менее г щп от ближайшей стенки (в противном случае следует провести необходимую коррекцию этого объёма, формулы которой зависят от этого расстояния; в настоящей работе этот случай не рассматривается).
Граничные условия. Условие на твёрдой стенке соответствует в принятой нами модели условию непротекания:
ип = 0 при хе да, (18)
где да — граница расчётной области, п — вектор нормали к границе в заданной точке (для условия (18) не важно, является нормаль внешней или внутренней). Условие (18) является достаточным, но для численной реализации модели может потребоваться следующее из него и уравнений (2) условие нулевой производной давления по нормали к стенке:
п- Ур = 0 при хеШ. (19)
Начальные условия. Для того, чтобы решить систему (1)-(3), необходимо знать состояние всех алгебраически независимых переменных в момент времени t = 0. Во всей расчётной области в этот момент задаются температура Тп и давление Рп. Начальная скорость полагается равной нулю. Начальные концентрации компонент задаются через их ненормализованные молярные доли (что для совершенных газов эквивалентно объёмным долям), которые обозначим Сщ,к, тем самым соотношение молярных плотностей [Н2] : [О2] : [N2] представляется как С ¡щ,н2 : Сшщ : С ¡щ,м2 (вместо номеров компонент удобно использовать их химические формулы как литералы). Пересчёт на парциальные плотности и на массовые доли проводится с учётом начального давления и температуры по формулам (различные варианты могут использоваться в различных реализациях):
СкЩк рт ,, РШСкЩк
р = , Рк = рУк = „ „ „ „ . (20)
к Е'СЩГ ЯоТш^кУк/Щк' ,к ' к ЯоТшЕ'С''
Расчёт параметров детонации Чепмена—Жуге, химически равновесного состава и параметров в пике Неймана
Параметры детонации Чепмена—Жуге и пика Неймана, который движется по теории (Зельдовича—Неймана—Дёринга) впереди зоны химических превращений, рассчитываются на основе теории скачков в газовой динамике, с привлечением только термодинамики компонент газовой смеси, и без привлечения данных о какой-либо химической кинетике; предполагается лишь само наличие химических превращений, переводящее систему в состояние химического равновесия.
Обозначим параметры с верхним индексом 0 как начальные, конечные параметры после ударного перехода (с возможным изменением состава смеси) — без верхнего индекса. Нижним индексом будем обозначать номера компонент и (при расчёте равновесного состояния) химических элементов.
Пусть волна детонации, представляющая сильный разрыв, движется по неподвижной среде со скоростью й в положительном направлении оси х. Сохранение массы, импульса и полной энтальпии среды при переходе через сильный разрыв приводят к следующим соотношениям в системе, связанной с разрывом:
- р0й = ри; (21)
р0й2 + р0 = р V2 + р; (22)
Н° + 1 й2= Н + 2 V2. (23)
Давление и энтальпия алгебраически зависят от таких параметров газа, как плотность , температура Т и молярные концентрации единицы массы Мк, 1 < к < N. Очевидно, уравнения (21)-(23) не дают однозначного ответа о параметрах за волной детонации без дополнительных условий, поскольку состав газа там меняется.
Дополнительным условием является химическое равновесие за волной детонации.
Изложим алгоритм, позволяющий найти скорость волны детонации и параметры позади неё (конкретно, позади зоны химической реакции). Алгоритм основан на минимизации скорости волны для заданного значения плотности за волной р. Алгоритм изложен в отчёте вАЬСТТ [11]. Заметим,
что алгоритм не претендует на нахождение «точных» условий позади волны детонации, равно как и точного значения скорости — дело в том, что установление химического равновесия предполагается по крайней мере до вхождения газовой смеси за волной детонации на звуковую линию, то есть достаточно быстро. Это верно для многих реакций, сопровождающих горение, но не для всех.
Первая часть алгоритма: расчёт О(р).
1. Ввести условия перед волной: р0, Т0, №; определить р° и Я0; ввести толерантности ец, г л и инкременты Ад, А т ; ввести р.
2. Ввести начальное приближение О = 2 000, Т = 2 000.
3. Найти равновесное состояние (ЛЬ) по начальному (Л0) для заданной Т в условиях постоянной плотности .
4. Вычислить (Н, Р)(О, Т):
О0 ( V
V = О^, Н = [н(Т Л) + -
И
Я0+?)
Р=
(р(Т,РЛ) + (*?)- (р° + р°О2) ■ (24)
5. Вычислить (Нт,Рт)(Т + Ат,О) и (Нд,Рд)(Т,О + Ад), затем пересчитать в компоненты якобиана
Нт =
Нт — Н
Нд =
Нд^Н
Рт =
РТ^Р
Рд =
Рв^Р
(25)
Дт ' д Ад ' 1 Ат ' д Ад
6. Решить линейную систему уравнений 2 х 2 (например, методом Крамера) относительно поправок 8т , 8д:
Нт$т + Нд8д = — Н,
РтЬ + Р^п = - Р■
Решение системы (26) выглядит как:
г =
Нт Нд Рт Рд
8Т = -г-
Н Нд Р Рд
бд = - г-
Нт Н Рт Р
(26)
(27)
7. Рассчитать поправку Т = Т + 8т, д = д + 8д. Определить соблюдение критерия сходимости
|<Г| < + £л,
< екд +
(28)
Если критерий (28) не соблюдается, идти на пункт 3.
8. Рассчитать равновесное состояние (Л) по начальному (№) для заданной Т в условиях постоянной плотности р. Найти р(Т,р,Л]) и Я(Т,Л).
Вторая часть алгоритма — минимизация д(р). Основана на методе золотого сечения (в отличие от громоздкого метода, предложенного в [11]). Эта часть включает в себя предыдущий алгоритм.
1. Ввести Хтщ(= 1^5) и Хтах(— 2^0), а также начальные данные перед волной и прочие исходные данные для алгоритма, вычисляющего д(р), кроме самой этой плотности.
2. Рассчитать дтт = д(р = р0Хтт), дтах = д^Хтах).
3. Рассчитать
= д(р0- (Х\ = 0■618Xmin + 0■382Xmax)); д2 = д^0- (Х2 = 0■382Xmin + 0^618Хтах))^
Если условие > Оm^n не выполнено, требуется уменьшить Хтщ. Если условие д2 > дтах не выполнено, требуется увеличить Хтах. После чего переходить к пункту 2, но если эти величины вышли за некоторые заданные пределы — это фатальная ошибка. Нельзя уменьшать Хтщ до единицы; относительно Хтах термодинамических ограничений нет.
4. Если < д2, то присвоить Хтах := Х2, дтах := д2, Х2 := Хь д2 := и рассчитать (О\,X\) по формуле (29). Иначе присвоить Хтщ := Хь Оmin := Х1 := Х2, := д2 и рассчитать ^2^2) по формуле (29).
5. Проверить критерий сходимости: если
(29)
Х2 < £п,
(30)
то вычислить как окончательный ответ д(р = р°(Х1 + Х2)/2) с финальным расчётом состояния позади волны. Иначе необходим переход к пункту 4.
Финальное состояние позади волны детонации рассчитывается следующим образом.
1. Плотность: р = p°(Xi + X2)/2;
2. Скорость детонации: D = D(p);
3. Прочие параметры (концентрации Nk, температура T и давление p) определяются в рамках алгоритма D(p).
Расчёт химического равновесия в условиях постоянной плотности и температуры. Для
того, чтобы рассчитать равновесное состояние, будем пользоваться алгоритмом NASA CEA, опубликованным в отчёте NASA [12] в нескольких вариантах. Нам требуется расчёт равновесного состояния для заданной температуры при постоянной плотности. Для этого минимизируется энергия Гельмгольца смеси f (T,Nj) с одновременным условием сохранения элементов в химических реакциях:
nc nc
bq = Y, nqkNk = J2 nqkNk = b°q, q = 1,... ,Ne . (31)
k=i k=i
Здесь nqk — матрица элементного состава компонент. Энергия Гельмгольца рассчитывается следующим образом:
f (T,Nj) = ¿Nkfk; fk = Ek(T) - TSk(T)+ RGT ln ^^fL. (32)
k=i PB
Вводятся множители Лагранжа по числу элементов Aq и составляется целевая функция минимизации:
ne nc ne / nc \
F = f + E Xq(bq - bq0) = £ fkNk + ^ Aq X>qkNk - b; . (33)
q=1 k=1 q=1 \k=1 )
Варьируя функцию (33), получим следующие условия равновесия:
ne nc
fk + £ Aqnqk =0; £ nqkNk - b°q = 0. (34)
q=1 k=1
После введения вариаций Д ln Nk, обозначения hj = — Aj/ (RgT), рассмотрения первого члена разложения в ряд Тейлора и деления на RgT, следуют уравнения относительно вариаций:
ne f
Aln Щ-^п^щ = - RT; (35)
j=1
nc
J^NknqkA ln Nk = b0 bq. (36)
k=1
Из уравнений (35)—(36) удобно выразить поправки логарифма концентраций через множители Лагранжа, сократив линейную систему до числа уравнений Ne :
ne f
Aln Nk = £ njk*j-RT; (37)
j=1 RgT
Ne
L
p=1
Nc
Nknpknqk
k=1
^P = bq + ^ZN^nqk. (38)
k=1
На каждой итерации алгоритма следует решать систему (38), затем из (37) выразить поправки концентраций, а затем умножить прежнее значение концентраций на эти поправки.
Равновесное состояние смеси весьма часто существенно отличается от состояния исходного. Поэтому для сходимости итерационного процесса имеет значение выбор начального приближения (не
путать с исходным состоянием!). В работе [12] предлагается в качестве начального приближения (в случае заданной температуры и минимизации энергии Гельмгольца) брать Лъ = 0,1/Лс.
Итерации, полученные с помощью метода, приведённого выше, весьма часто дают слишком сильные поправки на переменные. Это происходит в основном в двух ситуациях. Во-первых, в начале расчёта, когда приближения весьма далеки от равновесного состояния. Во-вторых, ближе к концу расчёта, если велики поправки на концентрации веществ, присутствующих в малых количествах. Далее молярная доля каждой из компонент сравнивается со значением е и 10-8. Если Лъ/Л = ^¡Л > е, то коррекция Лъ ограничена величиной е2 и 7,4. Тем самым фактор контроля для этих фракций вычисляется как:
А1—2/ шах{|Д1пЩ} . (39)
В том случае, когда Лъ /Л < е и при этом А 1п Лъ > 0, вводится ещё одно ограничение, определяемое фактором
1п Л /Л- 921
Л2 = тт
Д1п Л
(40)
Итоговый фактор коррекции определяется как:
А = шт{1,АЬА2} . (41)
Новые приближения итерационного процесса с учётом фактора коррекции вычисляются как:
Лъ(Ж) = Л«ехр(АД1п Лъ) ■ (42)
Условиями прекращения итерационного процесса являются следующие:
Л|Д 1пЛъ\ < 0,5- 10-5;
Л (43)
< 10^6-шахЬ°° для > 10_6, д = 1,...,ЛЕ,
Ц- Ьд
д
где шахд Ь^ — максимальное значение суммы (31) для химических элементов. Если процесс не сходится за 50 итераций, он аварийно прекращается. Как правило, при старте от случайного начального состояния требуется от 8 до 20 итераций, в большинстве прочих случаев — меньше.
Расчёт условий в пике Неймана. Пик Неймана — головная часть детонационной волны; на первичном сильном разрыве химические процессы ещё не инициированы, и состав смеси остаётся тем же, что и перед волной. Волна же движется с заданной скоростью д, которая перед расчётом пика Неймана получена при условии, что позади всего детонационного комплекса устанавливается химическое равновесие и что реакции достаточно интенсивны, чтобы закончиться до того, как газ позади волны перейдёт звуковую линию.
Расчёт пика Неймана содержит вспомогательную процедуру оценки дисбаланса энтальпии как функции плотности в пике Неймана и основную процедуру — обнуление этого дисбаланса методом Ньютона.
Расчёт дисбаланса энтальпии как функции плотности. Пусть заданы начальные условия р0, Т0, № и скорость волны д. Пусть при этом задано значение плотности на пике Неймана р. Из (21) следует
0
V = — д. (44)
Р
Подставляя (44) в (22), находим давление в пике Неймана как функцию :
р = pv2. (45)
По давлению, плотности и начальному составу находим температуру газа в пике Неймана как функцию плотности:
Т =-р-й ■ (46)
р*оТ,ъЛъ0
Дисбалансом энтальпии назовём дисбаланс уравнения (23); он рассчитывается на основе (44)-(46) как:
ДЯ(р) = h(T,Nj0) - h(T°,N0) + 1 (v?- D2) . (47)
Расчёт пика Неймана. Этот расчёт проводится после расчёта условий детонации Чепмена— Жуге; скорость детонации считается известной. В качестве начального приближения плотности за скачком берём плотность смеси за волной детонации Чепмена—Жуге. Далее запускается итерационный процесс метода Ньютона; на каждой итерации последовательно рассчитываем:
ДЯ1 =ДЯ(р + ер); ДЯ0 = ДЯ(р); Ар = - АЯ° ; р + = Др. (48)
¿лп\ — ¿лп°
Здесь производная от дисбаланса энтальпии вычисляется с помощью конечной разности, для чего приращение плотности ер должно быть достаточно мало. Условием прекращения итераций (48) можно взять, например,
|Др| <ер. (49)
После этого следует рассчитать условия в пике Неймана по формулам (44)-(47); сам дисбаланс энтальпии вычислять уже не обязательно (можно вычислить как меру погрешности расчёта).
Тестовая задача расчёта детонации в газовой смеси водорода, кислорода и азота
В качестве теста использовалась задача расчёта детонации в газовой смеси водорода, кислорода и азота (соотношение кислорода к азоту 1 : 3,75, что приблизительно соответствует таковому в воздухе). Зажигание смеси осуществлялось одновременным вводом в смесь извне за короткое время массы и энергии, что моделирует использование твёрдого детонатора; параметры зажигания были выбраны так, чтобы осуществлялось прямое инициирование детонации газовой смеси.
Геометрия системы представляет собой прямоугольный параллелепипед размером (50 х 5 х 5) см. Стенки замкнуты; начальный состав смеси [Я2] : [О2] : [N2] =2:1: 3.75. Начальное давление Pini = 1 бар, начальная температура Tini = 300 К.
Зажигание осуществляется подводом энергии и массы в область шарообразной формы радиусом rtgn = 2 мм в течение первых tgn = 2- 10_6 с после начала процесса. Подводимая масса равнялась ДМ = 0,335 мг, дополнительная энергия (кроме энергии, подведённой вместе с массой) ДQ = 134 мДж. Эти параметры могут соответствовать детонатору, например, 249,5 мкг азида свинца, разбавленного 85,5 мкг хлорида натрия [13, 14].
Параметры детонации Чепмена—Жуге для такой смеси: Dqj = 1956,7 м/с, Pcj = 15,31 бар.
Давление в пике Неймана, движущемся по теории Зельдовича—Неймана—Дёринга перед волной детонации: Pnp = 27,11 бар.
Список компонент и кинетический механизм. Водород, кислород и азот были лишь компонентами начальной смеси. При горении помимо окончательного продукта — водяного пара — образуются многочисленные промежуточные продукты (радикалы), которые при достаточно высокой температуре так и остаются в газовой смеси. Использовались следующие компоненты:
{Я2О, ОЯ, Я, O, ЯО2Я2О2; O2, Я2, N2} .
В расчётах используется кинетический механизм горения водорода, без учёта образования окислов азота (эти реакции достаточно медленные, чтобы влиять на детонацию и горение, и их обычно рассчитывают a posteriori). Нами взят за основу механизм Мааса—Варнаца из 19 обратимых реакций; интенсивность обратной реакции рассчитывалась по интенсивности прямой и термодинамическим данным участвующих компонент.
В первом столбце табл. 1 помещён номер реакции, в следующем формула реакции, далее — коэффициент скорости прямой реакции в виде обобщённого закона Аррениуса (конкретное значение, принятое в работе Мааса—Варнаца [15] (1989) и позже перепечатанное в составе более сложного механизма в работе Мааса и Поупа [16], (1992)). В последней работе существенно уточнены данные по механизму первой реакции О2 + Я ОЯ + О, имеющей, как показывают исследования [17], критическое влияние на весь процесс, поскольку представляет собой разветвление цепного механизма по наиболее быстро идущей ветви. К 19 реакциям добавлена ещё одна из механизма [18] (в таблице ниже
номер 20). Анализ и тестирование механизма приведены в [19, 20]. Единицами измерения в табл. 1 являются сантиметр, моль, секунда, градус Кельвина, килоджоуль.
В табл. 1 реакции 1-4 соответствуют обмену с участием «лёгких» радикалов О, Н, ОН; реакции 5-7 соответствуют рекомбинации «лёгких» радикалов; в реакциях 8-13 участвует ещё и «тяжёлый» радикал НО2; в реакциях 14-19 ещё один «тяжёлый» радикал Н2О2. Реакция 20 — ещё одна рекомбинация лёгких радикалов; её учёт несколько ускоряет образование конечного продукта и «выгорание» атомарных кислорода и водорода.
Механизм Мааса—Варнаца обладает ещё одним свойством, упрощающим расчёты: влияние третьего тела в реакциях рекомбинации определяется одними и теми же коэффициентами компонент для всех реакций (коэффициентами Шаперона). В кинетическом механизме в таких реакциях стоит символ «+М». Последняя реакция 20 имеет другие коэффициенты Шаперона [19]. Эти коэффициенты (выписанные лишь для компонент нашей смеси) приведены в табл. 2.
Таблица 1
Кинетический механизм Мааса—Варнаца
№ п/п Реакция Коэффициент скорости
1. О2 + Нт±ОН + О 2.00 • 1014 ^[^70.3/^7]
2. Н2 + О ОН + Н 5.06 • 104 • Т2-67 • exp [-26.3/ЯсТ]
3. Н2 + ОН Н2 О + Н 1.00- 108- Т160 ^[^13.8/^]
4. ОН + ОН Н2 О + О 1.50 • 109- Т LИ•expЬ0.4/ЯGT]
5. Н + Н + Мт±Н2 + м 1.80- 1018 • Т-100
6. О + О + Мт±О2+ м 2.90- 1017 • Т-100
7. Н + ОН + Мт±Н2 О + м 2.20 • 1022 • Т-2 00
8. Н + О2 + М Т± НО2 + м 2.30- 1018 • Т-0 80
9. НО2 + Нт±ОН + ОН 1.50- 1014 • exp [-4.2/ЯеТ]
10. НО2 + Нт±Н2 + О2 2.50- Ю^^^^/ЯсТ]
11. НО2 + Нт±Н2 О + О 3.00- 1013 • exp [^7.2/ЯсТ]
12. НО2 + От±ОН + О2 1.80- 1013 • exp [+1.7/ЯсТ]
13. НО2 + ОН <=> Н2 О + О2 6.00- 1013
14. НО2+ НО2 ^ Н2 О2 + О2 2.50- 1011 ^[+5.2/^]
15. ОН + ОН + Мт±Н2 О2 + м 3.25 • 1022 • Т-2 00
16. Н2 О2 + Н^Н2+ НО2 1.70- 1012 • exp [-15.7/ЯеТ]
17. Н2О2 + Нт±Н2О + ОН 1.00- 1013 • exp [-15.0/ЯоТ]
18. Н2 О2 + От±ОН + НО2 2.80- Ю^^^^/ЯсТ ]
19. Н2 О2 + ОНт±Н2 О + НО2 5.40- 1012 • exp [-4.2/ЯеТ]
20. О + Н + Мт±ОН + м 4.71 • 1018 • Т-1
При переходе к системе СИ, где единицей длины является метр, в кинетическом механизме из табл. 1 предэкспоненциальный множитель следует умножить на 10_6 для реакций обмена и на 10^12 для реакций рекомбинации или диссоциации.
Условия численного расчёта. При одних и тех же параметрах тестировались две различные численные схемы сквозного счёта:
1. Явная схема 2-го порядка точности по пространству и времени на основе МИ8СЬ-интерполяции переменных на грань при расчёте конвективных потоков [22]. Для выбора направления
Таблица 2
Коэффициенты Шаперона в механизме из табл. 1. Для прочих компонент коэффициенты
полагаются равными 1.0
Реакции H2 O2 H2O N2 H2 O2
1-19 1,0 0,4 6,5 0,4 1,0
20 2,5 1,0 12,0 1,0 12,0
интерполяции и для интерполяции давления применялся метод AUSMP [23]. Схема реализована на регулярной сетке из одинаковых элементов (прямоугольных параллелепипедов), состыковывающихся в произвольной, в общем случае, топологии. Текст программы на языке программирования C. Параллельность исполнения реализована методом OpenMP. Тестовый расчёт проведён на компьютерах ACER (4 процессора), для сетки размером 260 х 26 х 26, и на компьютере АПК-1 (144 процессора всего, 48 на одной плате с общей памятью), для сеток 500 х 50 х 50, 1 000 х 100 х 100 и 2 000 х 200 х 200.
2. Явная схема 2-го порядка точности по пространству и 2-го по времени на основе схемы Курга-нова—Тадмора [24, 25]. Использовалась регулярная сетка, состоящая из прямоугольных параллелепипедов. Программа написана на языке Фортран. Использовалась технология параллельного программирования CUDA [21] для работы с графическими сопроцессорами. Тестовые расчёты проводились на графической карте Nvidia Tesla k40 для сеток 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов.
В дальнейшем мы приведём детальные результаты, полученные с применением второй численной схемы, реализованной по технологии параллельного программирования CUDA. Результаты, полученные с применением первой схемы, будут использованы при сравнении работоспособности и быстродействия.
Физико-математический анализ результатов
Задача рассчитывалась в геометрии прямоугольного параллелепипеда размерами 500 х 50 х 50 мм. Применялись сетки размером 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов; также дополнительно было использовано по 3 фиктивные ячейки на каждой границе области для реализации граничных условий. Сетки состояли из одинаковых кубических элементов, длина стороны элемента была 1.92, 1.0 и 0.556 мм для каждой из сеток соответственно. Расчёты проводились на компьютере с двумя 6-ядерными процессорами и графической картой NVidia Tesla K40. Основные вычисления производились на графическом сопроцессоре. Часть расчётов, связанных с определением временного шага, производились на центральном процессоре.
На всех сетках расчёт производился до момента времени 250 мкс. Также производилось тестирование времени исполнения на первых 10 мкс для различных размеров сеток. Распараллеливание производилось с помощью технологии CUDA для работы с графическим сопроцессором.
На рис. 1-3 изображено давление в плоскости Oxy, проходящей через ось симметрии расчётной области, параллельную оси координат Ox. Давление в этой плоскости приведено для моментов времени 50, 100, 150, 200 и 250 мкс от начала инициирования. Ниже на рисунках изображено давление вдоль оси симметрии системы для тех же моментов времени. Рис. 1 соответствует сетке 260 26 26, рис. 2 сетке 500 50 50, рис. 3 — сетке 900 90 90. На графиках изображены также уровни теоретического значения давления Чепмена—Жуге и значения давления в пике Неймана.
Картина развития давления показывает, как после начального инициирования развивается детонационный процесс: образуется лидирующая волна с характерным для детонации профилем, которую догоняет отражённая от стенки ударная волна с большим давлением на фронте. Также наблюдаются многочисленные вторичные ударные волны и волны разрежения, как следствие многократных отражений от стенок. Со временем вторая ударная волна, распространяющаяся по уже прореагировавшей смеси газов, догоняет первую, и далее происходит распространение только одной волны чуть большей амплитуды. На более грубых сетках профили волн размыты сильнее, и их максимумы сглажены. На более подробных сетках заметно всё больше деталей волновой картины. Также заметно, что с увеличением разрешения сетки увеличивается по величине пик волны детонации, который на рассмотренных сетках ещё не достигает теоретического значения пика Неймана.
На рис. 4 изображено развитие поля плотности для моментов времени 50, 100, 150, 200 и 250
Рис. 1. Давление, рассчитанное на сетке 260 х 26 х 26. Вверху — сечение плоскостью Оху для моментов времени 50, 100, 150, 200 и 250 мкс. Внизу — давление на оси Ох по центру системы для этих же моментов времени
мкс в плоскости Оху для сетки 900 х 90 х 90. Картина плотности в целом повторяет картину давления.
На рис. 5 изображено развитие поля температуры в плоскости Оху для тех же моментов времени, что и для плотности на рис. 3 и для такой же расчётной сетки. После прохождения детонационной волны в системе устанавливается (приближённо вследствие динамики) химическое равновесие. Изменения температуры в пределах 2500-2800 К в этой области отслеживают волновую картину: в волнах сжатия температура смеси повышена, в волнах разрежения — понижена.
На рис. 6 изображено развитие поля модуля скорости в плоскости Оху для тех же моментов времени и сетки, что и рассмотренных ранее. Картина абсолютной скорости отслеживает картину давления (лидирующей волны и отражённых волн). Заметно, как после 100 мкс от начала процесса уменьшается величина скорости в отражённых волнах.
Рис. 2. Давление, рассчитанное на сетке 500 х 50 х 50. Вверху — сечение плоскостью Оху для моментов времени 50, 100, 150, 200 и 250 мкс. Внизу — давление на оси Ох по центру системы для этих же моментов времени
Чтобы лучше проследить за отражёнными волнами, следует вывести картину развития составляющей вектора Уу со временем. На рис. 7 изображено развитие поля модуля скорости в плоскости Оху для тех же моментов времени, что для рис. 4-6. Отрицательная величина вектора скорости в данном случае соответствует направлению составляющей вектора скорости ВНИЗ, а положительная ВВЕРХ.
Картина поперечной скорости показывает, как в системе развиваются и с течением времени затухают поперечные волны. Местоположение лидирующей волны детонации заметно на рис. 7 только для моментов 50 и в меньшей степени 100 мкс; в дальнейшем волна детонации становится практически плоской, и поперечная составляющая скорости становится достаточно малой. Заметно также затухание поперечных волн, возбуждённых первичным инициированием, с течением времени.
Рис. 3. Давление, рассчитанное на сетке 900 х 90 х 90. Вверху — сечение плоскостью Оху для моментов времени 50, 100, 150, 200 и 250 мкс. Внизу — давление на оси Ох по центру системы для этих же моментов времени
В используемом нами кинетическом механизме важную роль играет радикал пергидроксил НО2. Дело в том, что эта компонента играет важнейшую роль в первичной диссоциации исходных продуктов Н2 и О2, а также в том, что этот радикал при высокой температуре быстро разлагается. Тем самым концентрация пергидроксила хорошо подходит для отслеживания волн горения водорода в кислороде: на волне горения она на 2 порядка выше, чем в областях химического равновесия. В то же время и в области химического равновесия это вещество присутствует, отслеживая детали волновой картины (равно как температура, плотность и др.)
На рис. 8 изображено развитие поля массовой доли НО2 в плоскости Оху для тех же моментов времени, что для рис. 4-7.
Заметно, что форма лидирующей детонационной волны в плоскости Оху для момента времени
Рис. 4. Развитие плотности в плоскости Оху симметрии системы
Рис. 5. Развитие температуры в плоскости Оху симметрии системы
50 мкс ещё не является плоской. Заметно постепенное уменьшение доли НО2 позади волны детонации; это связано с тем, что в равновесном состоянии она уменьшается с уменьшением давления, последнее же постепенно падает за волной.
Рис. 9-12 посвящены сравнению полей давления, температуры, абсолютной скорости и массовой доли радикала НО2, рассчитанных для момента времени 50 мкс от начала процесса для различных сеток: 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов.
Из рис. 9 заметно, что по мере увеличения разрешимости сетки местоположение волны сме-
Рис. 6. Развитие поля модуля скорости в плоскости Оху симметрии системы
Рис. 7. Развитие поля поперечной составляющей скорости в плоскости Оху симметрии системы
щается влево; это, возможно, связано с увеличением скорости волны или с увеличенной величиной численной диффузии на более грубой сетке (размывание температуры и концентрации радикалов). Видно также, что для сеток 500 50 50 и 900 90 90 это явление уже незаметно, и результаты весьма похожи, и фронт волны практически совпадает, хотя разрешение волновых процессов на более тонкой сетке значительно лучшее, и заметно лучшее выпрямление фронта на более мелкой сетке.
На рис.10 поле температуры в целом отслеживает поле давления.
По изменению поля абсолютной скорости при изменении разрешения сетки хорошо заметно
Рис. 8. Развитие поля массовой доли радикала НО2 в плоскости Оху симметрии системы
Рис. 9. Поле давления в плоскости Оху в момент времени 50 мкс от начала процесса, для сеток 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов
уменьшение численного размывания по мере уменьшения шага сетки по пространству.
На картинах поля молярной доли радикала НО2 хорошо заметно выпрямление фронта волны.
На рис. 13 изображено состояние поля концентрации радикала НО2 для момента времени 100 мкс для тех же сеток, что и на рис. 12.
Можно заметить, что на рис. 13 хорошо прослеживается более высокая скорость волны на грубой сетке и слабое выпрямление фронта. Стоит отметить, что в целом, поведение волны на грубой сетке отличается от более медленного распространения на более мелких сетках и характеризуется более высокой скоростью. В результате чего, «вторая» волна не успевает догнать «первую» до конца расчётной области. Скорее всего, это связано с большой численной диффузией.
На рис. 14 изображено состояние профиля массовой доли всех химических компонентов вдоль центральной оси рабочей области при £ =100 мкс с начала процесса. Использована логарифмическая шкала для концентраций компонент. Рисунок содержит данные по расчётам на сетках 260 26 26, 500 х 50 х 50 и 900 х 90 х 90 элементов.
Рис. 14 показывает, как начальное состояние состава смеси (справа от скачкообразного перехо-
Т ; ( = 5.0306е-05 б
О 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
к. т
Т ; ( = 5.0977е-05 Б
Рис. 10. Поле температуры в плоскости Оху в момент времени 50 мкс от начала процесса, для сеток 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов
V ; ( = 5.0306е-05 б
Рис. 11. Поле модуля скорости в плоскости Оху в момент времени 50 мкс от начала процесса, для сеток 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов
Н02 ; (оо( = 5.0306е-05 б
Рис. 12. Поле молярной доли радикала НО2 (пергидроксила) в плоскости Оху в момент времени 50 мкс от начала процесса, для сеток 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов
да в окрестности 17-20 см от левого конца системы) переходит в состояние химического равновесия (слева от перехода), и как состояние равновесия изменяется за счёт колебаний температуры в поперечных волнах позади лидирующей детонационной волны, и в окрестности начального инициирования. Видно, что при переходе к равновесному состоянию:
- концентрация азота не меняется (практически; малые её отклонения, не заметные на данном масштабе, происходят за счёт различной численной диффузии различных компонент);
Н02 ; (01й = 0.00010044 Б
Рис. 13. Поле молярной доли радикала НО2 (пергидроксила) в плоскости Оху в момент времени 100 мкс от начала процесса, для сеток 260 х 26 х 26, 500 х 50 х 50 и 900 х 90 х 90 элементов
- концентрация кислорода, для данной (стехиометрической) смеси падает при сгорании значительно меньше (на порядок), чем концентрация водорода;
- концентрация лёгкого радикала ОН (гидроксила) практически совпадает в равновесном состоянии с концентрацией кислорода, она самая высокая среди концентраций радикалов, около 1 %;
- концентрация атомарного кислорода в равновесном состоянии на порядок превышает концентрацию атомарного водорода, что, скорее всего, объясняется соответствующим соотношением молекулярных компонент;
- концентрация тяжёлого радикала НО2 имеет резкий максимум на волне детонации, спадая сразу на 2 порядка, а затем спадая постепенно. Причина последнего спада — в изменении давления позади детонационной волны; чем меньше давление, тем меньше концентрация этого радикала в равновесном состоянии. Она в равновесном состоянии на 2 порядка ниже, чем концентрация атомарного водорода;
- концентрация Н2О2 в целом ведёт себя аналогично концентрации НО2, но на порядок ниже по значению.
Окрестность волны детонации на рис. 14 практически не разрешена; на рис. 15 также изображены массовые доли всех компонент, но уже в меньшей области на каждом из графиков в окрестности волны детонации при 100 мкс от начала процесса. Поскольку для различных сеток место волны детонации несколько смещено, на разных графиках этого рисунка изображено состояние концентраций в диапазонах от 19 до 22 см для самой грубой сетки и от 18 до 20 см для более мелких сеток.
На рис. 15 видны детали области перехода от неравновесного химически (начального) к равновесному состоянию в детонационной волне. В области перехода концентрация радикалов ОН, О, Н имеет небольшой максимум (для атомарного водорода и кислорода превышение по сравнению с равновесным на полпорядка, для ОН практически отсутствует. Концентрация тяжёлых радикалов Н2О2 и НО2 имеет резкий максимум — на 1,5-2 порядка выше равновесного. На одну-две ячейки впереди основного подъёма радикалов идут концентрации НО2, Н2О. С уменьшением размера ячейки область перехода пропорционально сужается, это означает, что с уменьшением размера ячейки пропорционально уменьшается зона лидирующей ударной волны, что характерно для схем сквозного счёта.
На следующем рис. 16 изображена динамика скорости волны детонации в зависимости от времени. Положение волны детонации определялось по точке максимального давления на пике детонационной волны; скорость определялась по отношению конечных разностей положения волны и времени.
Получено, что для грубой сетки скорость волны гораздо выше теоретической скорости волны Чепмена—Жуге (1 957 м/с). Для остальных сеток скорость волны немного выше теоретической скорости (на 14-16 м/с). На самом деле, определение скорости связано с большой погрешностью в фиксировании положения волны, а также в том, что волна проходит не целое число ячеек за одинаковые промежутки времени. Чтобы оценить действительную скорость волны для данной численной схемы, необходимо либо взять такой размер ячеек, чтобы волна проходила их одинаковое количество на каж-
260 х 26 х 26
( , = 0.0001 Б
0.05 0.1
0.15 0.2 0.25 0.3 0.35
I , т
Н02 Н202
0.4 0.45
500 х 50 х 50
__
----
/
.............
_^
-N2
-----\ —— ОН Н02 - Н202
900 х 90 х 90
Рис. 14. Профиль массовой доли всех химических компонентов вдоль центральной оси Ох рабочей области при £ =100 мкс с начала процесса
I , - 0.00010044 й
да
Ж -1-1-1-
г /^Ч :
-
;
1_______—" /^Д
; \\ — Н20 . —1—он н -+--о
—'—ь , \ \\ »—Н202
260 х 26 х 26
500 50 50
900 90 90
Рис. 15. Профиль массовой доли всех химических компонентов вдоль центральной оси Ох рабочей области при £ =100 мкс с начала процесса, в окрестности волны детонации
Рис. 16. Динамика скорости детонации, получае- Рис. 17. Динамика значения пика давления в волне мой при различных размерах сетки детонации для сеток с различным разрешением
дом выводимом временном шаге, либо выводить результаты в моменты времени, соответствующие прохождению целого числа ячеек для конкретной сетки. Для остальных случаев будет наблюдаться большая «болтанка» относительно некой средней скорости, близкой к теоретической скорости волны Чепмена—Жуге для данной смеси при выбранных условиях.
На рис. 17 изображена динамика пика давления в волне детонации, полученная для сеток с различным разрешением. Эта величина изображена как функция времени.
Получено, что примерно до времени 90-100 мкс пик давления увеличивается, затем его значение стабилизируется. Итоговая величина выше теоретического значения давления в волне Чепмена— Жуге, что говорит о том, что слабая детонация, вследствие опережающего ударную волну поджигания горючей смеси численной диффузией энергии или концентрации радикалов, не возникает даже на самой грубой сетке. С уменьшением пространственного шага сетки величина пика растёт каждый раз примерно на 1 бар (при уменьшении шага в 2 раза). В то же время полученный численно пик на этих сетках ещё примерно на 8 бар ниже теоретического значения в пике Неймана, также отмеченного на этом рисунке. Эта оценка говорит о том, что достижения близкого к пику Неймана значения при использовании этой сетки следует ожидать при дальнейшем уменьшении шага сетки не менее, чем на порядок. В то же время, для получения необходимой скорости детонации достаточно использовать «среднюю» по количеству элементов сетку.
Анализ эффективности использования технологии CUDA для графических сопроцессоров Далее проводилось сравнение времени вычислений с использованием технологии CUDA для проведения вычислений на графических сопроцессорах. Время вычислений оценивалось на всех рассмотренных сетках в течение первых 10 мкс от начала процесса. Время расчётов представлено в табл. 3.
Таблица 3
Время вычислений для различных размеров сеток на графическом процессоре NVidia Tesla k40
Размер сетки Время расчёта, сек
260x26x26 14,14
500x50x50 173,51
900x90x90 1906,11
Из табл. 3 видно, что компьютерное время расчёта tca\c до определённого момента физического времени процесса нелинейно зависит от количества ячеек расчётной сетки Nceii = Nx х Ny х Nz.
tcalc I "cell I ] , . П
1 1 1 < n < 2;
t0
calc
(NcellX {N0eu)
где верхним индексом ноль помечены значения калибровочных характеристик, которые получены для некоторого базового варианта.
На рис. 18 изображено, в логарифмическом масштабе, время исполнения указанной задачи на АПК-1 (универсальные процессоры) для 0-48 процессов ОрепМР
1 о -1-1-1-1-1-1-1-1-1-
I ю3
_i_I_I_i_i_i_I_i_i_
0 5 10 15 20 25 30 35 40 45 50 OpenMP threads
Рис. 18. Время исполнения тестовой задачи на компьютере АПК-1 для различного числа нитей исполнения при методе OpenMP на сетке 500 х 50 х 50. Нуль нитей исполнения соответствует серийному коду (без распараллеливания), одна нить — использование только одного процесса при задействовании метода OpenMP
Видно, что с ростом числа нитей исполнения на многопроцессорной машине время исполнения задачи в большинстве случаев уменьшается. Исключений два: серийный код (на рисунке время показано при нуле процессов) исполняется быстрее, чем параллельный код при одной нити исполнения. Это естественно, поскольку в серийном коде не тратится компьютерное время и прочие ресурсы на создание и обеспечение работы параллельных процессов. Второе исключение — немонотонный характер этой зависимости при числе процессов больше 40, то есть близкой к предельно возможному числу 48 для одной платы памяти. Причиной немонотонности может быть как характер пересылок, так и вмешательство конкурирующих системных процессов, замедляющих расчёт при дефиците свободных исполнителей.
Сравнение времени вычисления на первых 10 мкс для схемы MUSCL + AUSMP на сетке 500 х 50 х 50 на АПК-1 (рис. 18) и времени вычисления для схемы Курганова—Тадмора на NVidia Tesla показывает, что для того, чтобы добиться той же производительности на АПК-1, требуется не менее 40 нитей исполнения (процессоров), с поправкой, естественно, на различие численных схем. Следует при этом отметить, что по объёму памяти и тем самым по возможности проводить расчёт на более мелкой сетке АПК-1 значительно превосходит использовавшийся графический сопроцессор.
Выводы
Вычислительное моделирование процессов детонации иногда приводит к проявлению схемного эффекта — размыванию концентрации и температуры в окрестности зоны больших градиентов, что ведёт к ускорению распространения фронта до скоростей, превышающих значения Чепмена—Жуге, и, как следствие, понижению давления за фронтом до величин ниже значений Чепмена—Жуге. В результате такого схемного эффекта возникает режим «слабой детонации», который также является одним из возможных режимов распространения волны горения. Проведённые исследования и тестирование разработанных кодов показали, что слабой детонации, обусловленной размыванием концентрации радикала HO2 вперёд лидирующей ударной волной в результате «численной диффузии», не отмечено
даже на самой грубой сетке. Хотя небольшое размывание значений параметров впереди ударной волны в результате схемной «численной диффузии» наблюдается.
В результате схемной «численной диффузии» значения параметров вперёд лидирующей детонационной волной заметно изменяются не более чем в 2-х ячейках перед фронтом для сеток любого исследованного нами размера. Дальше всего вперёд проникают изменения концентрации атомарного водорода и радикалов HO2 и H2O2.
Концентрация тяжёлого радикала HO2 имеет резкий максимум на волне детонации, спадая сразу на 2 порядка, а затем спадая постепенно. Причина последнего спада — в изменении давления позади детонационной волны; чем меньше давление, тем меньше концентрация этого радикала в равновесном состоянии. Она в равновесном состоянии на 2 порядка ниже, чем концентрация атомарного водорода.
Концентрация H2O2 в целом ведёт себя аналогично концентрации HO2, но на порядок ниже по значению.
С уменьшением размера ячейки проявляется всё больше деталей волновой картины позади волны детонации. На грубых сетках многие детали размыты.
С уменьшением размера ячейки рассчитанная скорость детонационной волны всё точнее приближается к значению скорости Чепмена—Жуге.
Применение одного графического процессора NVidia Tesla K40 эквивалентно работе примерно 40 универсальных процессоров в режиме параллельных вычислений по системе OpenMP.
Работа выполнена при поддержке РФФИ грант 14-03-00171.
ЛИТЕРАТУРА
1. Mallard E., Le Chatelier H. L. Sur la vitesse de propagation de l'inflammation dans les melanges explosifso. Compt Rend Acad Sci. Paris, 1993. P. 145-148.
2. Berthelot M., Vieille P. Sur la vitesse de propagation des phenomenes explosifs dans les gaso. Compt Rend Acad Sci. Paris, 1993. P. 18.
3. Михельсон В. А. О нормальной скорости воспламенения гремучих газовых смесей // Уч. записки Императорского Московского ун-та. Отдел физико-математический, 1893. № 10, с. 1-92.
4. Rankine W. J. M. On the thermodynamic theory of waves of finite longitudinal disturbance. Phil Trans Roy Soc. London, 1870. P. 277-288.
5. Hugoniot H. Propagation des mouvements dans les corps et specialement dans les gaz parfaits. Journ Liouville. № 3, P. 477-492; № 4, P. 153-167.
6. Chapman D. L. On the role of explosion in gases. Phil Mag. № 47, 1990.
7. Jouget E. On the propagation of chemical reactions in gases. J. Mathematics, 1905. № 1. P. 347.
8. CHEMKIN. A software package for the analysis of gas-phase chemical and plasma kinetics. CHE-036-1. Chemkin collection release 3.6. Reaction Design, September 2000.
9. Marinov N. M., Pitz W. J., Westbrook C. K., Hori M., Matsunaga N. A. Experimental and Kinetic Calculation of the Promotion Effect of Hydrocarbons on the NO-NO2. Conversion in a Flow Reactor. // Proceedings of the Combustion Institute. Vol. 27, 1998. P. 389-396. (UCRL-JC-129372). UCRL-WEB-204236.
10. Kee R. J., Miller J. A., Jefferson T. H. Chemkin: a general-purpose, problem-independent, transportable Fortran chemical kinetics code package // Sandia National Laboratories Report SAND80-8003. 1980.
11. Browne S., Ziegler J., Shepherd J. E. Numerical Solution Methods for Shock and Detonation Jump Conditions // GALCIT Report FM2006.006 July 2004 — Revised August 29, 2008.
12. Gordon S., and McBride B. J. Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications I. Analysis. NASA RP-1311, October, 1994.
13. Поздняков З. Г., Росси Б. Д. Справочник по промышленным взрывчатым веществам и средствам взрывания. М. : Недра, 1977.
14. Орлова Е. Ю. Химия и технология бризантных взрывчатых веществ : учебник для вузов. 3-е изд., перераб. Л. : Химия, ленинградское отделение, 1981. 312 с.
15. Maas U., Warnatz J. Ignition process in hydrogen-oxygen mixtures // Combustion and Flame, 74, № 1, 1988. P. 53-69.
16. Maas U., Pope S. B. Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space // Combustion and Flame, 88. № 2. P. 239-264.
17. Li J., Zhao Z., Kazakov A., Dryer F. L. An Updated Comprehensive Kinetic Model of Hydrogen Combustion // International Journal of Chemical Kinetics. Vol. 36, 2004. P. 566-575.
18. Hong Z. An improved hydrogen/oxygen mechanism based on shock tube/laser absorption measurements. A dissertation submitted to the department of mechanical engineering and the committee on graduate studies of Stanford University in partial fulfillment of the requirements for the degree of doctor of philosophy. November, 2010.
19. Smirnov N. N, Nikitin V. F. Modeling and simulation of hydrogen combustion in engines // International Journal of Hydrogen Energy Vol. 39. Iss. 2, 2014. P. 1122-1136.
20. Smirnov N. N., Phylippov Yu. G., Nikitin V. F., Silnikov M. V. Modeling of combustion in engines fed by hydrogen. WSEAS Transactions on Fluid Mechanics. Vol. 9, 2014. P. 154-167.
21. NVIDIA CUDA. Programming Guide. 2014, http://developer.nvidia.com/cuda-downloads.
22. Van Leer B. Towards the Ultimate Conservative Difference Scheme. A Second Order Sequel to Godunov's Method, J. Com. Phys., Vol. 32, 1979. P. 101-136.
23. Liou M.-S. A Sequel to AUSM: AUSM+ J Comput Phys, Vol. 129, 1996. P. 364-382.
24. Jiang G.—S. Tadmor E. Nonoscillatory central schemes for multidimensional hyperbolic conservation laws, SIAM J SCI COMPUT, Vol. 19, 1998. P. 1892-1917.
25. Nessyahu H., Tadmor E. Non-oscillatory central differencing for hyperbolic conservation laws. J Comp Phys, Vol. 87, 1990. P. 408-463.
