CC BY
57
Вычислительные технологии
Том 12, № 2, 2007
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ФИЛЬТРАЦИОННОГО ГОРЕНИЯ ГАЗА
НА ОСНОВЕ ДВУХУРОВНЕВЫХ ПОЛУНЕЯВНЫХ РАЗНОСТНЫХ СХЕМ*
Ю.М. ЛАЕВСКИй, Л. В. ЯУШЕВА Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: laev@labchem.sscc.ru
In this paper we address a problem of numerical modeling on the propagation of a combustion front into an inert porous medium during the filtration of a combustible gaseous mixture. Possibilities of application of some new modifications of the explicit difference schemes to solution of strongly multi scale spatial-temporal problems along with the approaches to the space-time mesh adaptivity, etc. are discussed within the presented investigation. Presented numerical results demonstrate a robustness of this technique for the analysis of different physical phenomena in the processes of filtrating gas combustion.
Введение
Под фильтрационным горением газов (ФГГ) будем понимать процессы распространения зоны газофазной экзотермической реакции в химически инертной пористой среде при фильтрационном подводе газообразных реагентов к зоне химического превращения. Здесь мы приводим фрагменты работы [1, 2], в которых дан обзор по этому направлению исследований с достаточно подробным списком литературы. Дальнейшее изучение процессов ФГГ велось по различным направлениям: устойчивость стационарных режимов ФГГ [3], стационарные режимы в сферических волнах ФГГ [4], процессы с химически активной пористой средой [5] и т. д. Отметим также исследования белорусских ученых в области ФГГ (см. [6] и приведенную там библиографию).
Процессы ФГГ чрезвычайно трудны для численного моделирования. Во-первых, из-за традиционных трудностей, возникающих при численном решении задач газового горения и связанных с большой разномасштабностью протекающих при этом процессов. На эту тему имеется обширная литература (см., например, [7] и приведенную там библиографию). При газовом горении обычно выделяют три масштаба — внешний (газодинамический), средний (диффузионный) и внутренний (кинетический). Наиболее сложно для численного моделирования взаимодействие диффузионного и кинетического масштабов. Именно
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 07-01-00164а и № 06-03-32524а).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
это взаимодействие порождает такую важную макрохарактеристику, как нормальная скорость распространения пламени. Попытки численно определить эту величину не только качественно, но и количественно сталкиваются со значительными проблемами. Обсуждение этих проблем приводится во введении работы [8]. Во-вторых, возникают дополнительные трудности, связанные со спецификой собственно задач ФГГ, поскольку в этих задачах необходимо учитывать еще масштабы, задаваемые теплопереносом в пористой среде и тепловым взаимодействием газовой фазы и пористой среды. В отличие от газодинамического, эти масштабы самым существенным образом влияют на процесс ФГГ, обусловливая существование различных режимов распространения зоны горения (подробная классификация возможных режимов ФГГ приведена, например, в [2]).
В настоящей работе промоделирован режим низких скоростей (РНС), когда в основе процесса лежит эффективный перенос тепла по пористому телу из зоны тепловой релаксации (разогретые продукты газофазной реакции передают тепло в пористое тело) в зону подогрева (прогретое пористое тело нагревает горючую смесь). Таким образом, механизм распространения зоны горения обусловлен не только (и не столько) передачей тепла из зоны экзотермической реакции в зону подогрева, но и, главным образом, потоком тепла из пористого тела. В связи с этим в дальнейшем будем выделять зону подогрева, зону реакции и зону тепловой релаксации. Отметим, что в зоне тепловой релаксации горючая смесь отсутствует, но, в отличие от моделей ламинарного газового пламени, она играет чрезвычайно важную роль. Первые одномерные численные модели ФГГ были реализованы в работах [9, 10], ряд многомерных вычислительных моделей процессов, основанных на ФГГ, представлен в [11].
1. Постановка задачи
Приведем соображения (частично изложенные в [8]), определяющие цель данного исследования. Принято считать, что наиболее продуктивный инструмент численного моделирования в теории горения — балансные монотонные неявные разностные схемы с использованием сгущающихся в зоне горения сеток [12, 13]. Однако имеющая место локальная неустойчивость типа теплового взрыва делает использование неявных схем, на наш взгляд, нецелесообразным. Речь идет о локальной неустойчивости по времени в фиксированной пространственной точке решения уравнения для относительной концентрации реагирующего компонента горючей смеси п:
ж = <"-г)- (1)
Здесь и далее Ш<п,Т) — скорость химической реакции первого порядка, протекающей по закону Аррениуса,
Ш<п,Т) = ко пе-Е/КТ, (2)
где Т — температура газа; к0 — предэкспонент; Е — энергия активации; Я — универсальная газовая постоянная. В случае подобия полей концентрации и температуры (что вполне оправданно для газовых горючих смесей), т. е. при наличии равенства
О
Ье = — = 1, (3)
а
где а и О — коэффициенты температуропроводности и диффузии газа, имеет место связь между п и Т:
Т + % = То + ^ = Ть. (4)
Ср Ср
Здесь Т0 — температура "холодной" смеси; ср — удельная теплоемкость газа при постоянном давлении; Q — тепловой эффект реакции. Тогда Ш(п,Т) = и производная этой функции в точке п* имеет вид
Е
со е I 1 - —— (Ть - Т*)
I(п*) = ко 1 - ^ (Ть - Т*0
где согласно (4) Т* = Ть - (Ть - Т0)п*. При этом ^(п*) > 0 (область устойчивости), если
ап
{Е\2 Е Е
Т* ы + кТь - 2я-
Для типичных значений физических параметров горючей смеси Ть = 1500 К, Е/К = 2 • 104 К (что соответствует энергии активации Е ^ 40 ккал/моль) устойчивость имеет место только при Т* > 1402 К, что соответствует выгоранию горючей смеси до концентрации п* < 0.082. Таким образом, практически вся зона подогрева находится в области локальной неустойчивости. "Не помогает" и диффузионный перенос массы (более подробно см. в [8]). Из проведенного анализа следует важный вывод: применение неявной разностной схемы в .зоне подогрева нецелесообразно, поскольку неявные схемы абсолютно устойчивы при работе только с устойчивыми решениями. В противном случае они теряют свои преимущества по сравнению с явными схемами, оставаясь при этом значительно более трудоемкими с точки зрения реализации. С другой стороны, в зоне реакции необходимо существенное сгущение пространственной сетки для адекватного воспроизведения функции скорости реакции. При этом в режиме распространения фронта горения поведение решения аналогично поведению решений волновых задач, когда пространственная и временная переменные несут "пропорциональную нагрузку" (решение сохраняется вдоль характеристики). Это приводит к необходимости измельчения шага по времени, и лимитирующим фактором становится не устойчивость схемы, а точность вычислений. Но это означает нецелесообразность использования неявных схем и в зоне реакции. Совершенно иная ситуация возникает в зоне тепловой релаксации — здесь лимитирующим фактором при применении явной схемы является устойчивость, а следовательно, целесообразно в разностное уравнение для концентрации введение какой-либо регуляризации. Таким образом, эффективность используемого алгоритма определяется компромисом между этими, вообще говоря, противоречивыми требованиями.
При использовании явных схем (в той или иной редакции) в силу условий устойчивости шаг по времени ориентирован на самый мелкий пространственный шаг. В отличие от локального характера аппроксимации, устойчивость — это глобальное свойство сеточной системы (характеризуется спектром линеаризованного сеточного оператора). И если в зоне реакции такой шаг по времени оправдан (см. выше), то в зоне подогрева, где пространственный шаг значительно больше, устойчивость становится главным лимитирующим фактором, делающим явные алгоритмы крайне неэффективными. Указанное противоречие может быть устранено использованием многоуровневых явных схем, предложенных и исследованных в работах [14-16]. В частности, двухуровневые схемы организованы
таким образом, что расчет в зонах подогрева и тепловой релаксации ведется с шагом по времени а в зоне реакции — с шагом т ^ Главное свойство этих схем состоит в локализации условий устойчивости по группам переменных: если самосопряженный положительно полуопределенный пространственный оператор сеточной задачи представлен в блочном виде
А =( АО ' ЦАнИш <ЦЛ22Ц(2), (1.5) (5)
то устойчивость обеспечивается неравенствами
Д*||Ап||(1) < 1, Т||А22|(2) < 1. (6)
Отметим, что разработанная теория не касается имеющегося в нашем случае отсутствия самосопряженности и нарушения положительной полуопределенности, связанного с описанной выше возможностью возникновения локальной неустойчивости в уравнении для концентрации. Тем не менее данная методика хорошо себя зарекомендовала при расчетах динамики ламинарных газовых пламен [8]. Кроме того, в зоне реакции целесообразно использование легко обратимого регуляризатора. Такого типа модификация двухуровневых схем была рассмотрена в работе [17]. Эту модификацию мы и называем двухуровневой полунеявной разностной схемой.
Предлагаемая статья организована следующим образом. В разд. 2 формулируется задача, описываемая одномерной моделью ФГГ. В качестве исходной переменной принята температура газа Т. Но в случае горения с числом Льюиса для газа Ье = 1 имеет смысл перейти к энтальпийной постановке, используя вместо переменной Т переменную Н = Т + (Т — То)п. Отметим, что величина срН является полной энтальпией газовой смеси. Именно такая постановка автоматически обеспечивает баланс сеточной системы по отношению к аппроксимации функции скорости реакции. В случае же постановки задачи в терминах "температура—концентрация" и использования неявных схем "стандартная" аппроксимация функции скорости реакции (явная в уравнении для температуры газа: Ш(Тп,пп), и неявная по п в уравнении для концентрации: —Ш(Тп,пп+1)) приводит к нарушению энергетического баланса. В разд. 3 приводится двухуровневая полунеявная разностная схема по времени, лежащая в основе предлагаемой методики. Здесь же рассматриваются пространственная конечно-разностная аппроксимация и алгоритм адаптации сетки. В основе такой адаптации лежит воспроизводимость функции скорости реакции. Подчеркнем, что не большие градиенты температуры и концентрации требуют сгущения сетки, а большие градиенты потоков тепла и массы, величина которых и обусловлена скоростью химичесой реакции. В разд. 4 приводятся численные эксперименты, демонстрирующие как вычислительные характеристики алгоритма, так и физические результаты, полученные при варьировании параметров системы.
2. Одномерные уравнения ФГГ
Рассматрим процесс распространения зоны газофазной экзотермической реакции в химически инертной пористой среде при фильтрационном подводе газообразных реагентов к зоне химического превращения в одномерной постановке. При моделировании данного процесса гомогенизация пористой среды осуществляется таким образом, что в каждой точке отрезка длины Ь определены температуры пористой среды Т3 и газа Тд (см., например, [1]). Простейшая двухтемпературная модель ФГГ включает уравнение переноса
С'р-ж = А- их? + 1-—(Т» -Т-); (1)
и± д2Т дТ —
С»= А» -дф - СЯр»V + - (Т - т) + Яр»и-'(п. Т); (2)
£ = ^ - ^ - и(пт). (3)
тепла в пористой среде и уравнения переноса тепла и массы реагирующего компонента газовой смеси:
дТ8 Л д2Т3 а
—1 = А.,-- +--
дЬ дх2 1 - —
дТ д2Т дТ а
дТ» = а д Т - С р v —» + — дЬ » дх2 » » дх —
д'П я д2П дп дЬ дх2 дх
Здесь рг, Сг, Аг — плотность, удельная теплоемкость при постоянном давлении и коэффициент теплопроводности 1-й фазы (г = в, д); — — пористость (отношение всего объема к суммарному объему порового пространства); а — интенсивность межфазного теплообмена; V — скорость подачи горючей смеси (V > 0). Все физические параметры полагаем постоянными, т. е. в приведенной модели пренебрегается тепловым расширением газа и фильтрационными эффектами, вызывающими изменение давления. Укажем краевые условия для приведенной системы: — при х = 0
Т- = Т» = То, п =1; (4)
при х = Ь
дТ _
А» ^5 = 0, ^ = ° (5)
7Г = °, А?» = °, ^
дх дх дх
В дальнейшем предполагаем подобие полей концентрации и температуры в газе: А = СдРдЯ (Ье» = 1). С использованием этого равенства задача приводится к энтальпийной постановке. Вместо температуры газа Т» введем функцию
Н = Т» + (6)
С
Тогда система уравнений (1)-(3) принимает вид
дТ- д2Т- + ( Я
+ — Н - Т-- - - п ; (7)
дь дх2 \ с»
дН д2Н дН ( тт Я
дь - + —{Т-- Н + С» У ; (8)
дп д 2п дп т /0ч
дь = ^ - 'дх - ш(п>Н)- (9)
Здесь аг = Аг/сгрг — коэффициент температуропроводности г-й фазы, причем а» = Я; а8 = —/(1 - — )с-р-, а» = —/—^р». В соответствии с (6) в дальнейшем функция скорости реакции определяется как
и(п,Н) = и^п.Н - (10)
Краевые условия (4) и (5) принимают следующий вид: — при х = 0
Т- = То, Н = То + п = 1; (11)
С
— при х = Ь
дТ дН дп
—^ = 0, ай — = 0, = 0. (12)
дх дх дх
Главное преимущество системы (7)-(9) перед (1)-(3) состоит в том, что функция скорости реакции входит только в одно из уравнений и нет необходимости специально обеспечивать баланс между сеточными уравнениями для функций Н и п-
3. Двухуровневая полунеявная разностная схема
Основу приводимого ниже алгоритма составляет двухуровневая полунеявная разностная схема [17]- Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений:
^ + Аи = /(и), г> 0; (1)
(И
и = и0, г = 0.
(2)
Компоненты вектора и включают значения сеточных функций Т8, Н и п- Разобъем вектор и на два подвектора, содержащих две группы переменных: к первой группе отнесем значения, соответствующие сеточным узлам из зон подогрева и тепловой релаксации, а ко второй группе — из зоны реакции (с некоторым "запасом"). Тогда представлению вектор-столбца и = (м^, )т отвечают блочные представления оператора А и вектора правой части:
Ац А12
А21 А22
А
(3.3)
(3)
/(и) = (/Т(и),/2(и))
т
Здесь Т означает транспонирование. Отметим, что поскольку в правую часть не входят аппроксимации дифференциальных операторов, имеет место (/)г(и) = (/к)г((ик)г), к = 1, 2-Здесь (ик)г и (/к)г — компоненты соответствующих сеточных функций, отвечающих г-му сеточному узлу. Последнее означает, что матрица Якоби д//ди диагональная. Неравенство ||АИ||(1) ^ ||А22Н(2) из (5) соответствует сильному сгущению сетки в зоне реакции. Для введенных обозначений двухуровневая полунеявная разностная схема решения задачи (1), (2) имеет вид
п+1
дг
+ Ап< + А12иП = /1«);
иГР = < + к«+1 - О, к = 1,
Р
Е + тП
п+ ; \ и
— + 1
- и2
П+ Р
П+ Р
—
+ А^и'ГР + А22и2 Р = /2(и2"'Р), к = 0,...,Р - 1, (6)
,р - 1;
Р '
(4)
(5)
т
где Дг = рт (р — некоторое натуральное число); Е2 - единичная матрица соответствующей размерности;
п+ Р д/2 / п+ Р
П Р
ди2
и
(7)
По сути дела, речь идет о методе декомпозиции области в следующей редакции. Сначала в зонах подогрева и тепловой релаксации делается большой шаг Дг по явной схеме
и
2
2
Эйлера (равенство (4)). Затем на границе между этими зонами и расширенной зоной реакции решение линейно интерполируется на р вспомогательных слоев (равенство (5)). Далее, в соответствии с равенством (6) в зоне реакции по явной схеме Эйлера с регуляризацией делаются р шагов величины т = Д*/р. Последний шаг эквивалентен использованию неявной по функции скорости реакции схемы Эйлера с одной итерацией по методу Ньютона.
п+ к
Отметим, что матрица р диагональна, а следовательно, производительность схемы (6) не ухудшается по сравнению с явной схемой Эйлера. В этом случае отсутствует регуляризация в зоне тепловой релаксации. Ниже этот недостаток будет устранен.
Рассмотрим пространственную аппроксимацию уравнений (7)-(9). Для этой цели на отрезке [0, Ь] введем базовую равномерную сетку иь = (хг = гк}^ с шагом к = Ь/З. Далее, на интервале [xj, ж] С [0,Ь] введем, вообще говоря, неравномерную сетку иП0 = (¿^¿П-1. При этом 0 < ] < / < 3, ] = ] (п), / = /(п). Рассмотрим два дополнительных узла базовой сетки ¿о = ж^п)-1 и п = х^+ь Обозначим иП = иП,0 и (ж^п)-1, жг(п)+1} = {¿г}/^. Пусть кг = _ гг-1, г = 1,..., 3П, — шаги сетки и^. Интервал [х^)^, Жг(П)+1] будем называть расширенной зоной реакции. В соответствии с методом (4)-(6) для каждой из зон выпишем разностные соотношения, аппроксимирующие уравнения
(7)-(9).
В зонах подогрева и тепловой релаксации значения ТП+1 и ЯП+1 определяются равен-
ствами
Т™+1_т
П
а,г
Д*
к2 т,-, _ 2Т™ + ТаГг+О + аа ( Я? _ _ Я П
ЯГ1 _ Яп
Д
к2 (ЯГ-1 _ 2ЯГ + Я+1) _ к (ЯГ _ ЯТ-0 + а (т™ _ ЯП + Я пП
к
(8) (9)
1,...,^" (п) _ 1, /(п) + 1,...,3 _ 1,
и в соответствии с краевыми условиями (11), (12)
гр'П+1 _ гр
Т а,0 = Т0
Я
п+1
Я
Т П+1 Т
п
а,П а а
2Д*
к2 стт^ -1 _ тп + О1 _ тт^ _ пп),
Я
ЯП+1 _ ЯП а0 + ^к
(ЯП-1 _ ЯП) + ^ (ТПп _ ЯП + ^ пП
2Д*
к2
Я
В зоне подогрева концентрация на (п + 1)-м слое определяется равенством
пГ+1 - = к§ (пП-1 _ 2пП + пПи) _ к (пП _ пП-х) _ ^(пП, ЯП),
Д*
(10)
а в зоне тепловой релаксации
г = 1,---,Лп) _ 1,
(1 + Д*^) пГ+Д пГ = к- (пП-1 _ 2пП + пП+1) _ к (пП _ пП-1) _ ^(пП, ЯП), (11)
Д*
г = /(п) +1,..., з _ 1,
а
р
а
где
З'ц
(С,Н™) = ко е
—е/д(я?—ср п?)
1 - ЕО - -2 %
(12)
При этом краевые условия имеют вид
1,
„П+1
(1 + ДШп) ^
2Д£
Как уже говорилось, в зоне реакции и тем более в зоне тепловой релаксации (т/П мало) величины Д™ положительны, а следовательно, могут использоваться для регуляризации. Здесь мы модифицировали схему (4)-(6), осуществив регуляризацию явной схемы для концентрации в зоне тепловой релаксации.
Рассмотрим разностные соотношения в расширенной зоне реакции. Для значений сеточных функций , заданных на сетке введем обозначение: иг = г = 0, При этом Ио = и,(га)-1 и = иг(п)+1. Тогда
Г р
7
- Н
р
Т
а£
я*
т р - т
р
Лг
- .. Г р - г
8,г-1
г+1
р
Лг
+
+ «8
Н р
- н\
Лг
чп+— О —
Н" р - Т р - _ н Р
г /г
-га+ -р
п+ -
р
-га+- .
Нг+1Р - Н
п+ -р
Н р - Н
п+ -р
г- 1
Т
- ^ Н..+р - Н
Яг
1 + тД- р
п+ -р
г—1
+ а
о р
Яг
Лг+1 Лг
р - Нг Р , О и™+ + С9 'г ,
- Нг+1Р - Н Р Нг Р -
Л г+1 Лг
Н р '/г—1
(13)
(14)
V Лг
Т Р
и р
иг— 1
1,...,7га - 1,
- Ш(Иг к=0
п+ -р
-га+ -ч т р
НГ р), р - 1,
(15)
где = (Лг + Лг+1)/2. Соответствующая равенствам (3.5) интерполяция осуществляется в сеточных узлах ж^(га)—1 и Жг(га)+1. Напомним, что т = Д/р.
Рассмотрим вопрос об адаптации разностной сетки. В основе адаптации лежит воспроизводимость функции скорости реакции Ш(п, Н). Как функция пространственной сеточной переменной Ш(п, Н) имеет характерную дельтообразную форму с очень малым носителем. Воспроизведение этой формы (выбор сетки на интервале [ж^), Жг(га)]) может осуществляться различными способами (см., например, [8]). Не останавливаясь на этом вопросе и решив его экспериментально в следующем разделе, опишем алгоритм выбора интервала [ж^(га), ж^)].
В начальный момент времени на базовой сетке ^ по известным зачениям и Н0 найдем тах Ш(г^Н0), реализующийся в сеточном узле жг0. По заданному фиксированному
0<г< J
р
р
9
а
9
р
а
9
значению ^ = 0(Л) рассмотрим интервал [-¿О;-^'] С [0,Ь], ¿0 = — ^/2, ¿0' = + ^/2, на котором введем сетку ¡¿о с постоянным шагом Ло = ^/К, где К — заданное натуральное число. Затем рассмотрим минимальный интервал [ж^-0 , ж^0] С [0,Ь] такой, что [Хо,] Э [¿0,-0]. Пусть у(0) = уо — 1, 1(0) = /о + 1. На интервалах [ж^о),—] и [—'^(о)] введем сетки ¡0 и ¡0' с постоянными шагами Л0 = (¿0 — Х/(0))/К' и Л-0' = (жг(0) — ¿0')/К'' соответственно. Здесь и далее К' и К'' — заданные натуральные числа. Пусть
¡м = ¡0 и ¡0 и ¡0' = = ¡М и {ж^(0)-1,жг(0)+1}-
На построенной сетке определим значения Т0, Н0 и г—0 как результат интерполяции заданных на базовой сетке начальных сеточных функций Т°, Н0 и п°. Для определенных таким образом функций снова вычислим величину тах Ш(-т^0,-^0), реализующуюся
0<г<^
в сеточном узле € . Затем по описанному выше алгоритму повторяем построения, заменив ж0 на .
Теперь рассмотрим переход от сетки ¡П к сетке ¡П+1. Сетка ¡П характеризуется интервалом [¿П,ТП], индексами у(п), /(п) и шагами Лп, Л-П и Л-П. Аналогично предыдущему вычислим величину тах Ш(¿¡г+1, Дгп+1), реализующуюся в сеточном узле ^¿п+1 € ¡п. Пусть
тп = (¿п + ¿п)/2 — центр интервала [¿п, ¿п]. Если |г^п+1 — тп| < ^/6, то ¡п+1 = ¡п т. е. сетка остается без изменений. В противном случае определяется новый интервал [¿п+1, ¿¿^+1] длины где ¿+1 = ¿п + 81§п(^г„+1 — ¿¿п) К0Лп, Ко — заданное натуральное число. При этом ¿п+1 = %+1Эти равенства означают, что перемещение интервала, содержащего максимум функции скорости реакции, может происходить с ограниченной кусочно-постоянной скоростью. Тогда Лп+1 = ^/К = Лп = • • • = Ло. Таким образом, на интервале [¿п+1, ¿^+1] введена сетка сгп+1 с постоянным шагом Ло. Дальнейшие действия соответствуют алгоритму, описанному выше. Рассматривается минимальный интервал [ж?п+1 , Жгп+1 ] С [0,Ь]
такой, что [ж^п+1,жгп+1] Э [¿¿Л+и¿¿Л+х]- Пусть + 1) = уп+1 — 1, /(п + 1) = ¿п+1 + 1 На
интервалах [ж^(п+1), ^+1] и [¿п+1 , Хг(п+1)] введем сетки сГп+1 и ¡п+1 с постоянными шагами ¿п+1 = (¿п+1 — ж,(п+1))/К' и ¿п+1 = (жг(п+1) — ¿п+1)/К'' соответственно. Далее,
¡У = ¡п+1 и ¡п+1 и ¡п+1 = 1-1, ¡п+1 = ¡У и {жЛп+1)-1, жг(п+1)+1}.
На построенной сетке ¡п+1 интерполяцией переопределяем значения искомых сеточных функций Т^1, -Нп+1 и ¿п+1. Отметим, что /п+1 = К + К' + К'' + 2 = /п = • • • = /о.
4. Численные эксперименты
Изложенную выше методику применим для решения задачи (7)-(12). Прежде всего рассмотрим вопрос об экспериментальном выборе параметров дискретизации — параметров пространственной дискретизации Л, К, К', К'' и Ко и дискретизации по времени Д£ и т (целочисленного параметра р). Продемонстрируем расчеты, дающие представление о физике РНС. В частности, рассмотрим задачу о существовании режима низких скоростей в зависимости от соотношения скорости подачи горючей смеси V и интенсивности межфазного теплообмена а. В качестве начальных данных на интервале [0, Ь] заданы функции
Т0(ж) = То, Н0(ж) = То + Я п0(ж) = 1, ж < Ь/2,
сд
Т0(х) = Т1, Н0(х) = То + Q, п0(х) = 0, ж>Ь/2.
сд
Все расчеты проводятся при следующих значениях физических параметров: Ь = 0.1 м, Т0 = 300 К, Т1 = 1000 К, Q/cg = 1100 К,
р, = 103 кг/м3, с, = 103 Дж/(кг • К), Л, = 2 Вт/(м • К),
Рд = 1 кг/м3, Сд = 2 • 103 Дж/(кг • К), Лд = 0.1 Вт/(м • К),
т = 0.5 б/р, к0 = 1011 с-1, Е/К = 2 • 104 К.
По приведенным параметрам вычисляются коэффициенты температуропроводности а, и ад. Напомним, что Ьед = 1. Параметры V и а в дальнейшем будут варьироваться.
Выбор сеточных параметров. Здесь основная цель расчетов — подбор алгоритмических параметров, обеспечивающих устойчивость и воспроизводимость численных результатов. Используем значения скорости подачи горючей газовой смеси V = 1 м/с и интенсивности межфазного теплообмена а = 7.5-107 Дж/(м3- с- К). Приведем погрешности решения и функции скорости реакции в окрестности зоны реакции при варьировании некоторых сеточных параметров. В качестве "точного" решения, с которым осуществляется сравнение, используются результаты расчета при
к = Ь/32, 4 = к, К = 128, К' = К''К0 = 16, Д* = 10-5 с, р = 32. (1)
При этом "точное" решение обозначается как Т,ех, Нех и пех. На приведенных рисунках указаны абсолютные значения погрешности для температуры пористой среды Т,, характеризующей энтальпию газа функции Н и функции скорости реакции Ш(п, Н), полученные на момент времени = 10 с. Под погрешностью понимаются сеточные функции £у, £я, е^ со значениями
(е )п = (Та)? - (ТрП (е )п = (Н)? - (Нех)? (ет ^ (Т еж)? , (еЯ ^ (Неж)га ,
(е^)? = Ш[(п)?, (Н)?] - Ш[(пех)?, (Нех)?],
т. е. погрешности для функций Т и Н относительные, а для функции скорости реакции погрешность — абсолютная. Характерное значение максимума функции скорости реакции составляет примерно 5 • 104 с-1.
На рис. 1 показаны погрешности при варьировании количества внутренних временных слоев р в расширенной зоне реакции (величины внутреннего шага по времени т) — все параметры из (1) сохраняются, кроме р. Расчеты показывают одинаковую динамику положения фронта реакции для всех вариантов, т. е. сравнение результатов корректно в том смысле, что решения сравниваются в одних и тех же точках пространственной сетки. Аналогичная ситуация имеет место для экспериментов, проиллюстрированных на рис. 2, где показаны погрешности при варьировании шага по времени Д*; все параметры из (1) сохраняются, кроме Д* и р. В качестве "точного" использовалось решение при Д* = 2.5 • 10-6 с и р = 8. Параметр р подбирается таким образом, что величина внутреннего шага т не меняется, т. е. внешним шагам Д* = 5 • 10-6 с, Д* = 10-5 с и Д* = 2 • 10-5 с соответствуют значения р = 16, р =32 и р = 64, и для всех вариантов т = 3.125-7 с. Для "точного" решения используется это же значение т, т.е. для Д* = 5 • 10-6 с полагается р = 8. Из приведенных результатов можно сделать вывод о высокой точности "точного" решения.
Рис. 1. Погрешности ет, £и и е^ (а — в соответственно), варьирование т.
в
Моделирование РНС. Ниже приводятся расчеты, демонстрирующие физику процесса ФГГ в режиме низких скоростей [1, 2]. В качестве сеточных параметров использовались значения (1). На рис. 3 показана температура газовой фазы в моменты времени г = 3, 6 и 9с. Вычисления проводились при скорости подачи горючей смеси V = 3м/с и интенсивности межфазного теплообмена а = 3.75 • 107 Дж/(м3-с- К). Примечательная особенность моделируемого процесса — существенное превышение адиабатической температуры Ть = 1400 К во фронте горения: Тд,тах « 1600 К. Эта особенность ФГГ хорошо известна как из теории, так и из эксперимента [1, 2]. Отметим, что решение на рис. 3, а продемонстрировано только на характерном интервале длины — примерно 0.01 м. По мере распространения фронта горения температура пористой среды в зоне тепловой релаксации поднимается от начальной температуры Т1 = 1000 К до некоторой равновесной температуры Те. Именно этот процесс продемонстрирован на рис. 3, б. Как показано в [2], имеет место
гр пт . v - и Я (1 - ш)р^
Те = То +----—— •—, а —
V - (1 + а)и Од ШРд Од
В соответствии с приведенными выше значениями физических параметров а = 500. Далее, расчеты показывают, что при скорости подачи горючей смеси V = 3 м/с скорость фронта горения и ~ -1.6 • 10-3 м/с , и тогда Те « 1170 К. Примерно это значение и дает
1,0ОЕ-О5
9ДЮЕ 06 ■
8.00Е-06 \
7,ООЕ-О6 \
6.00Е-06 \
5.00Е 06 ■ \
4,ООЕ-Об \ —(И=2Е-05
З.ООЕ-Об ■ \ -о-1И=1Е-05
\ \ —1И=5Е-06
2,00Е-06 ■ Ч X
1,00Е 06 аооЕ+оо
4.05Е 02 4,06Е02
4,07Е-02 4.09Е02
а
4.10Е 02 4,12Е-02
4.06е-02 4,07е-02
в
Рис. 2. Погрешности ет, £и и е^ (а — в соответственно), варьирование ДЬ.
а
Рис. 3. Температура газа (а) и пористой среды (б): кривые 1-3 соответствуют Ь = 3, 6 и 9 с.
приведенный на рис. 3, б расчет.
На рис. 4 приведены зависимости скорости фронта и от скорости фильтрации V при различных значениях интенсивности теплообмена а. Характер кривых полностью соответствует теории процесса и экспериментальным данным, приведенным в [2]. При а =
аяиЕ+оо
-2.00Е M
-4.00Е M
-S00E-M
-В,ODE-Ю
-100E 03
-1.2DE 03
-140E03
геоЕ-оз
-180E03
—*-al=B.75e-»06
г.ооЕ оз
Ü.10 а4в азо 1.25 1.75 225 275 325 3L75 4.25 475
Рис. 4. Зависимость u = u(v) при различных а.
8.75 • 106 Дж/(м3^ с К) и v ~ 1.5 м/с происходит срыв РНС и горение переходит в режим высоких скоростей, т. е. фронт горения распространяется в режиме обычного газового горения, практически без теплового участия пористой среды. Данное явление ранее уже было описано как теоретически, так и экспериментально (см. [2]). В этом случае интенсивности межфазного теплообмена и расхода горючего не хватает для организации эффективного теплового взаимодействия газовой смеси и пористой среды.
Список литературы
[1] Бабкин В.С., ЛАЕвский Ю.М. Фильтрационное горение газов // Физика горения и взрыва. 1987. Т. 23, № 5. С. 27-44 [Combustion, Explosion and Shock Waves. 1987. Vol. 23. P. 531547].
[2] Лаевский Ю.М., Бабкин В.С. Фильтрационное горение газов // Распространение тепловых волн в гетерогенных средах / Под ред. Ю.Ш. Матроса. Новосибирск, 1988. С. 108-145.
[3] Минаев С.С., Потытняков С.И., Бабкин В.С. О неустойчивости фронта пламени при фильтрационном горении газов // Физика горения и взрыва. 1994. Т. 30, № 1. С. 49-54.
[4] Какуткина Н.А., Бабкин В.С. Характеристики стационарных сферических волн горения газа в инертных пористых средах // Физика горения и взрыва. 1998. Т. 34, № 2. С. 9-19.
[5] Laeysky Yu.M., Babkin V.S. On the theory of the travelling hybrid wave // Combust. Sci. and Technol. 2001. Vol. 164. P. 129-144.
[6] Добрего К.В., ЖдАнок С.А. Физика фильтрационного горения газов. Минск: Ин-т тепло-и массообмена НАНБ, 2002.
[7] Математическая теория горения и взрыва / Я.Б. Зельдович, Г.И. Баренблатт, В.Б. Ли-брович, Г.М. Махвеладзе. М.: Наука, 1980.
[8] ЗоткЕвич А.А., ЛАЕвский Ю.М. Численное моделирование распространения ламинарного пламени на основе двухуровневых явных разностных схем // Вычисл. технологии. 2006.
Т. 11, № 6. C. 31-43.
[9] ДрОБЫШЕВИЧ В.И. Математическое моделирование процесса формирования гибридной волны горения // Тепло- и массообмен в капиллярно-пористых телах. Минск: ИТМО, 1996. Т. 7. С. 146-150.
[10] DROBYSHEVICH V.I. Mathematical modeling of non-stationary hybrid combustion wave // Advanced Comp. and Analysis of Combustion / Eds G.D. Roy, S.M. Frolov, P. Givi. M.: ENAS Publ., 1997. P. 114-121.
[11] Рычков А.Д., ШокинА Н.Ю. Математические модели фильтрационного горения и их приложения // Вычисл. технологии. 2003. Т. 8. Спецвыпуск. Ч. 2. С. 124-144.
[12] ШКАДИНСКИЙ К.Г. О разностном счете задач зажигания и горения с учетом диффузии и гидродинамики // Физика горения и взрыва. 1969. Т. 5, № 2. С. 264-272.
[13] Мержанов А.Г., ХАЙкин Б.И., Шкадинский К.Г. Установление стационарного распространения пламени при зажигании газа накаленной поверхностью // Прикл. механика и техн. физика. 1969. № 5. С. 42-48.
[14] ЛАЕВСкиЙ Ю.М., Банушкина П.В. Составные явные схемы // Сиб. журн. вычисл. математики. 2000. Т. 3, № 2. С. 165-180.
[15] BANUSHKINA P.V., bAEVSKY Yu.M. Multi-level explicit schemes and their stability // Rus. J. Numer. Anal. Math. Model. 2001. Vol. 16, N 3. P. 215-233.
[16] Зоткевич А.А., Лаевский Ю.М. Об одном классе двухуровневых явных схем // Сибир. журн. вычисл. математики. 2002. Т. 5, № 2. C. 163-173.
[17] Banushkina P.V., bAEVSKY Yu.M., Zotkevich A.A. On local conditions of stability for multi-level difference schemes // Proc. Intern. Conf. Comp. Math. Pt 2 / Ed. G.A. Mikhailov, V. P. Il'in, Yu.M. Laevsky. Novosibirsk, 2002. P. 334-339.
Поступила в редакцию 10 ноября 2006 г.
