Научная статья на тему 'Численное решение задачи о возникновении верхового лесного пожара в трехмерной постановке'

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

CC BY
274
43
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ДИСКРЕТНЫЙ АНАЛОГ / КОНТРОЛЬНЫЙ ОБЪЕМ / ЛЕСНОЙ ПОЖАР / CONTROL VOLUME / DISCRETE ANALOGUE / FOREST FIRE

Аннотация научной статьи по физике, автор научной работы — Перминов Валерий Афанасьевич

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

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

Numerical solution of crown forest fire initiation problem in three dimensional setting

On the basis of laws of mechanics of continuous reacting media and experimental data the mathematical model of forest fires initiation is developed. The new technique of the numerical decision of spatial (three-dimensional) problems of the theory of forest fires is used. Distributions of velocities, temperatures, concentration a component gas (the oxygen, flying products of pyrolysis and burning) and condensed phases in space in a vicinity of the center of occurrence and distribution of forest fire are received.

Текст научной работы на тему «Численное решение задачи о возникновении верхового лесного пожара в трехмерной постановке»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2009 Управление, вычислительная техника и информатика № 1(6)

УДК 536.46:536

В.А. Перминов ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ О ВОЗНИКНОВЕНИИ ВЕРХОВОГО ЛЕСНОГО ПОЖАРА В ТРЕХМЕРНОЙ ПОСТАНОВКЕ1

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

Ключевые слова: дискретный аналог, контрольный объем, лесной пожар.

В связи с тем, что экспериментальные методы изучения лесных пожаров являются дорогостоящими и не позволяют проводить полное физической моделирование данного явления, представляют интерес теоретические методы исследования. Так, метод математического моделирования позволяет адекватно описывать состояние лесного биогеоценоза и приземного слоя атмосферы при лесных пожарах [1, 2]. С точки зрения механики сплошной среды лес представляет собой некоторый слой многокомпонентной многофазной реакционноспособной сплошной массы, обладающей неоднородными свойствами в вертикальном и горизонтальном направлениях [1]. В рамках физико-математической модели учитываются законы сохранения массы, импульса, энергии, а также физико-химические процессы при пожарах в лесу. Тепло из фронта пожара в несгоревшую зону передается конвекцией, турбулентной теплопроводностью и излучением. Следует отметить, что прямое решение задач по полной модели невозможно. Однако с помощью упрощений, обусловленных свойствами конкретной задачи, удалось получить математические модели, которые позволяют получить решение. Впервые такой подход осуществлен в работах А.М.Гришина и соавторов [1]. Обзор отечественной и зарубежной литературы показал, что в настоящее время по подобным тематикам интенсивно проводятся исследования в США (US Department of Agriculture в рамках федеральной программы по лесным пожарам и Environmental Protection Agency (EPA) по проблемам чистоты воздушной среды), Великобритании и Франции [3 - 7].

1. Физическая и математическая постановка задачи

Как правило, возгорание в лесах происходит в нижнем ярусе леса в напочвенном покрове (опавшая хвоя, мхи, лишайники, отмершая трава и т.д.), а затем огнем охватывается полог леса, то есть образование верхового лесного пожара происходит в результате перехода низового лесного пожара в верховой [8]. Предполагается, что: 1) течение носит развитый турбулентный характер и молекулярным

1 Работа выполнена при финансовой поддержке РФФИ (грант 07-01-96047).

переносом пренебрегаем по сравнению с турбулентным; 2) плотность газовой фазы не зависит от давления из-за малости скорости течения по сравнению со скоростью звука; 3) среда находится в локально-термодинамическом равновесии; 4) известна скорость ветра над напочвенным покровом в невозмущенных условиях; 5) газодисперсная смесь бинарна и состоит из частиц конденсированной фазы, а также газовой фазы - компонентов кислорода, газообразных горючих и инертных компонентов. Здесь очаг низового пожара имеет конечные размеры и над пологом леса задана скорость ветра.

Пусть начало системы координат х1, х2, х3 = 0 расположено в центре источника возникновения горения, в очаге низового лесного пожара. Ось 0х3 направлена вверх, а оси 0x1 и 0х2 - параллельно поверхности земли (ось х1 совпадает с направлением ветра) (рис. 1).

Рис. 1

Сформулированная выше задача сводится к решению следующей системы уравнений:

др + —— (ру,-) = т , 1 = 1,2,3, і = 1,2,3;

5 і д х} 1

скі д Рд _,_, .

Р— = ~^ + ^— (-РУУ1) - Р-С у Іу І -Р8, - ту ; аі д хі д ж,-

1 У

(1)

(2)

СТ д

рср — = — (-рсру’у’Т’) + д5Я5-а„ (Т - Т) + к& (сПк - 4ХаГ); (3)

Сс д -----------------------

Р~а = ^------------(-РУ’іС’а ) + К5а-тса , а = 1,5;

сі д хі

д ж,

(_с_ дП^ Л

3к д х,

V і

- ксПК + 4кГ оТ£ + 4кг аТ4 = 0, к = кг + кГ;

4 д Т

^Р,ср, Фі^Т = 43К3 - 42К2 + к*(сПя - 4аТ54) + ак(Т -^);

і=1

д і

(4)

(5)

(6)

р, =- *, • Р2 дт=- * • Р3 -5т=“с * -1Г * • Р4 %=0; ™

5 / д / д / М1 д /

X Са = 1 ^ =Р*ГXМ-, 17 = (У1,У2’У3^Я = (0,0,Я),

а=1 а=1 М а

Мс

/Й = (1 — ас )*, + *2 +-*3 + *54 + *55 ,

М1

*51 = —*3 — М7 *5, *52 = V Я (1 — ас )*, — *5, *53 = 0,

асУз

*54 =а4*,, *55 =--------------*3.

У3 + У3*

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

Е1 1 и =;, Р ™ Т’—10,5 Г Е2

R = Л1р1ф1ехр |^- —j , R2 = k2p2^2Ts ’ exp ^ r^

R=4зРфзs-0.ехр(-W, j• R>=‘>«2 (Mf) MfT~“5 exp(-RF

где R5a - массовые скорости пиролиза лесных горючих материалов, испарения влаги, горения конденсированных и летучих продуктов пиролиза, образования сажи и пепла и образования a-компонентов газодисперсной фазы; t0 - время зажигания; Opi, р, ф,- - удельные теплоемкости, истинные плотности и объемные доли i-й фазы (1 - сухое органическое вещество, 2 - вода в жидко-капельном состоянии, 3 - конденсированные продукты пиролиза, 4 - минеральная часть, 5 - газовая фаза); Т, Т, - температура газовой и конденсированной фаз; ca - массовые концентрации (a = 1 - кислород, 2 - горючие продукты пиролиза, 3 - инертные компоненты воздуха, 4 - сажа, 5 - пепел); p - давление; UR - плотность энергии излучения; ст -постоянная Стефана - Больцмана; к - коэффициент ослабления излучения; kg, к, -коэффициенты поглощения для газодисперсной и конденсированной фаз; aV - коэффициент обмена фаз; q, Е,, к - тепловые эффекты, энергии активации и пре-дэкспоненты реакций пиролиза, испарения, горения кокса и летучих продуктов пиролиза; , - удельная поверхность элемента лесных горючих материалов; Ма, Мс, М - молекулярные веса индивидуальных компонентов газовой фазы, углерода и воздушной смеси; s, od - удельная поверхность фитомассы и эмпирический коэффициент сопротивления полога леса; с - скорость света; v, - проекции скорости на оси x,; аС, v - коксовое число и массовая доля горючих газов в массе летучих продуктов пиролиза; m - массовая скорость образования газодисперсной фазы; v3* - характерная скорость вдува из очага лесного пожара; a4, a6 - эмпирические константы; g - ускорение свободного падения. Индексы «0» и «е» относятся к значениям функций в очаге горения и на большом расстоянии от зоны пожара соответственно. Верхний индекс «'» относится к пульсационной составляющей данной величины.

Термодинамические, теплофизические и структурные характеристики соответствуют лесным горючим материалам соснового леса и приведены в [2, 9]. Турбу-

лентные потоки тепла массы и количества движения записываются через градиенты среднего течения. Коэффициент межфазного (газ и конденсированная фаза) теплообмена определяется ау = а£ — уСРт, £ = 4ф£ / [5]. Здесь а = Ыык/ё£-

коэффициент теплообмена для элемента ЛГМ (например, хвоинки), Ыы - число Нуссельта для цилиндра, X - коэффициент теплопроводности для хвоинки; у - параметр, характеризующий отношение между молекулярной массой окружающих и вдуваемых газов. Таким образом, совокупность уравнений (1) - (7) с соответствующими начальными и граничными условиями, являются балансовыми соотношениями массы, энергии и количества движения, представляющие собой постановку задачи, решение которой позволяет определить характеристики сложного взаимосвязанного процесса возникновения лесного пожара. Начальные и граничные условия будут иметь следующий вид

t = °: у = ^ У2 = 0 у = 0 Т = Те ,т£ = Те, Са = Са « ,фг = фге ; (8)

1 =— *1« : V = V,, У2 = 0, = 0, Т = Те, Са= Са,, ——+ си*/2 = 0; (9)

д *1 3к д *1

д VI д у9 д У3 д Са д Т с д и* с

* = * :^ = 0, _2 = 0,—3 = 0, —“ = 0,-----------------------= 0,-------* + -и* = 0 ; (10)

д *1 д *1 д *1 д *1 д *1 3к д *1 2

= *20:^ = 0, = 0, ^ = 0,££а = ^ = 0,—_1 ££* + -и* = 0 ; (11)

А ¿.У) --V -'--V -'--V -'--V ^ Л 7 Л 7 V Г

д *2 д *2 д *2 д *2 д *2 3к д *2 2

*2 = *2,:-дЯ- = 0, ^ = 0,±1 = 0,££а = 0,^ + Си* = 0; (12)

¿е л -'-л -'-л -'-л -'-л ^ л / л \ у

д *2 д*2 д *2 д *2 д *2 3к д *2 2

*3 = 0 :у = 0, у2 = 0,дСа = 0, —ди* + -и* = 0 ,

д *3 3к д *3 2 (13)

К30

,Т = ТЯ при |*1 < Д,|*^ < Д и У3 = 0,Т = Те при |*1 > Д,|*2| > Д;

*3 = *3,:^ = 0, ^ = 0,= 0,^ 0^1 = 0,—ди*+-ий = 0. (14)

д *3 д*3 д *3 д *3 д *3 3к д *3 2

2. Методика решения

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

д-(рф) + д^(р«гФ) = д^^дф) + £ , ¿=1,2,3. (15)

^ д*г д*г У д*г )

где t, *г - временная и пространственные координаты, р - плотность, - компо-

ненты вектора скорости, Г - коэффициент переноса (например, Г — коэффициент диффузии), £ - источниковый член. В частности в £ может входить приток (сток) тепла за счет химических реакций в уравнении энергии или увеличение (уменьшение) концентраций компонент в результате химических реакций в уравнении диффузии. Конкретный вид Г и £ зависит от смысла переменной Ф (в действи-

тельности следовало бы использовать обозначения ГФ и 5ф, но это привело бы к слишком большому количеству нижних индексов в дальнейших выкладках). В (15) также подразумевается суммирование по индексу i. Начальные и граничные условия для уравнения (15) в общем виде записываются следующим образом:

дФ

Ф 1=0 = Фе , a1 X = x„ +а2Ф X=х„ = а3 . (16)

д xi e

Здесь а - функция аргументов t и xi, xie - значение аргумента xi на соответствующей границе расчетной области. Если в рассматриваемой математической модели учитываются химические реакции, то обычно дополнительно используются обыкновенные дифференциальные уравнения первого порядка или, как говорят, уравнения химической кинетики, которые в общем виде можно записать следующим образом:

дФ

— = -K(t)Ф + F(t), K(t) > 0, F(t) > 0, Ф t=o = Фe. (17)

dt

Построение дискретного аналога для уравнения вида (15) с соответствующими краевыми условиями типа (16) проведем на основе метода контрольного объема, подробное изложение которого дается в [10]. Основная идея метода поддается прямой физической интерпретации. Расчетная область разбивается на некоторое число непересекающихся контрольных объемов таким образом, что каждая узловая точка содержится в одном контрольном объеме. Дифференциальные уравнения интегрируются по каждому контрольному объему. Для вычисления интегралов используются кусочные профили, которые описывают изменение Ф между узловыми точками. Полученный таким образом дискретный аналог выражает закон сохранения Ф для конечного контрольного объема точно. Точно так же, как дифференциальное уравнение выражает закон сохранения для бесконечно малого контрольного объема. Важнейшее свойство метода контрольного объема состоит в том, что при его использовании точно выполняются интегральные законы сохранения таких величин, как масса, количество движения и энергия, в каждом контрольном объеме и для любой группы контрольных объемов и, следовательно, на всей расчетной области [10]. Таким образом, даже решение на грубой сетке удовлетворяет точным интегральным балансам во всей области, то есть дискретный аналог (разностная схема) Патанкара - Сполдинга является консервативным [10 - 12]. В результате дискретизации дифференциальной задачи получаем системы сеточных уравнений для каждого дифференциального уравнения. Причем матрица каждой из системы является семидиагональной. Для разрешения полученных систем алгебраических уравнений используем метод SIP [13].

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

функций, вычислялась невязка уравнений, которая затем трактовалась как фиктивный источник в каждом уравнении. Затем значения функций восстанавливались. Точность восстановления функций составляла не менее 0,5 %. Устойчивость и точность полученных решений проверялась также с помощью уменьшения шагов по времени и пространству. При расчетах использовался алгоритм автоматического выбора шага по времени. Согласование полей скорости и давления осуществлялось итерационным образом в рамках алгоритма SIMPLE [10].

3. Результаты расчетов

На основе представленной математической постановки (1) - (14) проводились численные расчеты по определению картины процесса возникновения верхового лесного пожара в результате зажигания полога леса от заданного очага горения. В результате численного интегрирования получены поля массовых концентраций компонент газовой фазы, температур, объемных долей компонентов твердой фазы в различные моменты времени. На основе полученных данных следует, что с течением времени возрастают температуры газовой и твердой фаз, происходит уменьшение массовой концентрации кислорода и изменение количества горючих продуктов пиролиза и объемных долей фаз на нижней границе полога леса вблизи очага горения. В результате воздействия очага повышенной температуры в его окрестности происходит прогрев полога леса, испарение влаги и разложение сухого ЛГМ. В результате этого в пологе леса выделяются летучие горючие продукты пиролиза. Во все время процесса температура газовой фазы выше температуры твердой фазы. Газообразные продукты пиролиза, выделившиеся в результате разложения ЛГМ, воспламеняются в пространстве между пологом леса и напочвенным покровом. При этом также происходит уменьшение концентрации кислорода. Начиная с момента зажигания, температуры газовой и конденсированной фаз становятся одинаковыми. То есть полученные результаты показывают, что процесс зажигания носит многостадийный характер. Таким образом, под воздействием очага повышенной температуры происходит воспламенение полога леса. В процессе прогрева под воздействием архимедовой силы и ветра происходит всплытие и перенос в горизонтальном направлении нагретых газов и продуктов пиролиза и горения из очага горения по направлению ветра. На рис. 2 представлены распределения поверхностей равной температуры (изоповерхностей) (1 - T =1,1; 2 - T = 1,5; 3 - T = 2; 4 - T = 3.) газовой фазы в различные моменты времени: а - t = 3 и б - t = 5 с соответственно), а также векторные поля скорости в области воспламенения полога леса от очага горения.

Из рисунков видно, что в результате воздействия очага низового лесного пожара произошло воспламенение полога леса и область горения продвигается в глубь. Кроме того, видно, что с течением времени за счет зоны горения происходит прогрев полога леса в вертикальном и в большей степени в горизонтальном направлении. В результате этого происходит сушка и начинается процесс пиролиза ЛГМ, тем самым будет обеспечиваться дальнейшее распространения лесного пожара по пологу леса.

На рис. 3 представлено распределение в области горения поверхностей равной концентрации кислорода в моменты времени t = 3 и 5 с соответственно (1 - C1 = 0,9;

2 - 0,5; 3 - 0,2). В результате возникновения горения происходит уменьшение массовой концентрации кислорода, который расходуется на окисление (горение) летучих горючих продуктов пиролиза.

Рис. 2. Распределение поля скорости и поверхностей равной температуры

х2, м

Рис. 3. Распределение поверхностей равной массовой концентрации кислорода

Распределение поверхностей равной массовой концентрации летучих продуктов пиролиза (в первую очередь это СО, которые распространяются по направлению ветра, представлено на рис. 4 (1 - С2 = 0,3; 2 - 0,5; 3 - 1) в различные моменты времени: а - t = 3 и б - t = 5 с соответственно).

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

х2, м

Рис. 4. Распределение поверхностей равной массовой концентрации

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

Заключение

Полученные численные характеристики процесса зажигания полога леса от очага низового лесного пожара согласуются с результатами экспериментальных исследований [1, 2, 14] и данными расчетов, представленных в [9].

Автор выражают благодарность профессору Томского госуниверситета д.ф.-м.н. А.М. Гришину за разработку общей математической модели лесных пожаров, на основе которой получена постановка задачи представленная в данной работе.

ЛИТЕРАТУРА

1. Конев Э.В.Физические основы горения растительных материалов. Новосибирск: Наука, 1977.

2. Гришин А.М. Математические модели лесных пожаров и новые способы борьбы с ними. Новосибирск: Наука, 1992.

3. Albini F. An overview of research on wildland fire // Proc. of the 5th Intern. Symp. on Fire Safety Science, Intern. Associat. for Fire Safety Science, Melbourn, Australia, 3 - 7 March, 1997. P. 59 - 74.

4. Adams J.S., Williams D. W., Tregellas-Williams J. Air velocity, temperature and radiant -heat measurements within and around a large free burning fire // 14th Intern. Symp. on Comb. Pittsburgh, 1972. P. 1045 - 1052.

5. Palmer T.Y. Convection columns above large experimental fires // Fire Technology. 1975. V. 11. No. 2. P. 11 - 118.

6. Stocks B.J., Alexander M.E., Lanoville R.A. The Intern. Crown Fire Modeling Experiment -an update // International Forest Fire News. 1999. V. 21. P. 89 - 90.

7. Alexander M.E., Stocks B.J., Wotton, B.M., et al. The international crown fire modeling experiment: an overview and progress report // Symp. on Fire and Forest Meterology, Phoenix, AZ, 11-16 January 1998. American Meteorological Society, Boston, MA, 1998.

8. Гришин А.М., Зятнин В.И., Перминов В.А. Экспериментальное исследование перехода низового лесного пожара в верховой. Томск: Томский госуниверситет / Деп. ВИНИТИ № 982-91 от 06.03.91. 22 с.

9. Гришин А.М., Перминов В.А. Переход низового лесного пожара в верховой // Физика горения и взрыва. 1990. Т. 26. № 6. С. 27 - 35.

10. Патанкар С.В. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.

11. Harlow F.N., Welch J.E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface // Phys. Fluids. 1965. V. 8. P. 2182.

12. Patankar S. V., Spalding D.B. A calculation procedure for heat, mass and momentum transfer in three dimensional parabolic flows // Intern. J. Heat Mass Transfer. 1972. V. 15. P. 1787.

13. Stone H.L. Iterative solution of implicit approximations of multi-dimensional partial differential equations // SIAM J. Num. Anal. 1968. No. 5. P. 530 - 558.

14. Леонтьев А.К., Моршин В.Н. Метод расчета воспламенения тонкой растительной частицы в конвективном потоке газа // Интенсификация лесозаготовительных и лесохозяйственных производств. Л.: ЛТА, 1989. С. 59 - 67.

Статья представлена кафедрой прикладной математики факультета прикладной математики и кибернетики Томского государственного университета и оргкомитетом VII Всероссийской научно-практической конференции с международным участием «Информационные технологии и математическое моделирование». Поступила в редакцию 21 ноября 2008 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.