CC BY
47
Вычислительные технологии
Том 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 г.
