Тестирование метода адаптивных сеток на расчетах одномерных детонационных волн Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 14, № 3, 2009

Тестирование метода адаптивных сеток

*

на расчетах одномерных детонационных волн

И. А. Бедарев, А. В. Федоров Учреждение Российской академии наук Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]

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

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

Введение

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

Особенностями детонационных течений является то, что в течении присутствуют области с большими и малыми градиентами параметров (плотности, температуры, концентрации компонентов и т.д.). При этом протекают быстрые процессы, связанные с кинетикой химических реакций (время ~ 10-6 с), на фоне сравнительно медленных газодинамических процессов (время ~ 10-3 с). Разрешение мелкомасштабных процессов требует очень подробной сетки. Поэтому для реальных расчетов детонационных течений, очевидно, актуально использование метода адаптивных сеток.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 08-01-92010-ННС-а), гранта Президента РФ МК-2539.2007.1 и АВЦП "Развитие научного потенциала высшей школы" на 2009 год, код проекта 2.1.1/4674.

© ИВТ СО РАН, 2009.

Существует множество способов реализации этого метода [2]. В [3] предложен алгоритм адаптации расчетной сетки к градиенту решения, который позволяет проводить расчеты нестационарных течений для модельного уравнения газовой динамики — уравнения Бюргерса. Этот алгоритм был протестирован в [3] на непрерывных решениях с большими градиентами в некоторых областях пространства. В настоящей работе алгоритм обобщен для решения задач с разрывами параметров течения (ударные и детонационные волны, контактные разрывы) для системы уравнений гиперболического типа, которая описывает газовую динамику неравновесной среды, претерпевающей химические превращения. При этом упор будет сделан в основном на численных аспектах проблемы математического моделирования изучаемого явления.

Физико-математическая модель детонационного течения

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

др + д (ри) = 0 дЬ дх

д (Ри) . д (Ри2 + Р) =0 (1)

дЬ + дх ( )

д (рЕ) + д [(рЕ + р) и]

dt дх

u2

где р — средняя плотность; p — давление; u — скорость; E = e + —--полная энергия;

м м

e = cv(T — T0) + J2 £aeJ (0) — внутренняя энергия; cv = J2 cVa£a — теплоемкость смеси;

a=1 a=1

e{ (0)— энергия образования а-го компонента при T0 = 0 K, e£ (0) = Hf (0). Предполагается калорическая идеальность, т. е. парциальные теплоемкости компонентов cVa = const (выбиралась cVa при T=300 K). Уравнение состояния смеси в соответствии с законом Дальтона имеет вид

м £ м 1м

Р = PRTY,MT ' Y = cPa£a cVa (2)

a=1 M« a=1 / a=1

Здесь Ma, £a — молекулярный вес и относительная массовая концентрация соответствующего компонента. Данная система уравнений должна быть дополнена уравнениями кинетики неравновесных химических реакций, на чем мы и остановимся ниже.

Уравнения кинетики газофазных реакций. Для описания кинетики воспламенения и горения реагирующих компонентов смеси — водорода и его производных — использована детальная кинетическая схема [1], которая учитывает 38 прямых и обратных реакций для восьми компонентов (H2, O2, H2O, OH, O, H, HO2, H2O2). В [1] приведены константы скоростей прямых и обратных реакций. В работе [4] данная кинетическая схема сравнивалась с двумя другими и была верифицирована в широком диапазоне параметров течения по зависимостям времени задержки воспламенения от температуры за фронтом инициирующей ударной волны.

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

Здесь тг — порядок г-й реакции; I — количество реакций; а = 1, ..., ^ — количество компонентов смеси; kfr, кьг — скорости прямых и обратных реакций; иа — стехиомет-рические коэффициенты, штрих относится к продуктам реакции.

Тем самым система законов сохранения (1), дополненная уравнением состояния (2) и уравнениями неравновесной кинетики химических превращений (3), является замкнутой относительно искомых динамических и термодинамических параметров газа, а также концентраций компонентов смеси, возникающей в результате ударно-волнового воздействия.

О выборе кинетической схемы воспламенения и горения водорода в кислороде. Предварительно остановимся на задаче о воспламенении смеси под действием проходящей ударной волны и интерпретации расчетных данных. Важным моментом при этом становится выбор критерия, по которому определяется момент воспламенения ¿¿5П. В качестве критериев могут использоваться максимум концентрации радикала ОН, максимум роста концентрации радикала ОН, максимум роста температуры. Предыдущие наши расчеты [4] показали, что результаты определения ^дп различными способами могут существенно (в 2-3 раза) отличаться друг от друга. Для верификации математической модели по времени задержки воспламенения нами были выбраны экспериментальные данные, в которых явно указывался критерий, по которому определялся момент воспламенения ^дп.

На рис. 1 приведено сравнение расчетных данных по кинетической схеме [1] и экспериментов [5, 6] по зависимости времени задержки воспламенения (tign) (^о2 — начальная молярно-объемная концентрация О2) от температуры за ударной волной. Видно, что кинетическая схема правильно предсказывает величину tign в диапазоне температур

ОА 1 1

-а = -мвх; ртг К г - ^) к/г

°ь р г=1

/а \ У/Зг

О 1 Д 2 - 3

800

1300

1800

2300

2800 Т

Рис. 1. Расчетное время задержки воспламенения (3) в зависимости от температуры за ударной волной в сравнении с экспериментом [5] (1), [6] (2)

1000 < Т < 2500 К и занижает при Т < 1000 К. Очевидно, что в диапазоне температур Т < 1000 К для лучшего соответствия эксперименту [6] необходима корректировка констант скорости.

Численный метод. При решении системы одномерных нестационарных уравнений неравновесной газовой динамики (1)-(3) для аппроксимации по времени применен конечно-разностный метод типа универсального алгоритма [7]. Для пространственной аппроксимации использована ТУБ-схема третьего порядка точности с расщеплением вектора потоков по Ван Лиру [8]. Детали алгоритма подробно изложены в [9].

После определения динамических и термодинамических параметров смеси на (п + 1)-м шаге по времени решалась задача определения относительных массовых концентраций компонентов смеси. Для этого на каждом шаге по времени решалась задача Коши для обыкновенных дифференциальных уравнений химической кинетики, записанных вдоль траекторий газовых частиц. Реализация осуществлялась решателем жестких систем И,АВАи5. Для интерполяции значений относительной массовой концентрации в узлы сетки использовалась кубическая сплайн-интерполяция. Оказалось, что использование кубического сплайна для интерполяции разрывных решений, возникающих на контактных разрывах, разделяющих различные газы, приводит к появлению осцилляций. Эти осцилляции удалось подавить с помощью сглаживающего фильтра. На рис. 2 приведен пример интерполяции разрывного решения без применения (кривая 1) и с применением (кривая 2) сглаживающего фильтра.

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

ди „. . ди

дн + ' (и) дХ = а (4)

Введем подвижную сетку (х0, Х\, ... , хп}, в которой узлы уплотняются или разрежаются в зависимости от градиента и. Назовем х физической, а £ расчетной координатой. Узлы сетки в плоскости £ располагаются равномерно. При этом х = х(£,£), £ £ [0,1],

0.8

0.4

о 1 - 2

0 ООССССССССООООООо

".........о '

0.4 0.5

Рис. 2. Пример интерполяции разрывного решения

х(0,Ь) = 0, х(1,Ь) = 1. Перепишем уравнение (4) в расчетной координате £ следующим образом:

ди (дх Л ди д£

дь "Iдь - ;(и)) д£дХ = 0'

дх д дх

где — = — выражает изменение сетки; — монитор, характеризующий сте-

дь дМ д£

\

'ди^ 2

1 + И ) ; а — коэффи-

пень сгущения сетки в зависимости от градиента и, циент, с помощью которого можно регулировать степень сжатия сетки

Тестирование алгоритма адаптации сетки

Тест 1. Нелинейное уравнение переноса. Рассмотрим применение алгоритма адаптации к решению данного уравнения. В качестве начальных данных ( в отличие от [3]) примем разрывное распределение в виде

, пч / 1, х < 0.5,

и (х, 0) = ^ ' ' [0, х > 0.5.

На рис. 3 приведено сравнение результатов, полученных без сгущения (кривая 1) и со сгущением (кривая 2) сетки. Видно, что более качественным получается решение с использованием алгоритма адаптации. Отметим, что в обоих расчетах бралось 200 узлов в пространственном направлении. Это значит, что шаг равномерной сетки равнялся 5 • 10-3, а в расчете с неравномерной сеткой параметр а подбирался так, чтобы минимальный шаг равнялся 5 • 10-4. Таким образом, чтобы провести расчет с шагом 5 • 10 без использования алгоритма адаптации, потребовалось бы в 10 раз больше узлов.

-4

Рис. 3. Сравнение расчетов, полученных без сгущения (1) и со сгущением (2) сетки

а ЛЖт1п • 10-4 Лхтах • 10 ' N ЛЬ • 10-3

0.04 5 5.1 3 2.2

0.06 4 5.2 3 1.8

0.08 3.6 5.3 2 1.2

0.12 3.0 5.3 2 1.0

0.15 2.65 5.3 2 0.8

,-4

1

0.02 0.06 0.10 0.14 0.18 0.22 0.26 0.30 °

Рис. 4. Минимальный шаг сетки в зависимости от параметра а

В таблице приведены данные параметрического анализа алгоритма адаптации в зависимости от параметра а. Здесь Ахт;п, Ахтах — минимальный и максимальный шаг соответственно; N — число узлов, укладывающихся на разрыв; АЬ — ширина разрыва.

Из рис. 4 видно, что при увеличении параметра а максимальный шаг расчетов меняется слабо. В то же время минимальный шаг с увеличением а уменьшается более заметно, имея, вообще говоря, тенденцию к стабилизации. Кроме того, рост параметра а приводит и к уменьшению ширины волны примерно в 3 раза при заданном количестве узлов. Отметим, что значения а зависят от количества используемых в сетке узлов.

Тест 2. Задача о распаде произвольного разрыва. Покажем преимущества применения механизма адаптации на примере задачи о распаде произвольного разрыва в замороженной смеси (без учета химических реакций). Начальные данные в этой задаче выбирались в соответствии с известным тестом Сода:

При решении уравнений (1) в качестве величины, в зависимости от градиента которой сгущалась сетка, выбирали плотность. На рис. 5 показано распределение плотности, полученное на равномерной (а) и неравномерной (б) сетках со 100 узлами в пространственном направлении. Видно, что и в случае математической модели газовой динамики решение, полученное с использованием алгоритма адаптации, позволяет существенно лучше разрешить волновую картину течения. Фронты контактного разрыва и ударной волны получаются более крутыми. Кроме того, читатель может увидеть, как узлы сетки растягиваются в областях с малым или отсутствующим градиентом и сгущаются при больших градиентах решения. В этом расчете параметр а, ответственный за степень сгущения сетки, выбирался равным 0.05. При таком значении параметра а максимальный шаг сетки был в 10 раз больше минимального, а минимальный шаг адаптивной сетки (Ахт;п = 0.008) был примерно в 5 раз меньше, чем в расчете с равномерной сеткой (Ах = 0.04). Соответственно расчет с шагом 0.008 без сгущения сетки занял бы в 5 раз больше времени, чем со сгущением.

х < 0,

х > 0,

и (х, 0) = 0, х € [-2, 2].

PtrÖbÜJÜJJjUjUjUÜJÜJl^i

о

0.8

о-о о о о о о<

0.6

0.4

0.2

Охрхсагазхсшшгь

-2

0.8

0.6

0.4

0.2 -2

юоооооооо

-1

X

Рис. 5. Распределение плотности, полученное на равномерной (Дж = 0.04) (а) и адаптивной (Джтш = 0.008) (б) сетках

а

Рис. 6. Распределение узлов в расчетной области

На рис. 6 показано распределение узлов в расчетной области при 100 узлах сетки и параметре а = 0.05. Меньший угол наклона графика свидетельствует о большем количестве узлов в этой области. Видно, как сетка сгущается к особенностям течения: ударной волне (участок А), контактному разрыву (Б) и волне разрежения (О).

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

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

Итак, рассмотрим пространство, заполненное смесью кислорода и водорода. По смеси распространяется плоская детонационная волна за замороженным фронтом ударной волны, в которой водород воспламеняется и начинает гореть (приближение Зельдовича, Неймана, Деринга — ЗНД). Уравнения, описывающие течение смеси газов в системе координат, сопутствующей ДВ, имеют вид трех законов сохранения: массы, импульса и энергии, дополненных кинетическими уравнениями химических процессов (3) и соответствующим УРС:

pu = pouo = C1, p + Ciu = Po + CiUo = C2, (5)

e + p/p + u2! 2 = eo + po/po + u2oj2 = C3.

Посредством выражения искомых функций через массовую скорость и относительные массовые концентрации компонентов известным образом система уравнений (3), (4) может быть сведена к системе ß +1 обыкновенных дифференциальных уравнений — для описания изменений скорости смеси и кинетики процесса:

du

dz

R

cv M,

ß

E

m a=1

1

(cVaT + hoa — cPa Too) - rTJ2tF

a=1 M«,

u2 — cf

Wa

(6)

ddta = p Ma E pmr dt p Г=1

(var Var )

ß ff \Vßr ß ^Л Мв) — kbrpmr-mrä

Мв

ßr

Wa

Здесь Мт и Cf — молекулярный вес и замороженная скорость звука смеси, Я — универсальная газовая постоянная. Выражения остальных газодинамических параметров через и, £а в зоне химической релаксации находятся из законов сохранения массы, импульса и энергии (4).

Сформулируем краевую задачу для системы (5). Найти разрывную вектор-функцию Ф(и, р,р, е, ), которая удовлетворяет системе (5) в области (-то, и число Маха детонационной волны — собственное число задачи М0, такие, что

и,£ ^ 0, и ^ Щк, ^ £ак < 1 при X ^ + ТО,

и = ио, £а = £а0 при X ^ -ТО,

в виде решения Чепмена—Жуге (СЛ) (Mfk < 1, Мек = 1) или решения типа пересжатой детонации (Мек < 1) или недосжатой детонации (Мек > 1). Здесь Мек > 1 — равновесное число Маха. Индекс f указывает на замороженную, а индекс е — на равновесную величины. Индексом к обозначается конечное состояние, а индексом 0 — состояние перед фронтом ДВ.

При х = +0, т. е. на замороженной ударной волне, состояние определяется из условий сохранения массы, импульса и энергии: и = Uf, р = pf, р = pf, е = ef.

Стационарное решение для детонации Чепмена—Жуге в стехиометрической смеси водорода с кислородом, разбавленной на 50 % аргоном, было получено на основе

численного решения данной краевой задачи на собственные значения. В детонационной трубе перед ДВ задавались, кроме того, следующие термодинамические параметры смеси газов: Т0 = 297 К, р0 = 0.2 атм. Не останавливаясь на физической интерпретации полученных данных, рассмотрим лишь методические аспекты исследования. Решение, достигнутое численным методом пакетом И,АВАи5 в стационарной постановке, принималось в качестве начальных данных (задача Коши) для нестационарной модели. О правильности построения численного алгоритма можно судить по сохранению стационарности распространения ДВ во времени. На рис. 7 показано начальное положение ДВ, полученное при решении стационарной задачи (кривая 1), которое служит данными Коши. Здесь же приведены кривые 2 и 3, позволяющие сопоставить расчетные концентрации радикала ОН, полученные в рамках нестационарного подхода при описании распространения ДВ без использования (кривая 2) и с использованием (кривая 3) алгоритма адаптации сетки. Как видно, в расчетах без сгущения и со сгущением сетки в направлении пространственной координаты х взято одинаковое количество узлов: 2.5 • 103. Параметр а в расчете со сгущением был подобран так, чтобы минимальный

Рис. 7. Сравнение расчетных концентраций радикала ОН а б

Рис. 8. Распределение плотности (а) и температуры (б)

шаг расчетной сетки Axmin был равен примерно 1 • 10-5. Чтобы добиться такого шага в расчете без адаптации сетки, потребовалось бы взять в направлении пространственной координаты x примерно 4-104 узлов, т. е. почти в 20 раз больше. Из рис. 7 видно, что расчет без адаптации сетки неправильно воспроизводит профиль стационарного решения. Это очевидно связано с тем, что шаг расчетной сетки в данном случае недостаточно мал, чтобы разрешить зону химической реакции.

Далее, на рис. 8, показано решение задачи Коши для уравнений неравновесной газовой динамики (3) с такими параметрами, как плотность и температура в различные моменты времени. Видно, что:

— стационарность решения для данных параметров практически не нарушается;

— скорость движения детонационной волны Чепмена—Жуге в расчете равна примерно 2150 м/с, что примерно на 2 % отличается от скорости, известной из эксперимента [10], и не отличается от численных данных стационарного расчета.

Структура решения показывает типичные особенности ДВ в рамках подхода ЗНД: фронт замороженной ударной волны, зона протекания неравновесных химических реакций. Отметим, что в литературе описывается продольная неустойчивость детонационной волны, которая получается в модельных одностадийных и двухстадийных кинетических законах окисления реагирующего компонента [11, 12]. Данный тип неустойчивости интерпретируется как отражение ячеистой структуры детонации в одномерном нестационарном подходе. Он возникает при некоторых условиях на значение энергии активации и теплоподвода, при решении задачи об ослаблении пересжатой детонационной волны волной разрежения. Оказалось, что возможен, в зависимости от значения энергии активации, выход как на режим Чепмена—Жуге, так и на колебательный аттрактор. Отметим, что в наших расчетах колебательный характер профилей решения не наблюдался.

Список литературы

[1] TlEN J.H., Stalker R.J. Release of chemical energy by combustion in a supersonic mixing layer of hydrogen and air // Comb. Flame. 2002. N 130. P. 329-348.

[2] Шокин Ю.В., ЛисЕйкин Л.Д. и др. Методы римановой геометрии в задачах построения разностных сеток. Новосибирск: Наука, 2005. 254 с.

[3] Ceniceros H.D., Hou T.Y. An efficient dynamically adaptive mesh for potentially singular solutions // J. of Comput. Phys. 2001. N 172. P. 609-639.

[4] Бедарев И.А., Федоров А.В. Сравнительный анализ трех математических моделей воспламенения водорода // Физика горения и взрыва. 2006. Т. 42, № 1. С. 26-33.

[5] Schott G.L., Kinsey J.L. Kinetic studies of hydroxyl radicals in shock waves II: Induction times in the hydrogen-oxygen reaction // J. Chem. Phys. 1958. N 29. P. 1177-1188.

[6] Skinner G.B., Ringrose G.H. Ignition delays of a hydrogen-oxygen-argon mixture at relatively low temperatures // J. Chem. Phys. 1965. N 42. P. 2190-2204.

[7] Ковеня В.Ы., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, СО, 1981. 384 с.

[8] Van Leer B. Flux-vector splitting for the Euler equations // Lecture Notes in Physics. 1982. Vol. 170. P. 507-512.

[9] Бедарев И.А., Федорова Н.Н. Исследование факторов, влияющих на качество предсказания турбулентных отрывных течений // Вычисл. технологии. 1999. Т. 4, № 1. C. 14-33.

[10] Akbar R. Mach Reflection of Gaseous Detonations. PhD thesis, Rensselaer Polytechnic Institute. Troy. N. Y., 1997.

[11] Медведев С.А. Об ослаблении пересжатых детонационных волн с конечной скоростью реакции // Изв. АН СССР. Механика жидкости и газа. 1969. № 3. C. 23-30.

[12] Левин В.А., Марков В.В. Возникновение детонации при концентрированном подводе энергии // Физика горения и взрыва. 1975. Т. 11, № 4. С. 623-633.

Поступила в редакцию 24 октября 2008 г., в переработанном виде —16 января 2009 г.