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

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

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 с соответственно).

х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 г.