Научная статья на тему 'Математическое моделирование распространения плоского фронта лесного пожара'

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

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

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

We propose a mathematical formulation of the problem of the forest fire initiation and propagation based on the averaged three dimensional non-stationary mathematical model of forest fire. A method of control volume is used to obtain the discrete analogies. The boundary-value problem is solved numerically. The results of mathematical modeling allow investigation of the initiation and spreading dynamics of the forest fire.

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

Текст научной работы на тему «Математическое моделирование распространения плоского фронта лесного пожара»

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

Том 11, Специальный выпуск, 2006

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ПЛОСКОГО ФРОНТА

ЛЕСНОГО ПОЖАРА

В. А. Перминов Беловский филиал Кемеровского государственного университета, Россия e-mail: p_valer@mail.ru

We propose a mathematical formulation of the problem of the forest fire initiation and propagation based on the averaged three dimensional non-stationary mathematical model of forest fire. A method of control volume is used to obtain the discrete analogies. The boundary-value problem is solved numerically. The results of mathematical modeling allow investigation of the initiation and spreading dynamics of the forest fire.

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

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

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

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

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

к

J ф ¿х3 = ф к. о

Здесь ф — некоторое среднее значение величины ф.

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

§7 + 1Г-(РЧ)=™> 3 = 1Л г = 1,2; (1)

дЪ дхз

¿V, дР д - -

Р~ГГ = - п--^ ~ р8СаУг\у\ ~ р^Л ~ ГПУЦ (2)

аъ дх,- дхз

¿Т д _

рсР-^ =-^{-рсру]Т') +ау{Т-Т8) + - + кд{с11К - 4стТ4); (3)

з

_

= + Къа ~ + ^ + а = 1,5; ^

д

дх \3к дх

- » - ксик + 4квоТ! + АкдаТА + - д+)/И = 0, к = кд + к8] (5)

4 дТ^

^Г РгСргРг— = фДз ~ + ^(сГ/д - 4(тТ^) + 0!у(Т - Т<?); (6)

¿=1

= = РЗ^ = асД1 - = 0; (7)

5 5

С,

$> = 1,рв = = ^ = (0,0,0),

а=1 а=1

ш = (1 - ас)Д1 + Д2 + ^Дз + Д54 + Д55,

М1

Дб1 = —-ЙЗ — „ л ! Дб; Дб2 = — СКс)Д1 — Дб, ДбЗ = О, Д54 = СК4Д1, Д55 = -—-—Дз-

2М2 Уз + Уз*

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

Д1 = кгрг^г ехр ( ) , Д2 = к2р2<-р2Т~0-5 ехр (— Е2

ДТЧ / ' ' 5 1 I ДТ

/ЕЛ , 025 с2м _2 25 / е5

Дз = ехр (-— ^ , Д5 = к5м2 ^ —Т- • ехр (- —

°.25

Начальные и граничные условия имеют следующий вид:

г = 0 : У1 =0 , У2 =0 , Т = Те , Са = Сае , Т = Те , ^ = ^ (8)

= -Ж1е : VI = Уе, у2 = 0 , Т = Те, са = сае, - + = 0; (9)

Х1 Х1е : д у1 дж1 = 0, дь2 дж1 = 0, ду3 дж1 д Са дж1 дТ дж1 с диД стт = 0, + = 3к дж1 2

ж2 = ж2 0 : д у1 дж2 = 0, дь2 дж2 0, дь3 дж2 дЖ2 дЖ2 0, + = 0 3к дж2 2

Х2 = Ж2 е д у1 ' дх2 = 0, ду2 дж2 = 0, <9У3 дж2 д Са дж2 о дт дх2 ^ с д11я ^ с ф ' 3к дж2 2

(12)

Значения функций в очаге зажигания на напочвенном покрове (|ж11 < Ах, |ж2| < Ау) задаются в зависимости от времени:

Те + ^ (То - те), г < ¿0,

рУз = т , Т = Т = ' ^

Те + (Т° - Те) ехр

-м - -1

, г > г°.

В приведенных уравнениях Я\ — Я5, Я5а — массовые скорости пиролиза лесных горючих материалов, испарения влаги, горения конденсированных и летучих продуктов пиролиза и образования «-компонентов газодисперсной фазы; £0 — время зажигания; срг, рг, фг — удельные теплоемкости, истинные плотности и объемные доли г-й фазы (1 — сухое органическое вещество, 2 — вода в жидко капельном состоянии, 3 — конденсированные продукты пиролиза, 4 — минеральная часть, 5 — газовая фаза); Т, Т8 — температуры газовой и конденсированной фаз; са — массовые концентрации (а =1 — кислород, 2 — горючие продукты пиролиза, 3 — инертные компоненты воздуха, 4 — сажа, 5 — пепел); р — давление; ид — плотность энергии излучения; а — постоянная Стефана — Больцмана; к — коэффициент ослабления излучения; кд,к8 — коэффициенты поглощения для газодисперсной и конденсированной фаз; ау — коэффициент обмепа фаз; дг, Ег, кг — тепловые эффекты, энергии активации и предэкспоненты реакций пиролиза, испарения, горения кокса и летучих продуктов пиролиза; 5 — удельная поверхность элемента лесных горючих материалов; Ма, Мс, М — молекулярные веса индивидуальных компонентов газовой фазы, углерода и воздушной смеси; 5, са — удельная поверхность фитомассы и эмпирический коэффициент сопротивления полога леса; с — скорость света; уг — проекции скорости на оси жг; ас, V — коксовое число и массовая доля горючих газов в массе летучих продуктов пиролиза; т — массовая скорость образования газодисперсной фазы; — характерная скорость вдува из очага лесного пожара; а4, а5 — эмпирические конетанты; д — ускорение свободного падения. Индексы "О" и "в" относятся к значениям функций в очаге горения и на большом расстоянии от зоны пожара соответственно. Штрих относится к пульсационной составляющей данной величины. Термодинамические, теплофизические и структурные характеристики соответствуют лесным горючим материалам соснового леса и приведены в [2], Турбулентные потоки тепла массы и количества движения записываются через градиенты среднего течения на основе подхода, предложенного в [2], В системе (1)-(7) величины С, тг, За., характеризуют обмен импульсом, массой а-компонента и энергией как с приземным слоем атмосферы, так и с нижним ярусом леса и определяются соответствующими граничными условиями.

Таким образом, совокупность уравнений (1)-(7) является балансовыми соотношениями массы, энергии и количества движения и представляет собой постановку сопряженной задачи, решение которой позволяет определить характеристики сложного процесса распространения лесного пожара. Решение сформулированной выше задачи связано со значительными математическими трудностями. Упомянутые в нем величины С, тг, За, вообще говоря, должны определяться в процессе решения сопряженной задачи. Предполагается, что тепло- и массообмен фронта пожара со слоем атмосферы описывается законом Ньютона, и соответствующие величины в (1)-(7) определяются подобным образом, как и

в [2].

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

кя = ая<ря/¿я , кд = ао^о/dо + аzу 2/¿2, = + + ^з + •

Здесь ая — коэффициент поглощения элемента лесных горючих материалов (ЛГМ); ¿0, — диаметры элемента ЛГМ, частиц сажи и дыма (в виде сферы); а0, — коэффициенты поглощения отдельных частиц сажи и дыма.

Коэффициент межфазного (газ и конденсированная фаза) теплообмена определяется как ау = аБ — 7СРт, Б = Здесь а = 1NuX/ds — коэффициент теплообмена для

элемента ЛГМ (например, хвоинки); Nu — число Нуссельта для цилиндра; А — коэффициент теплопроводности для хвоинки; y — параметр, характеризующий отношение между молекулярной массой окружающих и вдуваемых газов.

Система уравнений (1)-(7) с начальными и граничными условиями (8)—(12) для численного интегрирования редуцирована к дискретной форме с помощью метода контрольного объема Патанкара — Сполдинга [6]. При получении дискретного аналога использовалась так называемая шахматная сетка, т, е, компоненты скорости рассчитывались в точках, расположенных на гранях, а скалярные функции — в центрах контрольных объемов, В этом случае только физичные поля скорости могут удовлетворять уравнению неразрывности, Другое важное преимущество шахматной сетки заключается в том, что разность давлений между двумя соседними узловыми точками определяет составляющую скорости в точке, расположенной между этими узловыми точками [6],

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

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

На основе изложенной математической постановки (1)—(12) проводились численные расчеты по определению картины процесса возникновения верхового лесного пожара в результате зажигания полога леса от заданного очага горения. Инициирование горения осуществлялось с помощью задания в прямоугольной области температуры горения в зависимости от времени, которая характерна для фронта горения. Данный очаг существовал некоторое время до момента зажигания полога леса, а затем отключался, В результате численного интегрирования получены поля массовых концентраций компонентов газовой фазы, температур, объемных долей компонентов твердой фазы.

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

меняются в пространство между пологом леса и напочвенным покровом (приблизительно в момент времени 2,5 с), При этом также происходит уменьшение концентрации кислорода. Начиная с момента зажигания температуры газовой и конденсированной фаз становятся одинаковыми (рис, 2,а).

Рис. 2. Распределение с течением времени температур газовой и конденсированной фаз, концентрации кислорода и продуктов пиролиза и объемных долей фаз: а — кривая 1 — T = T/Te, кривая 2 —Ts = Ts/Te,Te = 300К; кривая 1 — ci, кривая 2 — c2, ca = ca/cie; в — кривая 1 — pi = pi/pie, кривая 2 — ф2 = Р2^2/Рс, кривая 3 — фз = рзфз/acPiVie-

Рис. 3. Распределение изотерм и поля скорости в различные моменты времени: Ь = 4.3с (а), Ь = 6с (б), Ь = 8с (в); V = 5м/с; Т = 2, 2.5, 3, 4 (кривые 1-4 соответственно); Т = Т/Те, Те = 300 К.

Таким образом, под воздействием очага повышенной температуры воспламеняется полог леса. При дальнейшем прогреве полога леса происходит продвижение зоны воспламенения в глубь полога леса в направлении ветра. На рис, 3-5 представлено распределение температур, массовых концентраций кислорода и летучих продуктов пиролиза в различные моменты времени в процессе распространения фронта лесного пожара. Из рисунков следует, что формируется фронт горения, который распространяется по лесному массиву, На рис, 3 показано распределение изотерм газовой фазы в три различных момента времени.

Перед областью горения формируется фронт пиролиза (см. рис. 5) в виде зоны с повышенным содержанием летучих горючих продуктов пиролиза. Под влиянием горения происходит также выгорание последних, что вызывает уменьшение концентрации кисло-

го'-'-'-1---1-1-^-1-

О 10 20 30 40 30 40 50

Рис. 4. Распределение кислорода С1 в моменты времени: £ = 4.3с (а), £ = 6с (б), £ = 8с (в); Уе = 5м/с; С1 = 0.9, 0.8, 0.5, 0.4 (кривые 1-4 соответственно).

О 10 го 20 Щ 10 30 ДО 50 60 АрМ

Рис. 5. Распределение газообразных продуктов пиролиза С2 в моменты времени: £ = 4.3 с (а) £ = 6с (б), £ = 8с (в); Уе = 5м/с; с2 = 0.1, 0.5, 1 (кривые 1-3 соответственно).

рода (см, рис, 4), При дальнейшем расчете наблюдается перемещение изотерм и изолиний компонентов газовой фазы, что указывает на перемещение и расширение фронта лесного пожара. Расчеты показывают, что с увеличением скорости ветра возрастает скорость распространения горения. Увеличение влагосодержания ЛГМ ведет к уменьшению скорости распространения фронта пожара, а затем и к его прекращению, что соответствует результатам работ [1, 2] и результатам экспериментальных исследований, представленных

в [7-9].

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

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

[2] Гришин A.M., Грузин А.Д., Зверев В.Г. Математическая теория верховых лесных пожаров // Теплофизика лесных пожаров / НТФ СО РАН СССР. Новосибирск, 1984. С. 38-75.

[3] Перминов В. А. Математическое моделирование возникновения массовых и верховых лесных пожаров с учетом радиационно-конвективного тепломассопереноса и двухтемпературности среды: Дис. ... канд. физ.-мат. наук. Томск, 1995.

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

[5] Perminov V. A. Mathematical modeling of crown forest fire initiation // III Intern. Conf. on Forest Fire Research and 14th Conf. on Fire and Forest Meteorology. Portugal, Luso, 1998. P. 419-431.

[6] Патанкар C.B. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.

[7] Исаков Р.В. Воспламенение хвои при развитии низовых пожаров в верховые: Дис. ... канд. техн. наук: Красноярск, 1985.

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

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

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

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