А. Н. СУББОТИН, д-р физ.-мат. наук, профессор Энергетического института Национального исследовательского Томского политехнического университета, г. Томск, Россия
УДК 536.24:532.546.2
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ НИЗОВОГО ЛЕСНОГО ПОЖАРА ПО ОПАДУ ХВОИ С ПРОТИВОПОЖАРНЫМИ РАЗРЫВАМИ
Рассматривается математическая модель, описывающая процессы зажигания и распространения волны горения по пористому слою лесных горючих материалов, имеющих противопожарные разрывы. Процессы тепломассообмена, протекающие в факеле пламени, и его взаимодействие с горящим слоем рассчитываются по эмпирическим зависимостям. Показано, что предложенная математическая модель с достаточной для практики точностью моделирует распространение волны горения по плотному слою лесных горючих материалов (опад хвои, подстилка и т. д.) и позволяет рассчитывать размеры противопожарных разрывов. Получены зависимости длины факела и ширины противопожарного разрыва от запаса горючего материала, его влаго-содержания и условий тепломассообмена лесного горючего материала с приземным слоем воздуха.
Ключевые слова: опад хвои; низовой пожар; противопожарный разрыв; факел пламени; вла-госодержание.
Физико-математическая постановка задачи
Рассмотрим зажигание и горение слоя лесного горючего материала (ЛГМ). Будем считать, что слой ЛГМ пористый и имеет достаточно высокую плотность, так что можно использовать уравнение движения в форме квадратичного закона фильтрации. При исследовании процесса распространения волны горения будем учитывать испарение влаги, пиролиз ЛГМ, горение кокса, образование золы и газообразных продуктов пиролиза и горения. Предположим, что пожар инициируется плоским источником тепла, а возникший очаг распространяется как плоская волна горения в направлении оси х, совпадающей с направлением ветра.
Рис. 1. Схема низового пожара: 1 — слой ЛГМ; 2 — факел пламени; х, хp — координаты противопожарной полосы разрыва
Схематически картина процесса после образования фронта горения изображена на рис. 1. Согласно опытным данным следует учитывать процессы, протекающие внутри слоя ЛГМ, такие как сушка, пиролиз и гетерогенная реакция окисления кокса. Процессы же, связанные с газообразными продуктами пиролиза, не учитываем, поскольку они сгорают в факеле над слоем ЛГМ. Значения параметров, характеризующих структуру слоя, таковы, что характерное расстояние между элементами мало по сравнению с размерами горящего слоя, поэтому для описания переноса энергии излучением внутри слоя ЛГМ использовалось приближение лучистой теплопроводности [1].
Усредняя систему уравнений [2] по вертикальной переменной 2, с учетом граничных условий на верхней и нижней поверхностях ЛГМ, получим одномерную математическую модель, описывающую тепло-и массообмен в слое в случае низового пожара:
Р1
5ф1
= -Яг
Р 5ф2 = Я Р 2 = 'Я
21;
дг р' Г2 дг Яр = кр Р1Ф1 ехР(- Ер/(ЯТ));
Я21 = к2* Р 2Ф2 ехР (- Е2*1(ЯТ));
дф з
Рз = ъсЯ„ - БЯз*;
дг
Р 4
дф 4 чт-
= п*яР; ефг =1;
5
Е
г = 1
(1)
(2)
© Субботин А. Н, 2012
(Р5 Ф5) + Ф5 и) =
= + я2р + 5я3р - (Р5Ф5")№ ; (3)
н
о о
(Р5 Ф5 Са ) + ОХ (Р5 Ф5 иС« ) =
о Г ^ Зс„ ^ Я5аw ...
= ЗХ 1 Р5 Ф5 Яр + 5Я5а,--; (4)
а
= !>6; X са = !;
а =1
Я5аw = а1 Ф5 (са Сае VСр5 ; яз, = (М3/М51) Р5 Ф5 С1 к35 Ф3 ехр(- Е3,/(ЯТ));
ОТ ОТ О Л. , .ОТ^
X РФ. СР1 + Р5 Ф5 Ср5 и = ^ I^* + ^г) зх I+
О £ р Ох Ох\
+ 935 - 92,Р2 5Я2, -
-^1(Т - Те)(Т - Т,);
ч+ = (1 -ф5) е/ аТ/ х/; ч- = (1 -ф 5) ^ а Т1; к
ц + кРр51 и | Ох ' М = 1
Ор Р5 ЯТ
г ■ р = 5 ■
М
X —
а = 1 М„
(5)
(6)
Здесь р15 р2, Р3, р4, р5 — истинные плотности соответственно ЛГМ, воды, кокса, золы и газовой фазы, кг/м3; ф1, ф2, ф3, ф4 — объемные доли соответственно ЛГМ, влаги, кокса и золы; £ — время, с; Яр — массовая скорость реакции пиролиза ЛГМ, кг/(м3 с); Яъ, Я3, — массовые скорости соответственно испарения воды и окисления кокса, кг/(м2-с); кр — пред-экспоненциальный множитель реакции пиролиза ЛГМ, с-1; к2,, к3Б — предэкспоненциальные множители соответственно испарения влаги из слоя и реакции окисления коксового остатка, м/с; Ер, Еъ, Е3Б — энергии активации соответственно реакции пиролиза ЛГМ, испарения влаги из слоя и реакции окисления коксового остатка, Дж/моль; Я — универсальная газовая постоянная, Дж/(моль К); Т — температура, К; п яС, я,, п^ — массовые доли соответственно газообразных продуктов, кокса, золы и а-компоненты газовой фазы, образующихся при разложении ЛГМ; 5 — удельная площадь поверхности пор, м-1; ф5 — объемная доля газовой фазы (пористость); и, V — компоненты вектора скорости, м/с; Са — массовая концентрация а-компоненты газовой фазы (индекс а = 1 соответствует 02, а = 2
— СО, а = 3 —С02, а = 4 — Н2, а = 5 —СН4, а = 6
— К2, а = 7 — Н20); Ба — эффективный коэффициент диффузии а-компоненты газовой фазы, м2/с;
Я5а, — массовая скорость образования (исчезновения) а-компоненты газовой фазы вследствие протекания гетерогенных реакций внутри слоя ЛГМ, кг/(м2-с); Я5аи, — скорость изменения массы а-компоненты на верхней поверхности слоя ЛГМ, кг/(м2-с); Н — высота слоя ЛГМ, м; а1, а2 — коэффициенты теплоотдачи, Вт/(м2-К); а 1 = а 0 [1 + к0(р5ф5V)] ; а0 — коэффициент теплоотдачи, определяющий теплообмен непроницаемой стенки с потоком газа; к0(р5ф^)„, учитывает изменение теплообмена на верхней поверхности вследствие выдува образующихся газообразных продуктов через пористую стенку; к0 — эмпирическая константа, м2-с/кг; Сае—массовая концентрация а-компоненты в окружающей среде; Ср1, Ср2, Ср3, Ср4, Ср5 — удельные теплоемкости соответственно ЛГМ, воды, кокса, древесной золы и газовой фазы, Дж/(кг К); М3 — атомарная масса углерода; М51 — молекулярная масса кислорода, кг/кмоль; С1 — массовая концентрация кислорода; X,, Хг — коэффициенты соответственно кондуктивной и лучистой теплопроводности, Вт/(м-К); X г = 16а Т3/5; а — постоянная Стефана - Больцмана, Вт/(м2-К4); X — эффективный коэффициент теплопроводности
слоя ЛГМ; X =Х, +Х г = X X.
ф. + Х г ; X1, X2, Х3'
. = 1
X4, X5 — коэффициенты теплопроводности соответственно ЛГМ, воды, кокса, золы и газовой фазы;
X 5 = X X
5а са; ^а
коэффициент теплопровод-
а=1
ности а-компоненты газовой фазы; 4 чъ, Чъ — тепловые эффекты реакций соответственно разложения ЛГМ, испарения влаги и окисления кокса, Дж/кг; ч+, 9- — лучистый тепловой поток от факела пламени, соответственно падающий на верхнюю поверхность и исходящий от поверхности ЛГМ, Вт/м2; Т, Те — температура соответственно почвы и внешней среды над слоем ЛГМ, К; е^ — излуча-тельная способность факела; Т^ — средняя температура в факеле пламени, К; — угловой коэффициент, который определяет долю энергии от факела пламени, приносимую в рассматриваемую точку; е„, — интегральная степень черноты слоя горючего материала; к — коэффициент проницаемости пористой среды, м2; ц — коэффициент динамической вязкости, Пас; Р — эмпирический коэффициент в квадратичном законе сопротивления; р — давление в макропорах ЛГМ, Па; М — молекулярная масса газовой фазы, кг/кмоль; М5а — молекулярная масса а-компоненты газовой фазы, кг/кмоль; индексы "е" и характеризуют соответственно внешнюю среду и верхнюю поверхность слоя лесного горючего материала.
Предполагая, что теплопроводность непроницаемого инертного источника зажигания значительно больше теплопроводности слоя ЛГМ, пренебрегаем распределением тепла в источнике. Тогда началь-
ные (7) и граничные (8) условия для системы уравнений (1)-(6) можно записать в виде:
г = 0:
Т|;
= Т
фг
0:
и = 0
са сан> а = 1, 6, + |х>х0 --н>
Т1х = х0 = Т0, Р 5 = Р 5н, и = 0
фг = фгр для XI < х < Хр, _
фгн для XI < X < Хр, г = 1, 4 ;
(7)
х = Ь:
Здесь са
д са
дх
са са
„ , дТ = ¿Т
= 0, ~дх _Р о Со хо "¿7 '
Т = Тн.
р = Рн.
(8)
_ан - начальное значение массовой концентрации а-компоненты газовой фазы в порах ЛГМ; х0 — ширина источника зажигания, м; Тн — начальное значение температуры ЛГМ, К; Р5н — начальное значение плотности газовой фазы в порах, кг/м3; фгр (г = 1,4) — значения объемных долей соответственно ЛГМ, влаги, кокса и золы в зоне разрыва; х1, хр — левая и правая координаты противопожарного разрыва, м; фгн (г = 1,4) — начальные значения объемных долей соответственно ЛГМ, влаги, кокса и золы; р0 —плотность источника зажигания, кг/м3; с0 — теплоемкость источника зажигания, Дж/(кг-К); Ь — длина исследуемой полосы лесного горючего материала, м;рн — начальное значение давления газовой фазы в порах, Па.
Известно, что лучистый поток от факела пламени является определяющим в процессе горения, поэтому его необходимо учитывать. Для упрощения задачи процессы, протекающие в газовой фазе над слоем ЛГМ (в пограничном слое), рассчитывались не по математической модели, а по известным экспериментальным зависимостям. В частности, высота факела при отсутствии ветра согласно [3] определялась из зависимости
= ^
тз и„
(9)
где В — эмпирический множитель, равный: 10 — для живой хвои, 9 — для опада хвои, 8 — для листьев
и травы, м0'5- с
°'5/кг°'5;
тз — запас ЛГМ (масса 1 м площади ЛГМ), кг;
— скорость распространения фронта лесного пожара, м/с.
При наличии ветра вместо высоты факела пламени необходимо использовать понятие длины факела пламени, а формулу (9) переписать в виде:
=^л/т
(10)
где Ь/ — длина факела пламени, м.
В работе [4] предложена простая зависимость для определения угла наклона факела пламени к подстилающей поверхности:
0,25
^у = [еД(4и2 )]
(11)
где у — угол между факелом пламени и горизонтальной плоскостью, рад; е — ускорение свободного падения, м/с2; Д—ширина фронта пожара (зоны разложения), м; ие — невозмущенная скорость ветра на высоте факела, м/с.
Ширину факела пламени согласно [4] можно вычислить по формуле
Д7 = Н
(л 2 ^0,2
4ие
еН
(12)
Длина факела зависит от запаса ЛГМ, поэтому его также необходимо знать. Если известна парциальная (насыпная) плотность слоя ЛГМ, то запас можно найти по формуле
тз = £ 1НР1 >
(13)
где £ — единичная площадь поверхности ЛГМ, м р° — парциальная плотность слоя ЛГМ, кг/м3;
Р0 = Р1 ф1 .
Зная угол наклона факела пламени к плоскости слоя ЛГМ и длину факела и используя зависимость
для углового коэффициента X/ из работы [1]:
( \
х /(х)=2
1+
Ь/ СОБ у - х
х 2 - 2 хЬ/ СОБ у + Ь2
(14)
можно определить значение лучистого теплового потока от факела в точке с координатой х по формуле
чЯ (х ) = (1 -ф5 ) Б/ О Т/ х/ (х ).
(15)
Здесь х —расстояние от основания факела пламени до нагреваемой точки, м; х = х - х/ ; Б/ — излу-чательная способность факела, которая согласно работе [5] определяется по формуле
б/ = б 0[1 - ехр (- А/ д/)], (16)
где б0 и А/ — эмпирические константы, значения которых для ряда углеводородных топлив приведены в работе [6].
Физический смысл углового коэффициента—это доля светового потока, излучаемого факелом, которая попадает в точку на поверхности с координатой х. Численно угловой коэффициент равен телесному углу, под которым виден факел пламени из точки х.
Скорость свободно-конвективного потока в факеле можно определить по эмпирической зависимости [7]:
V = 14,/ехр
11,5 (х - хт )
2
(17)
где ре — плотность газовой фазы в невозмущенной среде, кг/м3;
2
г
2 — вертикальная координата, в которой определяется скорость свободно-конвективного потока; хт — координата максимума температуры, м. Согласно [7] зависимость (17) соответствует установившемуся режиму течения в плоском двумерном факеле, возникающем от линейного источника тепла при отсутствии ветра.
Анализ численных расчетов
Краевая задача (1)-(8) с учетом экспериментальных зависимостей (10)-(17) решалась численно с применением итерационно-интерполяционного метода (ИИМ). Разностная схема ИИМ приведена в [8]. Теплофизические, кинетические и другие параметры, используемые в расчетах, были взяты из работ [2,6,7,9-12]: е„, = 0,9; р0 = 75 кг/м3; р1 = 1250 кг/м3; Ср1 = 1400 Дж/(кг-К); X1 = 0,102 Вт/(м-К); р2 = = 1000 кг/м3; Ср2 = 4190 Дж/(кг-К); X2 = 0,610 Вт/(мК); к0 = 0,4; р3 = 1360 кг/м3; Ср3 = 898 Дж/(кг-К); X3 = = 0,186 Вт/(м-К); р4 = 1300 кг/м3; Ср4 = 790 Дж/(кг-К); X4 = 0,137 Вт/(м-К); Ер = 35,97 кДж/моль; кр =1,54с-1; Чр = 0; к2з = 380 м/с; Е1з = 34,09 кДж/моль; ч2з = = 2,18106 Дж/кг; Е3, = 83,8 кДж/моль; к3, = 547 м/с
Таблица 1. Сравнение расчетных и экспериментальных данных
ч3
= 3,29 107 Дж/кг; кн = 8,5-10-7 м2; 5н = 1000 м-1; Л/=
= 0,018-0,026; е0 = 0,2-0,25; Т* = 573 К; Тн = 283 К; Тг = 1100 К; яС = 0,373; я, = = 0,382; я^ = 0,245 (где 5н — начальное значение удельной площади поверхности пор; Т* — характерная температура).
В процессе решения задачи необходимо рассчитывать лучистый теплообмен между слоем ЛГМ и факелом пламени, который зависит от размеров факела, невозмущенной скорости ветра, скорости распространения фронта пожара и запаса ЛГМ. Величины и Ду являются функциями процесса, поэтому при расчете длины факела скорость горения находили численно, а толщину факела Ду определяли по ширине зоны горения слоя, согласовывая расчетную и экспериментальную [13] скорости горения. Следует заметить, что численные расчеты хорошо согласуются с экспериментальными результатами, если ширину факела определять из условия Ду = хк -хн, где координаты хк и хн находились из условий:
Ф5(£> хн) = °>°5ф5н; Ф5(£> хк) = 0,95Ф5к> (18)
а параметры Лу и е0 принимать равными соответственно 0,018 и 0,25.
Скорость распространения очага горения определялась из уравнения
= ё х»/ё £
(19)
где х* — координата, удовлетворяющая условию Ф5(£, х*) = ВФ5н; 0< В < 1. Поскольку координаты границ зоны разложения вычисляются, для нахождения скорости горения по формуле (19) можно использовать одно из усло-
Параметр Значение параметра Значение параметра при Ж =12%
Расчет Ж = 1,5 % Эксперимент [9] Ж = 1-2 % Расчет Эксперимент [9]
Ь, м 0,175 0,180 0,114 0,110
и^, м/с 0,00273 0,00286 0,00129 0,00130
У, град 84,9 82,2 87,1 84,3
Примечание. Влагосодержание слоя ЛГМ вычислялось по формуле Ж = р2Ф2н ■ 100(р1 Ф1н).
вий (18). В рассматриваемых ниже расчетах использовалось первое условие (18). Константы 0,05 и 0,95 были подобраны из условия наилучшего согласования результатов численных расчетов и экспериментальных исследований. В частности, в табл. 1 приведены расчетные и экспериментальные данные основных параметров, определяющих процесс, при влагосодержании Ж, равном 1,5 и 12 %, и запасе опада хвои сосны тз = 0,3 кг. Из таблицы видно, что отличие расчетных данных от экспериментальных значений не превышает 5 %. В работе [7] константа Лу изменяется от 0,018 до 0,026, а е0— от 0,25 до 0,20. Если параметры Лу и е0 взять равными соответственно 0,026 и 0,20, то найденные численно длина факела пламени и скорость горения отличаются от экспериментальных значений на 16 %.
После отработки методики расчета радиационного теплообмена факела пламени со слоем ЛГМ (в отличие от работы [14], в которой исследуется механизм распространения пожара, влияние влаго-содержания и условий теплообмена с приземным слоем воздуха на скорость распространения пламени по слою опада хвои и газовый состав, выбрасываемый в приземный слой атмосферы) были проведены параметрические исследования влияния определяющих параметров на факел пламени и ширину противопожарного разрыва, который может преодолеть фронт горения.
В работе [7] приведена экспериментальная зависимость, которая определяет условия теплообмена слоя ЛГМ с приземным слоем воздуха при низовом пожаре при различной скорости ветра на высоте факела пламени. Согласно этой формуле при изменении скорости ветра на уровне высоты пламени от 0,01 до 9 м/с коэффициент теплоотдачи слоя ЛГМ с приземным слоем воздуха изменяется от 0,22 до 21 Вт/(м2 К). Поэтому в численных исследованиях было принято, что коэффициент теплоотдачи изменяется в указанном выше диапазоне.
На рис. 2 представлены зависимости длины факела пламени от условий тепломассообмена ЛГМ (при расчетах предполагалось подобие процессов тепло-
Рис. 2. Зависимость длины факела пламени при Ш = 3,5 %от условий теплообмена при запасе хвои сосны 0,75 кг (1)и3кг(2)
и массообмена) с приземным слоем воздуха при запасе массы 0,75 и 3 кг.
Проанализировав рис. 2, приходим к выводу, что при 0,22 < а < 21 и запасе хвои сосны 0,75 кг (соответствующем толщине слоя 0,01 м) длина факела пламени составит не более 0,72 м, а при запасе 3 кг (толщина слоя 0,04 м) может превысить 1 м.
Зависимости длины факела пламени и относительной ширины противопожарного разрыва от запаса массы хвои сосны приведены на рис. 3, из которого следует, что с увеличением запаса массы длина факела растет, а относительная ширина противопожарного разрыва Ьр/ Ь^ уменьшается.
Анализ зависимостей 1 и 2 на рис. 3 показывает, что при малом запасе ЛГМ длина факела увеличивается значительно быстрее ширины противопожарного разрыва по сравнению с большим запасом лесного горючего материала (кривая 1 на рис. 3). Сростом запаса эта зависимость выходит на асимптоту.
Данные исследования зависимости ширины противопожарного разрыва от влагосодержания ЛГМ представлены на рис. 4. Из рисунка видно, что при увеличении влагосодержания ЛГМ уменьшается вместе с длиной факела пламени и ширина противопожарного разрыва, так как отношение ширины противопожарного разрыва к длине факела (кривая 2 на рис. 4) остается постоянной величиной.
Изменения ширины противопожарного разрыва от условий конвективного теплообмена слоя ЛГМ с приземным слоем воздуха представлены на рис. 5. Из анализа данного рисунка следует, что при сравнительно малом конвективном тепломассообмене с возрастанием а ширина противопожарного разрыва растет быстрее длины факела пламени, а при а >3 Вт/(м2 К) длина факела пламени по мере роста а увеличивается быстрее ширины противопожарного разрыва. Из рис. 5 также следует, что отношение ширины противопожарного разрыва к длине факела пламени уменьшается с ростом запаса горючего материала.
Рис. 3. Зависимость относительной ширины противопожарного разрыва (1) и длины факела пламени (2) от запаса массы хвои сосны при Ш = 3,5 %, а = 2,21 Вт/(м -К)
ьр!Ч> Ч>м
Рис. 4. Зависимость длины факела пламени (1) и относительной ширины противопожарного разрыва (2) от влагосодержания хвои сосны при а = 2,21 Вт/(м -К), тз = 1,5 кг
Рис. 5. Зависимость относительной ширины противопожарного разрыва опада хвои сосны от условий конвективного теплообмена с приземным слоем воздуха при Ш = 3,5 %: 1 — тз = 0,75 кг; 2 — тз = 3 кг
На рис. 6 приведены зависимости угла наклона факела пламени к поверхности земли и скорости распространения низового пожара от скорости ветра.
Из рис. 6 следует, что чем больше скорость ветра, тем сильнее пламя отклоняется от перпендикуляра к поверхности слоя, а скорость перемещения очага горения возрастает.
Следует также отметить, что с увеличением запаса массы горючего вещества, несмотря на возрастание длины факела пламени (см. рис. 3), скорость распространения низового пожара по слою хвои сни-
Рис. 6. Зависимость угла наклона факела пламени к поверхности слоя опада ЛГМ (а) и скорости распространения низового пожара по слою хвои сосны (б) от скорости ветра на высоте пламени при Ш = 3,5 %: 1 — тз = 0,75 кг; 2 — тз = 3 кг
жается (см. рис. 6, б). Такой результат объясняется тем, что скорость свободно-конвективного течения в факеле с увеличением длины факела возрастает, и это приводит к повышению скорости подсоса в зону горения, т. е. образуется более сильное фильтрационное течение холодного воздуха в направлении, противоположном направлению распространения пожара.
Выводы
1. Предложенная математическая модель, в которой влияние приземного слоя воздуха на процессы горения рассчитывается по известным эмпири-
ческим зависимостям и проведено усреднение по высоте слоя ЛГМ, описывает процессы распространения низового пожара, согласующиеся с натурными экспериментами, например с работами [12, 13]. Данных по ширине противопожарного разрыва обнаружить не удалось.
2. Полученные зависимости можно использовать для оценки скорости распространения низового пожара и ширины противопожарных полос, размер которых зависит от запаса лесного горючего материала, влагосодержания и скорости ветра на высоте факела пламени.
СПИСОК ЛИТЕРАТУРЫ
1. Зигель P., Хауэлл Дж. Теплообмен излучением. — М. : Мир, 1975. — 934 с.
2. Гришин А. М. Математическое моделирование лесных пожаров и новые способы борьбы с ними.
— Новосибирск : Наука, 1992. — 407 с.
3. Телицын Г. П. Лесные пожары, их предупреждение и тушение в Хабаровском крае. — Хабаровск : ДальНИИЛХ, 1988. — 95 с.
4. Albini J. A. Physical model for fire spread in brush // 11th Symposium (Int.) Combustion. Berkley, Calif.
— Pittsburg, 1967. — P. 553-560.
5. Блох А. Г. Теплообмен в топках паровых котлов. — Л. : Энергоатомиздат, 1984. — 240 с.
6. Блох А. Г. и др. Теплообмен излучением. — М. : Энергоатомиздат, 1991. — 431 с.
7. Джалурия Й. Естественная конвекция. Тепло- и массообмен. — М. : Мир, 1983. — 399 с.
8. Гришин А. М., Зинченко В. И., Ефимов К. Н., Субботин А. Н., Якимов А. С. Итерационно-интерполяционный метод и его приложения. — Томск : Изд-во ТГУ, 2004. — 320 с.
9. Гришин А. М., Синицын С. П., Акимова И. В. Сравнительный анализ термохимических постоянных сушки и пиролиза лесных горючих материалов // ФГВ. —1991. — Т. 27, № 6. — С. 17-23.
10. Теплотехнический справочник / Под ред. В. Н. Юренева, П. Д. Лебедева. — М. : Энергия, 1976.
— Т. 2.— 896 с.
11. Кошкин В. Б., Синицын С. П. Определение эффективных кинетических постоянных пиролиза хвои сосны // Механика реагирующих сред и ее приложения. — Новосибирск : Наука, 1989. — С. 47-57.
12. Гришин А. М., Зятнин В. И., Перминов В. А. Экспериментальное исследование перехода низового лесного пожара в верховой. —Томск : ТГУ, 1991. —22 с. Деп. ВИНИТИ№ 982-13-91 от 06.03.91 г.
13. Конев Э. В. Физические основы горения растительных материалов. — Новосибирск : Наука, 1977.—239 с.
14. Субботин А. Н. Математическая модель распространения низового лесного пожара // Пожарная безопасность. —2008. — № 1. — С. 109-116.
Материал поступил в редакцию 15 июня 2012 г. Электронный адрес автора: [email protected].