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

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

Том 11, № 6, 2006

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ЛАМИНАРНОГО ПЛАМЕНИ НА ОСНОВЕ ДВУХУРОВНЕВЫХ ЯВНЫХ РАЗНОСТНЫХ СХЕМ*

А. А. ЗОТКЕВИЧ, Ю. М. ЛАЕВСКИй Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: laev@labchem.sscc.ru

The problem of numerical modeling of a laminar flame front propagation is addressed. Possibility of application of new variants of explicit difference schemes to solution of problems with significant variations of time and space scales, the ways of mesh adaptivity, etc., are discussed. Presented numerical results demonstrate a robustness of this technique for the analysis of multi dimensional effects in the combustion theory.

Введение

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-00171а) и программы NWO-RFBR (грант № 047.008.007).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

сеток под структуру волны горения, достаточно популярны и в настоящее время (см., например, [5] и приведенную там библиографию). Сразу отметим одну неточность в [5] при описании проблем численного моделирования. Узкая зона горения характеризуется не большими градиентами температуры и концентрации, а большими градиентами потоков тепла и массы. Большими градиентами температуры и концентрации характеризуется диффузионная зона (зона прогрева газа перед фронтом горения), которая не вносит трудностей в численное моделирование. Это замечание играет самую существенную роль при построении алгоритмов адаптации. Значительно усложняется ситуация при рассмотрении многомерных моделей распространения ламинарного пламени. При этом только такие модели позволяют исследовать вопросы геометрии фронта пламени, его устойчивости к пространственным возмущениям и пр. К исследованиям по численным методикам решения многомерных задач о распространении пламени и их применению следует отнести [6-9] и многие другие. Во всех этих работах используется адаптация разностной сетки под структуру фронта горения (см., например, [10] и приведенную там библиографию).

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

дТ 9 2Т + (Т ) дП пд 2П и (Т ) (01)

где Т — температура газа; п — относительная концентрация реагирующего компонента горючей смеси; ер — удельная теплоемкость газа при постоянном давлении; а — коэффициент температуропроводности; О — коэффициент диффузии; Q — тепловой эффект реакции, проходящей в соответствии с законом Аррениуса

и(Т, п) = ко пе-Е/КТ; (0.2)

к0 — предэкспонент; Е — энергия активации; Я — универсальная газовая постоянная. Чтобы не уводить читателя в сторону от основной проблемы, мы преднамеренно пренебрегаем тепловым расширением газа за фронтом реакции и как следствие неоднородностью полей скорости и давления, зависимостью параметров от температуры, используем суммарную кинетику первого порядка и т. д. Более того, в дальнейшем будем предполагать подобие полей концентрации и температуры (что вполне оправданно для газовых горючих смесей):

О , Л

Ье = — = 1. (0.3)

а

В этом случае первое из уравнений (0.1) может быть переписано в виде

9? = аЦ- © = Т + ? п. (0.4)

Отметим, что величина ср© является полной энтальпией горючей смеси. Оператор правой части этого линейного уравнения имеет отрицательный спектр, а следовательно, все его решения устойчивы. В частности, решением является адиабатическая температура Ть = То + Q/cp:

© = Ть.

Здесь T0 — температура "холодной" смеси. Для второго из уравнений (0.1) рассмотрим задачу Дирихле (другие краевые условия менее устойчивы) на характерной длине L. При этом минимальное по модулю собственное число одномерного оператора Лапласа, соответствующее самой гладкой гармонике, есть Amin = — (n/L)2. Линеаризуем функцию скорости

реакции (0.2) около значений T* = Tb — Q п* и п*. Тогда оператор линеаризованного урав-

Cp

нения для концентрации имеет собственное число вида

.(T*) = — + ko e-E/RT*(RT* (Tb — T*) — 0 ' (0.5)

v

Локальная неустойчивость имеет место при ^т;п > 0. Отметим, что Т* < ТЬ и реакция полностью завершается по достижении адиабатической температуры (п = 0). Из формулы (0.5) немедленно следует, что при Т*, близких к Ть, имеет место локальная устойчивость: ^т;п(Т*) < 0. Приведем значения ^т;п(Т*) при некоторых Т* для следующих значений физических параметров: Т0 = 300 К, Ть = 1500 К (что соответствует значению Ц/ор = 1200 К), Е/К = 2 • 104 К (что соответствует энергии активации Е « 40 ккал/моль), к0 = 1014 с-1, Б = 10-4 м2/с, Ь = 0.1 м. Вычисления показывают существование двух значений — Т' « 516 К и Т" ~ 1401 К, таких, что ^т;п(Т*) > 0 в диапазоне Т' < Т* < Т'', т. е. при этих значениях Т* решения локально неустойчивые. Отметим, что в указанный диапазон входит практически вся зона подогрева смеси перед фронтом горения.

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

Разрешение указанного противоречия может быть осуществлено использованием многоуровневых явных схем, предложенных и исследованных в работах [11-13]. В частности, двухуровневые схемы организованы таким образом, что расчет в диффузионной зоне ведется с шагом по времени Д£, а в зоне горения — с шагом т ^ Д£. Главное свойство этих схем состоит в локализации условий устойчивости по группам переменных: если самосопряженный положительно полуопределенный пространственный оператор сеточной

задачи представлен в блочном виде

А

А

А12 А22

11 А12

А

(0.6)

то устойчивость обеспечивается неравенствами

Л£||Аи||(1) < 1, Т||А22||(2) < 1.

Отметим, что разработанная теория не описывает имеющееся в нашем случае нарушение положительной полуопределенности (^т;п > 0) и эффекты подвижности сетки. Тем не менее численные расчеты распространения фронта ламинарного пламени показали хорошую работоспособность двухуровневых схем [14]. Собственно, главная цель данной работы — продемонстрировать эти расчеты.

Предлагаемая статья организована следующим образом. В разд. 1 приводится простейшая двумерная модель распространения пламени в энтальпийной постановке (аналогично (0.4)). Следует отметить, что в случае газового горения с числом Льюиса Ье = 1 именно такая постановка автоматически обеспечивает баланс сеточной системы по отношению к аппроксимации функции скорости реакции. В случае же постановки задачи в терминах "температура — концентрация" и использования неявных схем "стандартная" аппроксимация функции скорости реакции (явная в уравнении для температуры Ш(Тп, пп) и неявная по п в уравнении для концентрации — Ш(Тп,пп+1)) вызывает нарушение энергетического баланса. Приводится двухуровневая явная разностная схема (типа предиктор-корректор в диффузионной зоне и схема Эйлера в зоне горения) по времени, лежащая в основе предлагаемой методики. Пространственной аппроксимации, осуществляемой по методу конечных элементов, и алгоритмам адаптации конечно-элементной сетки посвящен разд. 2. В разд. 3 приводится ряд численных результатов, демонстрирующих главным образом неодномерные эффекты искривления фронта горения, обусловленные характером зажигания, внешнего теплообмена и пр. Там же приводятся сравнительные результаты расчетов нормальной скорости пламени ип с асимптотической формулой Зельдовича для реакции первого порядка (см., например, [1, с. 219]), демонстрирующие не только качественную, но и количественную корректность этих численных результатов. И, наконец, делаются краткие выводы и рекомендации для дальнейших исследований.

1. Уравнения ламинарного пламени и двухуровневая явная схема

Рассматривается задача моделирования движения фронта ламинарного пламени между двумя пластинами, бесконечными в направлении одной из осей. Область решения представляет собой произвольное прямоугольное сечение О = (0,Ь) х (0,/), / ^ Ь. Простейшая математическая модель данного процесса описывается системой уравнений

(1.1)

(1.2)

где T = в — Q a — коэффициент температуропроводности (напомним, что рассматри-

вается газовая смесь при Ье = 1); функция скорости реакции первого порядка Ш (Т, п) задана равенством (0.2). Физический смысл остальных обозначений указан во введении. Формально первое уравнение этой системы не зависит от второго. Фактическая зависимость в от п реализуется через краевые условия: — условия при х = 0

© = ТО + Q, п = 1;

(1.3)

— условия при x = L

условия при y

де

dy

ao

д п З'ц о

dx dx '

П Q т\ dv

п — С"v — То) = дУ cp / ду

(1.4)

(1.5)

условия при y

д П

+ ai

дУ

П Q Т\ дп

п — С"n — Tv = дУ

0.

(1.6)

Здесь а0 и щ — положительные коэффициенты теплообмена с внешней средой (с температурой Т0) через "боковые" стенки (у = 0, у = /). Кроме того, заданы начальные распределения температуры и концентрации в момент времени £ = 0:

Т = T0(x, y), n = П0 (x, y), П = T0(x, y) + Qn0(x,y)

(1.7)

Напомним, что мы преднамеренно рассматриваем простейшую модельную задачу, не учитывающую неоднородность поля скоростей.

Пространственной конечно-элементной аппроксимации приведенной задачи и способу задания триангулированной адаптивной сетки посвящен следующий раздел. Здесь лишь отметим, что краевые условия (1.3) удовлетворяются непосредственно для сеточных функций (главные краевые условия), а (1.4)—( 1.6) входят в оператор сеточной задачи и в правую часть (естественные краевые условия). Аппроксимированная по пространственным переменным задача (1.1)—(1.7) может быть записана как задача Коши для системы обыкновенных дифференциальных уравнений вида

, ^ dv . „. . M— + Av = f(v), t> 0;

v = v0, t = 0.

(1.8) (1.9)

Здесь A — сеточный оператор диффузии; M — диагональная матрица масс. Компоненты вектора v включают значения сеточных функций П и п Разобьем вектор v на два подвектора, содержащих две группы переменных: к первой группе отнесем значения, соответствующие сеточным узлам из диффузионной зоны, а ко второй группе — из зоны горения. Тогда представлению вектора-столбца v = (vf, v;f )T (T означает транспонирование) отвечают блочные представления диффузионного оператора (0.6) и вектора правой части

"

0

0

l

"

/ (V) = (/^(V), /2^(1>))т. Отметим, что поскольку в правую часть не входят аппроксимации дифференциальных операторов, /(и) = /(г^), к = 1, 2. Неравенство |А11|(1) ^ ||Н(2) из (0.6) соответствует сильному сгущению сетки в зоне горения. Для введенных обозначений двухуровневая разностная схема решения задачи (1.8), 1.9) имеет вид

M

w?+1 - v?

At

+ An w? + Ai2vn = Л«);

fe+1

M

— V.

w

-

p

p

k

w? + -(w?+1 - w?), k = 1,.. P

p

p

T

+ A>1 p + A22V2 p = /2(V2 p), k = 0,... ,p - 1;

p

(1.10)

(1.11)

(1.12)

? П+1 _ 1 1 1

м 1 Д 1 + 2 Аи«+1 + <) + 2 А^1 + V?) = ¿(ЛК+1) + /1«)), (1.13)

где Д£ = рт, р — некоторое натуральное число; = V0. Как нетрудно видеть, для первой группы переменных (1.10), (1.13) является схемой предиктор-корректор с одним умножением оператора на вектор и одним вычислением правой части на шаге интегрирования. По сути дела речь идет о методе декомпозиции области в следующей редакции. Сначала в диффузионной области делается большой шаг Д£ по явной схеме Эйлера (равенство (1.10)). Затем на границе между зонами диффузии и горения решение линейно интерполируется на р вспомогательных слоев (равенство (1.11)). Далее, в соответствии с равенством (1.12) по явной схеме Эйлера делаются р шагов величины т = Д^/р в зоне горения. И, наконец, в диффузионной зоне решение корректируется в соответствии с равенством (1.13).

1

p

V

2

2. Конечно-элементная дискретизация и стратегия сеточной адаптации

Для пространственной дискретизации задачи (1.1)—(1.7) используется стандартный конформный метод конечных элементов с кусочно-линейной аппроксимацией на триангулированной сетке [15]. При этом осуществляется диагонализация матрицы масс (элементами диагонализованной матрицы являются суммы элементов строки исходной матрицы). Как известно (см., например, [16]), такая процедура не снижает точности вычислений. Отметим, что без диагонализации матрицы масс исчезают все алгоритмические преимущества явных схем. Аналогичным образом аппроксимируется нелинейное слагаемое:

У W(Th, dxdy « mes ^ W(Tj, ni),

Si

где — кусочно-линейная базисная функция — "крышка" с носителем sj С П; ^ — барицентрическое множество [16] сеточного узла (x^y); Tj = Th(xi,yi); n = nh(xj,yj). Здесь и далее Th(x,y) и nh(x,y) — кусочно-линейные функции. В приведенных терминах диагонализованная матрица масс есть M = diag{mes wj}.

Рассмотрим вопрос о построении адаптивной триангулированной сетки. Пусть в области П введена сетка, состоящая из квадратных замкнутых ячеек a(fc) с длинами сторон = 2-kh, где k = 0,... , K. Положим

nhk) = U k = 0,...,K.

Рис. 1. Типы ячеек из множеств

Множество О^ будем называть сеточной областью (или сеткой) к-го уровня. Обозначим = П О®. Пусть введенные множества удовлетворяют следующим требованиям:

к (к)

а) У О^ к=0

О;

б) шв82 (4*'°) =0, к = /;

в) ^ = 0, |к — > 1;

г) содержит не более двух, причем смежных, сторон ячейки а(к) С О^ .

Рассмотрим ячейку а(к) С О(к) такую, что 5а(к) П О(к+1) = 0, т.е. содержащую более четырех сеточных узлов. Множество таких ячеек, образующих переходную зону от множества О(к+1) к множеству О(к), обозначим через О(к'к+1). Из условия г) следует, что

триангуляция ячеек из множества О^ порождает ячейки восьми типов, приведенных

к-1

на рис. 1. Ячейки, не принадлежащие множеству и О(к'к+1), триангулируются диагона-

к=0

лью.

Теперь рассмотрим механизм, в соответствии с которым осуществляется выбор параметра К и реализуется представление а). В основе такого механизма лежит функция плотности сетки Ф(ж,у). Рассмотрим некоторую ячейку а(к) с локально пронумерованными вершинами ^ = (ж^, у^), г = 1,... , 4, и пусть

1

£

г=1

Ф = 1

г=1

Критерием разбиения ячейки а(к) на четыре ячейки (к + 1)-го уровня является условие

|Фа — )| > £1.

(2.1)

Для сетки, которая пересчитывается во времени, есть необходимость объединения четырех ячеек с общей вершиной в одну ячейку (к — 1)-го уровня. Критерием такого объединения является условие

|Фа — Ф(^7)| <£2, 3 = 1,

4.

(2.2)

Экспериментальным путем было найдено оптимальное соотношение параметров £1 и £2, позволяющее наилучшим образом отслеживать при изменении сетки вид функции плотности: £1/£2 = 4. Приведем пример сетки в единичном квадрате О = (0,1)2, построенной описанным выше образом для функции плотности

Ф(ж,у) = шах{е100((х-0-5)2+(у-0-5)2), 0.1 + 0.3 8т4ж} .

(2.3)

4

4

2

а

4

Рис. 2. Триангуляция для функции (2.3). Рис. 3. Функция плотности (2.3).

При этом полагалось к = 1/24, £1 = 10-3. Отметим, что в этом примере построения стационарной сетки мы используем критерий (2.1) только для разбиения ячеек. На рис. 2 приведена полученная сетка, а на рис. 3 — функция плотности. Отметим, что алгоритм "среагировал" на разрыв производной функции (2.3).

Сделаем некоторые замечания относительно построения триангулированной сетки конкретно для задачи (1.1)—(1.7). В качестве функции плотности используем функцию скорости реакции (0.2) на п-м временном слое:

Ф(х,у) = Ж (ТП(х,у),пП(х,у)). (2.4)

При образовании из ячейки а(к) новых четырех ячеек (к + 1)-го уровня ^¿к+1), 3 = 1,... , 4, возникает новый сеточный узел

4

;<к+1) = П <г]к+1).

п

¿=1

Тогда необходимые для дальнейших расчетов значения Т,(г(к+1)) и %(^(к+1)) в новом сеточном узле определяются как средние по значениям этих функций в вершине ячейки а(к). Подчеркнем, что критерием сгущения сетки являются не большой градиент концентрации Пн(х, у), а градиент потока массы, который и определяется главным образом функцией скорости реакции.

И, наконец, укажем способ разбиения области П на подобласти, в соответствии с которым реализуется двухуровневая схема (1.10)—(1.13). Для некоторого К1 < К полагаем

К1

П = иПк), П = П \ п

к=0

В соответствии с указанной декомпозицией области осуществляется разбиение искомого вектора на два подвектора (см. разд. 1).

3. Численные результаты

Изложенную выше методику применим для решения некоторых задач, формулируемых в рамках постановки (1.1)—(1.7). Во-первых, рассмотрим вопрос об экспериментальном выборе параметров дискретизации — параметров пространственной дискретизации к, К и

К1 и дискретизации по времени ДЬ и т (целочисленного параметра р). Подбор параметров будет осуществляться в задаче о распространении фронта пламени без теплопотерь (а0 = а = 0), что фактически соответствует одномерной постановке. Для этой задачи будет приведено сравнение численных результатов по вычислению нормальной скорости пламени ип с результатами, полученными по хорошо известной аналитической формуле Зельдовича. Затем будет рассмотрена задача о зажигании, включающая процесс формирования фронта пламени. И, наконец, будут представлены результаты для распространения фронта неадиабатического горения при наличии теплопотерь, включая критический случай гашения пламени. Все расчеты проводятся при следующих значениях физических параметров: Ь = 0.1 м, / = 0.004 м, а = 8 ■ 10-5 м2/с, Т0 = 300 К, к0 = 1012 с-1, Е/Е = 2 ■ 104 К. Адиабатическая температура ТЬ и коэффициенты теплообмена а0 и а в дальнейшем будут варьироваться.

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

00(х,у) = Ть, п0(х,у) = 1, х < Ь/2, п0(х,у) = 0, х > Ь/2.

В расчетах полагалось К = К(0) = //2. В табл. 1 приведены значения нормальной скорости пламени при К = 6 и К1 = 1 (минимальный пространственный шаг К(6) = / ■ 2-7) и для различных значений шага по времени ДЬ в диффузионной области и величины р, определяющей количество вспомогательных временных слоев в зоне горения. Как видно из таблицы, граница устойчивости близка к границе точности алгоритма: при трех верных знаках в решении при ДЬ = 10-5 с увеличение ДЬ до значения 4 ■ 10-5 с и уменьшение р приводят к потере точности до 5 % с одновременной потерей устойчивости (ДЬ = 4 ■ 10-5 с, р = 8). В табл. 2 приведены значения скорости для ДЬ = 2 ■ 10-5 си р =16 в зависимости от числа уровней пространственной сетки. Очевидно, что пяти пространственных уровней недостаточно для корректного моделирования распространения пламени. В то же время разница между использованием шести и семи пространственных уровней несущественна (~ 1.5%). В итоге устойчивость результатов для заданных параметров физического процесса была достигнута при ДЬ = 2 ■ 10-5 с, р = 16, К = 6. Вообще говоря, можно было бы увеличить значение К1 с целью повышения эффективности алгоритма, но выбор К1 = 1 обеспечивает устойчивость во всем цикле приводимых ниже экспериментов. С этими сеточными параметрами выполнялись все дальнейшие численные расчеты.

Таблица 1. Скорость фронта (зависимость от сетки по времени), м/с

д* 4■10-5 2■10-5 10-5

р = 8 — 0.752 0.794

р = 16 0.761 0.794 0.794

р = 32 0.785 0.794 0.794

Таблица 2. Скорость фронта (зависимость от К), м/с

К 2 3 4 5 6 7

ип 0.230 0.345 0.460 0.503 0.794 0.808

Таблица 3. Нормальная скорость пламени, м/с

Ть 1100 1200 1300 1400 1500 1600 1700

ип 0.102 0.231 0.459 0.794 1.381 2.047 2.986

Пп (3.1) 0.108 0.243 0.488 0.891 1.509 2.404 3.640

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

пп = ^ . (3.1)

Ть — То Е

Отметим, что в отличие от [1] в данной формуле пренебрегается тепловым расширением газа в зоне горения. Сделано это для корректности сравнения с численным расчетом, поскольку в исходных уравнениях (1.1), (1.2) коэффициент температуропроводности предполагается постоянным. Асимптотический анализ этой формулы (малым параметром является величина ДТЬ/Е) показывает, что она дает чуть завышенный результат. Данные численного расчета в зависимости от адиабатической температуры приведены в табл. 3. Как и предполагалось, данные, полученные по формуле (3.1), несколько превышают соответствующие результаты прямого численного моделирования.

Формирование фронта пламени. Целью приводимого ниже численного эксперимента является демонстрация формирования геометрии фронта пламени на начальной стадии при зажигании источником квадратной формы в момент времени £ = 0. При этом полагается ТЬ = 1400К и а0 = а = 5м-1, т.е. рассматривается задача с теплоотводом через "боковые" стенки (у = 0, у = /). Область начального источника тепла есть квадрат

По = [хо — /о, Хо + /о] X [уо — /о, Уо + /о]

с центром в точке жо = 0.09 м, уо = 0.002 м и со стороной 2/о = 0.001 м. При этом начальная сетка адаптировалась к размерам источника в соответствии с принципами, изложенными в разд. 2. В качестве начальных данных рассматриваются функции

0о(ж,у) = Ть, п°(х,у) = 1, (ж,у) € П \ По, П°(х,у) = 0, (ж, у) € По.

При этом поскольку Т = 0 — (Ть — То)п, распределение начальной температуры имеет вид

То(ж, у) = То, (ж, у) € П \ По, То(ж,у) = Ть, (ж, у) € По.

Поле температуры в моменты времени ¿1 = 0, ¿2 = 8 • 10-4 с, ¿3 = 1.6 • 10-3 с, ¿4 = 2.4 • 10-3 с, ¿5 = 4• 10-3 с представлено на рис. 4, где отражен качественный характер распространения пламени. Как хорошо видно, фронт пламени до заполнения всего пространства щели принимает форму окружности, что находится в полном соответствии с принципом Гюйгенса для газофазного горения [1]. После заполнения щели начинается формирование стационарного фронта. При этом пламя распространяется в обе стороны от источника. Влияние теплопотерь сказывается на форме фронта (охлаждение около стенок и как следствие запаздывание зоны химической реакции) и температуре в области прореагировавшего газа

(серый цвет указывает на охлаждение продуктов горения). Отметим, что в области продуктов горения п = 0 и охлаждение происходит за счет уменьшения в по отношению к адиабатической температуре Т,.

Неадиабатическое горение. Дальнейшее развитие описанного выше процесса приводит к стационарному распространению пламени при наличии теплопотерь. В этом случае уже нельзя пользоваться формулой (3.1). Результаты расчетов представлены на рис. 5. Положение фронта приведено в моменты времени ¿6 = 8• 10-3 с, ¿7 = 1.2• 10-2 с, ¿8 = 1.6• 10-2 с, ¿9 = 2 • 10-2 с, ¿ю = 2.4 • 10-2 с. Как уже говорилось, фронт пламени имеет характерный искривленный профиль, обусловленный захолаживанием химической реакции около стенок, через которые в окружающую среду выходит поток тепла. Отчетливо видно понижение

Рис. 4. Формирование фронта пламени. Рис. 5. Распространение фронта пламени.

Рис. 6. Гашение пламени.

температуры в зоне продуктов реакции.

Как известно [1], при повышении интенсивности внешних теплопотерь (увеличении коэффициентов а° и/или од) начиная с некоторых критических значений наступает гашение пламени (так называемый тепловой предел). Аналогичная ситуация имеет место при понижении теплотворной способности химической реакции Q, которое может быть достигнуто уменьшением концентрации горючего реагента в газовой смеси (концентрационный предел). Результаты такого типа вычислительного эксперимента продемонстрированы на рис. 6 в моменты времени ¿1 = 2 • 10-5 c, ¿2 = 4 • 10-3 c, ¿3 = 8 • 10-3 c, ¿4 = 1.6 • 10-2 c, ¿5 = 2.8 • 10-2 c.

В приведенном расчете полагалось а° = а^ = 14 м-1. Для "облегчения" гашения пламени полагалось Tb = 1100 К. При этом инициировался плоский фронт горения

0°(x,y) = Tb, n0(x,y) = 1, x < 0.09м, n0(x,y) = 0, x> 0.09м.

Полное затухание (остаточная температура не превышает 310 К) происходит за 0.08 c.

Таким образом, предложенная численная методика "справилась" со всеми основными процессами, описываемыми модельной задачей (1.1)—(1.7).

Заключение

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

Кратко коснемся возможностей дальнейшего использования многоуровневых явных схем в задачах теории горения. Во-первых, поскольку в энтальпийной постановке энергетический баланс учитывается автоматически, уравнение (1.1) для определения энтальпии (функции cp0) можно аппроксимировать на грубой сетке. При этом можно использовать и неявную схему, соблюдая "осторожность" при реализации краевых условий, зависящих от концентрации. Правда, ограничение на временной шаг At все равно будет определяться локальным условием устойчивости явной схемы для концентрации в диффузионной зоне. Далее можно было бы рассмотреть многозонную задачу (с переходными подобластями от диффузионной зоны к зоне горения). При этом можно использовать уже не двухуровневые, а многоуровневые явные схемы [12]. Но главное — это применение предложенной методики в реальных двух- и трехмерных задачах с учетом внешнего гидродинамического масштаба. И, наконец, отметим широкие возможности для распараллеливания многоуровневых явных схем на многопроцессорных комплексах при решении больших трехмерных задач. Однако следует учесть, что разработка соответствующего высокоэффективного программного обеспечения требует специальных и довольно значительных усилий.

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

[1] ЗЕльдович Я.Б., Баренблатт Г.И., Ливрович В.Б., Махвеладзе Г.М. Математическая теория горения и взрыва. М.: Наука, 1980.

[2] Spalding D.B. Theoretical aspect of flame stabilization // Aircraft. Eng. 1953. Vol. 25, N 295. P. 264-276.

[3] Мержанов А.Г., ХАйкин Б.И., ШкАдинский К.Г. Установление стационарного распространения пламени при зажигании газа накаленной поверхностью // Прикл. механика и техн. физика. 1969. № 5. С. 42-48.

[4] ШкАдинский К.Г. О разностном счете задач зажигания и горения с учетом диффузии и гидродинамики // Физика горения и взрыва. 1969. Т. 5, № 2. С. 264-272.

[5] Демин М.М., Мажукин В.И., Шапранов А.В. Метод динамической адаптации в проблеме ламинарного горения // Журн. вычисл. математики и мат. физики. 2001. Т. 41, № 4. С. 648-661.

[6] Aly S.L., Simpson R.B., Hermance C.E. Numerical solution of the two-dimensional premixed laminar flame equations // AIAA J. 1979. Vol. 17. P. 56-63.

[7] Winters K.H., Rae J., Jackson C.P., Cliffe K.A. The finite element method for laminar flow with chemical reactions // Intern. J. Numer. Methods Eng. 1981. Vol. 17. P. 239-253.

[8] Braack M., Becker R., Rannacher R., Warnatz J. An adaptive finite element method for combustion problems // Proc. 3d Summer Conf. Numer. Modelling in Cont. Mechanics. Prague, 1998. P. 91-100.

[9] Рычков А.Д., ШокинА Н.Ю. Математические модели фильтрационного горения и их приложения // Вычисл. технологии. 2003. Т. 8, спецвыпуск. Ч. 2. С. 124-144.

[10] Braack M. An adaptive finite element method for reactive-flow problems // Ph.D. Thesis. Ruprecht-Karls Univ. of Heidelberg, 1998.

[11] ЛАЕвский Ю.М., Банушкина П.В. Составные явные схемы // Сиб. журн. вычисл. математики. 2000. Т. 3, № 2. С. 165-180.

[12] Banushkina P.V., Laevsky Yu.M. Multi-level explicit schemes and their stability // Russ. J. Numer. Anal. Math. Model. 2001. Vol. 16, N 3. P. 215-233.

[13] ЗоткЕвич А.А., ЛАЕвский Ю.М. Об одном классе двухуровневых явных схем // Сиб. журн. вычисл. математики. 2002. Т. 5, № 2. C. 163-173.

[14] ЗоткЕвич А.А. Численное моделирование распространения фронта пламени в двумерном случае // Тр. конф. молодых ученых. Новосибирск: ИВМиМГ СО РАН, 2004. C. 61-67.

[15] СьярлЕ Ф. Метод конечных элементов для эллиптических задач. М.: Мир, 1980.

[16] ЛАЕвский Ю.М. Метод конечных элементов решения многомерных параболических уравнений. Новосибирск: Изд-во НГУ, 1993.

Поступила в редакцию 1 июня 2006 г., в переработанном виде — 28 сентября 2006 г.