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

МОДЕЛИРОВАНИЕ ПОЖАРОВ

' I

А. Ю. Катаева

д-р физ.-мат. наук, профессор Нижегородского государственного технического университета, г. Нижний Новгород, Россия

И. Е. Белоцерковская

ассистент Нижегородского государственного технического университета, г. Нижний Новгород, Россия

Д. А. Масленников

аспирант Нижегородского государственного технического университета, г. Нижний Новгород, Россия

А. А. Куркин

д-р физ.-мат. наук, профессор Нижегородского государственного технического университета, г. Нижний Новгород, Россия

УДК 519.6

СРАВНЕНИЕ АНАЛИТИЧЕСКОГО И ЧИСЛЕННОГО РЕШЕНИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ НИЗОВОГО ПОЖАРА С УЧЕТОМ ВЛИЯНИЯ УГЛА НАКЛОНА ПОДСТИЛАЮЩЕЙ ПОВЕРХНОСТИ

Целью данной работы является получение численного и аналитического решения задачи о распространении низового пожара с учетом влияния угла наклона подстилающей поверхности. Результаты, полученные для ровной поверхности, согласуются с известными данными, кактеоретическими [1, 5, 6], так и экспериментальными [3]. Ключевые слова: математическое моделирование; лесной пожар; подстилающая поверхность; численное и аналитическое решения.

Введение

Существует две основные модели, учитывающие влияние наклона в полуэмпирической модели распространения лесных пожаров: Доррера [5] и Ротермела [6]. В модели Доррера [5] учитываются уравнения теплового баланса в твердой фазе горючего в каждой точке слоя горючих материалов. Тепловое воздействие локального пламени на близлежащие слои горючих вводится на основе экспериментальных данных. При описании взаимодействия пожара с атмосферой он рассматривается как точечный очаг. Модель Ротермела [6] используется для оценки скорости распространения пожара, а также для вычисления его интенсивности, которая может быть интерпретирована как интенсивность в очаге пожара, и ширины фронта огня.

Постановка задачи

Физически задача ставится следующим образом: рассматривается движение плоского фронта пожара по склону (рис. 1); из решения гидродинамической задачи по обтеканию препятствия находится поле скоростей; известны температура окружающей среды, геометрические, структурные и

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

Для простоты анализа делаются следующие допущения:

1) среда является серой и пятифазной, включающей в себя сухое органическое вещество, воду в жидко-капельном состоянии, коксик, золу и газовую фазу;

2) газовая фаза состоит из кислорода, горючих компонентов продуктов пиролиза, инертных компонентов воздуха, а также водяного пара и инертных продуктов горения;

3) градиент температуры поперек полога леса мал по сравнению с градиентом температуры в продольном направлении.

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

© Катаева Л. Ю., Белоцерковская И. Е., Масленников Д. А., Куркин А. А., 2010

Рис. 1. Косоугольная система координат

Начальные и граничные условия для слоя горючего материала: пусть Ог — область фронта пожара; В (х, у, 2) — любая точка зоны пожара (рассматриваемого лесного массива); х0, у0 — центральная точка фронта пожара по осям Ох и Оу соответственно.

Тогда имеем: • при t = 0 и В е Ог

: ехр

х ехр

Т = (Ге - Те ) х

( 2 ^ 2 ^ (х - хо )2 - (У - Уо )2

V ¿Г

Гу

Т, = (Те - Те ) X

( 2 2 ^ --^ (х - хо )2 - (У - Уо )2

V ¿Г

Гу

+ Те; (1)

+ Те; (2)

(

са = (cаg сае ) ехр

2 ( А

(у - уо)

4у

2 ( А -~^(х -хо) -V А2гх

+ с а = 12'

ик = ике; Ф1 = Фlg; ф2 = Ф2g;

Ф з = Ф Зg; Ф4 = Ф4g;

и = • к = • ж = М ^ )Ш И2 ; К И2 ; ^ 1п(А)

• при t = 0 и В £ ОГ

Т = Т • Т = Т • с = с

ик = иКе ; Ф1 = Ф1е ; ф2 = Ф2е ;

Ф3 = ФЗе ; Ф4 = Ф4е ;

1П( 2 )Ше

и =

2 2и 2 2У

е ; V = е ; Ш =_

И2 ; И2 ; 1п(И)

(3)

(4)

(5)

(6)

(7)

(8)

Уравнение сохранения объемной доли компонентов газовой фазы имеет вид:

4

Ф 5 = 1 -ЕФ''.

1=1

Список обозначений, использованных в формулах, приведен в конце статьи.

Учет рельефа местности при численной реализации

Для учета рельефа местности в данной постановке был осуществлен переход в косоугольную систему координат. Пусть будет Ох' ортонормирован-

, к

ная система координат, а Лг — косоугольная система координат с началом в точке А, координаты

которой обозначеныX' (см. рис. 1). Пусть для изме-

и . к

рения длин вдоль осей Л2 установлены некоторые масштабы Вк. Возьмем какую-нибудь точку Б, координаты которой в системе Ох1 обозначены через х. Разложим вектор г = ЛБ на составляющие по осям х'; концы этих составляющих обозначим через Лгк, а длину отрезков ЛБк, измеренную единицей длины, обозначим через . Результат построения длин тех же отрезков при помощи выбранных масштабов Вк обозначим через гк и назовем косоугольными координатами точки Б.

Таким образом, мы имеем 2а = Ва2а. Если мы припишем величинам 2а и Ва размерность длины, то косоугольные координаты 2а будут безразмерными величинами. В ортонормированной системе не имеет смысла различать координаты и составляющие, поэтому условимся, что координаты точек в системе Ох1 имеют размерность длины. Тогда и составляющие вектора г в системе Ох1, равныех' -X', будут иметь размерность длины. Обозначим направляющие косинусы осей Лгк относительно Ох' через а к. Тогда можно записать:

х' - X'' = а

2к =

= а1 В1 г1

а

' 73 2 2 В 2 2

а

3 В 3 2

Введем объект

Рк =

В1 а 1

В1 а^ В2а2 В3а3

В2 а2

В 3 а 3

В1 а3

В2 а3

В 3 а 3

(9)

(10)

Если ввести в рассмотрение базисные векторы Вк, как это часто делают, то '-я строка объекта Р к представляет собой проекции базисных векторов на ось х .

Итак, при помощи объекта Рк можно записать:

' V о' к х - Х =Р к2

(11)

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

что декартовы координаты х' выражены в виде функ-

к

ций от трех параметров д :

х' = /' (д1

д2 д3).

(12)

Условимся, что /' однозначны, непрерывны и имеют частные производные всех порядков, которые

нам понадобятся. Будем считать, что уравнения (12)

к

разрешены относительно д , отсюда вытекают равенства:

дк =фк (х1, х2, х3). (13)

Причем функции фк однозначны, непрерывны и допускают все нужные нам частные производные. При этих условиях каждая точка определяется как тремя числами х', так и тремя числами дк. Последние называются криволинейными координатами.

Чтобы пояснить смысл этого термина, фиксируем некоторую точку А (рис. 2). Пусть ее декартовы координаты будут X', а криволинейные — Qk, причем Qk = фк (X1, X 2 X3). Если в уравнениях (12) зафиксировать какие-нибудь два параметра, то координаты х станут функциями только одного параметра, вследствие чего уравнения (12) определят некоторую кривую. Очевидно, что через точку А проходят три таких кривых:

кривая д :

кривая д :

кривая д :

х'1 = Г (д1, Q2 Q3);

х''2 = Г (Q1, д2 Q3); х''3 = Г(Q1, Q2 д3).

(14)

В криволинейных координатах уравнения этих кривых будут соответственно иметь вид:

кривая д :

кривая д :

кривая д :

2 г\2

д = Q :

1 /->1

д = Q >

1

д = Q >

3 г,3

д = Q ; д = Q ;

2 ^>2 д = Q .

(15)

Эти кривые называются координатными. Поскольку через каждую точку пространства проходят три координатные кривые (см. рис. 2), можно сказать, что точка определяется как пересечение трех координатных кривых.

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

Рис. 2. Криволинейная система координат

Возьмем точку В (см. рис. 2), бесконечно близкую к А, и обозначим расстояние АВ через йт. Соответствующие дифференциалы декартовых координат обозначим через йхБ:

йх5 =

ддк

йдк = $кйдк.

(16)

У А

Здесь индекс "А" у производных означает, что после дифференцирования нужно подставить дк = Qk, а

(дГ V ддк) а =Р к.

Запишем частные производные по декартовым координатам через криволинейные координаты (воспользуемся оператором свертки для более короткой записи):

д_ дх

д дд1 дд1 дх

д дд 2 д дд 3

д дд к

3

=

дд2 дх дд3 дх к = 1 ддк дх

д = — д дд ду к = 1 ддк ду

д = — д дд дх к = 1 ддк дх

(17)

Учитывая знания о том, как сформирована система криволинейных координат (а именно: для любой точки А, принадлежащей рассматриваемой области, ось Ог всегда совпадает по направлению с осью Ад3; ось Ох всегда перпендикулярна Ад2, так как последняя лежит в перпендикулярной к Ох плоскости, параллельной Оуг; ось Оу всегда перпендикулярна Ад1, так как последняя лежит в перпендикулярной к Оу плоскости, параллельной Охг), получим:

дх2

д(д1 )2

дд1

дх

д д 2 д1

1 ;

дд дх

2

2

2

dy

d( q 2)2

С Л 2\2

dq dy

э д2 2

d d q dq2 dy 2

3 f 3

^r -S(S-i=Aj -

dq dqj

dz 2 k - 1I j -1 dqk dq j dz

dz

я д2 k ^

d d q dqk dz2

3

-s

k -1

d(qk )2

k2

dq

dz

d d2qk dqk dz2

dq dqj

2 3 d k-i j -k +1 dq dq j dz dz

(18)

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

Методика получения аналитического решения

Используя систему обыкновенных дифференциальных уравнений [1], повернем ось Ox на угол Р:

cos Р - x ^ x - x cos Р;

х

d d 1 d

dx d(x cos Р) cos Р dx' и введем следующие обозначения:

u® - ■

cos Р

. _ DT ; _ 3 т DT 2 а ; т 2

cos2 Р

cos2 Р

(19)

Переходя к декартовой координате в подвижной системе отсчета, связанной с положением максимального профиля температуры ^ = х' - ю^ предположим, что весь фронт движется с одной и той же скоростью ю. Получим систему уравнений:

3ф1 -

Р1 = к1' К1 = к01Р1Ф1е ш ;

Р 2 »^ - R

_К2

dm 3

Р3лт - - аcR1

2> R2 - k02Р 2 m2e M„

(20)

d %

R3 - k03 ^аР 5 cl93 e

M

R

E3

R3. dLL - 0;

-Rf. dh d %

d

d %

dca

d%

dca

— |p 5 D'T I -Р 5 (U® + R 5a - Ca Q =

d%

a

- - (Ca - Ca® A a = 3 2,3;

Cpa

T^T Ц)- [Р 5(U»-®)CP5 -

-®SPimiCpi]dT%-a (т - T®) /-1

(21)

p!Jd % h + q3R3 - q2R2 + q5R5 - 0;

(22)

M1

R51 - -R3 R5; R52 - (1 -aC)vrR1 -R5;

2M

R53 - 0; R5 - M2k05

Р 5 -

P®

c1 M v,2j f c2M

M

0,25

M

T-2,25 eRT ;

RT® S -Ma

z -1 a

; S Ca -1;

i-1

Q - (1 -a c R + R2 + M^R3.

Граничные условия к системе (20) - (22):

% - ±®:

T Tl ; ca C

a®' T(

m z - m i

dc a

d%

- 0; dT ; d%

- 0;

. dm z

(23) - 0;

■ да 3 ^

Тда < Т < Т*:

Л, = 0; Я5а = 0; Я5 = 0; а =1, 2; ' = 1, 2, 3. (24)

Как частный случай рассматривалось распространение плоской волны горения по наклонной поверхности. Такое предположение позволило детально изучить влияние наклона подстилающей поверхности на скорость распространения пожара. Распределение температуры газовой фазы при различных моментах времени и различных углах наклона подстилающей поверхности представлено на рис. 3. При х = 0 задавалось распределение температуры в экспоненциальной форме с максимумом 1200 К. Из рисунка отчетливо видно, что скорость распространения плоской волны горения растет с увеличением угла наклона подстилающей поверхности.

Проанализируем поведение решений уравнений для температуры Т и концентрации компонентов са при ^ ^ ± да. Представим решение в виде:

Т = Т + Тс = с + с' ;

-1 1 да ^ 1 > а сада а'

, (25)

Ф'' = Ф''да + Ф'' >

где штрихом обозначена малая величина возмущения, характеризуемая отклонением параметров от равновесных значений. Подставляя (25) в (21), (22)

2

2

2

2

E

u

f = 0 с

0,68

i = 50 с

Х 1500 2000 ^

Рис. 3. Динамика распределения температуры плоской волны горения вдоль оси х при углах наклона а, равных 0, 15, 30 и 45°, в моменты времени 0, 50 и 100 с

и отбрасывая малые величины более высокого порядка, получаем:

Р 5» D

/ d ca

d i2

,, d2T'

di

/ ' \ dca a ' rv

-Р 5» (u ---ca =

d i cpa

a =1,2, 3;

4

Р 5» ( u» -ю) cp 5 -®ZP i Ф icp

(26)

x ^-a t ' = o. d i h

(27)

Решением уравнений (26) и (27) при ^ ^ ± да будет:

Т'= «1,2е; са = 61,2а ^. (28)

Для диффузного режима горения с хорошей степенью точности можно полагать:

Т(0) = Т0, с1 = 0 = с2 = 0 = 0;

dc1

1=o

dc2

= 0; ф1(0) = 0; ф2(0) = 0. (29)

1=0

Пусть p 5 (и» - ю) = p 5» (и» - ю) = const приближенно выражает закон сохранения массы при лесном пожаре. При 0 < i < » профиль температуры определяется формулами (25) и (28), где a1 = Т0 -— постоянная, определяемая из условий (29).

Интегрируя уравнение (22) по 0 < i < », используя уравнения (20) и (21) для вычисления интегралов и обозначив

q =

Чз M1 a c

M„

+ Ч 5 (1 -a c И

c1 Ф1Н + Ч2Р 2Ф2H j (P 5cp5T» ) 1;

_ ю

ZP i Ф icpi

© 0 =7»; ю = —; л =

i=1

p 5» cp 5

. 0,38 -

й 0

5 0,63

6

0,505

0,349

0,698

0,38

0,349 Угол, град.

0,698

Рис. 4. Численное (у) и аналитическое (ю(Р)) решения зависимости скорости распространения пожара от угла наклона Р при ©„ = 0,5 (а)и ©„ = 0,9 (б):----у;--ю(Р)

получим:

Чю + [1 -ю (1 + л)](©0 - 1) = 2e(©0 -1)

(30)

[1 - ю(1 + л)]2 + 4e - 1 + ю(1 + л) Решая (12) относительно ю, получаем

ю(Р) = ■

Ч + q2 - 4 e (©0 - 1)1 q(1 + л)

(©0 - 1)

2 cos(P) I q(1 + л)

ч 2

(©0 -1)

(31)

На рис. 4 представлена зависимость скорости распространения пожара, полученная численно (у) и аналитически (ю(Р)) для значений 00 = 0,5 (а) и 00 = 0,9 (б), от угла наклона Р (рад). Из рисунка видно, что с ростом угла наклона скорость распространения пожара также возрастает. Здесь под численным решением подразумевается решение задачи в полной постановке. Как видно из анализа полученных результатов, упрощенное (аналитическое) решение дает оценку сверху для значений 00 < 0,5 и оценку снизу для значений 00 > 0,9.

На рис. 5 показано сравнение моделей Доррера [5], Ротермела [6] и Гришина [1]. При использовании эмпирических данных модели Ротермела в модели Доррера была определена зависимость скорости распространения от угла наклона и показано, что модель Ротермела при углах наклона до 23° дает завышенные значения скорости распространения пожара, а модель Доррера дает значения, близкие к модели Гришина, однако при углах наклона более 23° модель Доррера имеет ограничения, модель Рот-термела дает заниженные значения скорости распро-

2

и

0

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

странения пожара, в то время как даже упрощенная модель Гришина при этом не имеет ограничений.

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

Представленные результаты поисковой научно-исследовательской работы получены в рамках реализации мероприятия 1.2.1 "Проведение научных исследований научными группами под руководством докторов наук" ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009- 2013 гг.

Условные обозначения

t — время;

Т — температура газовой фазы; Тц — температура твердой фазы; Т0 — максимальная температура во фронте лесного пожара;

Тх — температура окружающей среды вдали от фронта пожара, принимается равной 300 °С; са — массовая концентрация компонентов газовой фазы (а = 1 — кислород; а = 2 — горючие компоненты продукта пиролиза; а = 3 -инертные компоненты газовой фазы, не реагирующие компоненты продукта пиролиза и водяного пара; с1 = 0,23; с2 = 0; с3 = 0,769);

и, V, Ж — проекции скорости над пологом леса по

осям соответственно х, у, 2; ф1 — объемная доля газовой фазы; ф1 = р/р1 (р = 0,3 - плотность слоя; р1 = 360 — плотность древесины);

ф2 — объемная доля горючих компонентов продуктов пиролиза; Ф2 =р(1 -ОЖигм/Р2 (Р2 = = 1000 — плотность древесины; £ = 0,3 — процентное содержание золы в лесных горючих материалах (ЛГМ); ЖЛГМ = 0,35 — влагосодержание ЛГМ); ф3 — объемная доля инертных компонентов воздуха; Ф3 = 0,01; ф4 — объемная доля водяного пара и инертных продуктов реакций окисления, пиролиза и горения кокса; ф4 = 0; И — высота слоя горючих материалов; Р — угол наклона;

х — декартова координата в системе координат, связанной с поворотом оси Ох на угол Р (ось х направлена в сторону невозмущенной скорости ветра, параллельной горизонтальной подстилающей поверхности); и — равновесная скорость ветра; и = 0,2; ида — компонента скорости в выбранной системе

координат, вычисляется согласно (19); р1 — плотность сухого органического вещества; р2 — плотность воды в жидкокапельном состоянии; р3 — плотность конденсированных продуктов пиролиза; р3 = 200; р4 — плотность минеральной части (золы); р4 = 200; р5 — плотность газовой фазы (вычисляется согласно [5]); Р — давление в потоке; р да — давление; рсв = 1,01105; Q — массовая скорость образования газовой фазы; Я1 — массовая скорость реакции пиролиза сухого

органического вещества ЛГМ; Л2 — массовая скорость реакции испарения воды из ЛГМ;

Л3 — массовая скорость реакции горения коксового остатка;

Л5 — массовая скорость реакции горения газообразных горючих продуктов пиролиза; Л5а — массовая скорость образования компонентов газовой фазы; Л51, Л52, Л53 — массовые скорости образования компонентов газовой фазы; М1 — молекулярная масса газовой компоненты кислорода; М1 = 32; М2 — молекулярная масса газовой компоненты горючих компонентов пиролиза; М2 = 28; М3 — молекулярная масса газовой компоненты

инертных компонентов воздуха; М3 = 28; Ыа — молярная масса индивидуальных компонентов;

Мс — молярная масса углерода; Мс = 12; М — молярная масса смеси газов в целом; Х'Т — коэффициент турбулентной теплопроводности; Х'т =ц т Ср 5; — коэффициент динамической вязкости; ^1>2, й1>2а — корни соответствующих характеристических уравнений; а12, Ь1 2а — произвольные постоянные. ср1 — теплоемкость фазы сухого органического

вещества; ср2 — теплоемкость фазы воды; ср3 — теплоемкость фазы конденсированных продуктов пиролиза угля; ср4 — теплоемкость фазы минеральной части золы; ср5 — теплоемкость газовой фазы; БТ — коэффициент турбулентной диффузии; ХТ — коэффициент турбулентной теплопроводности;

аТ — коэффициент внутреннего теплообмена; Я — универсальная газовая постоянная; д2 — тепловой эффект процесса испарения; q3 — тепловой эффект процесса горения конденсированного горючего; q5 — тепловой эффект процесса окисления газообразного горючего продукта пиролиза; к01 — предэкспонента реакции пиролиза; к02 — предэкспонента процесса испарения; к03 — предэкспонента процесса горения конденсированного продукта пиролиза (коксика); Е1 — энергия активации реакции пиролиза; Е2 — энергия активации процесса испарения;

Е3 — энергия активации процесса горения конденсированного продукта пиролиза (коксика); 5С — удельная площадь поверхности конденсированного продукта пиролиза; ас — коксовое число ЛГМ;

уг — доля горючего газа в газообразных продуктах пиролиза; иЯ — плотность потока излучения; к^ — коэффициент ослабления; к, — спектральный коэффициент поглощения; с — скорость света; а — постоянная Стефана-Больцмана; Л, — коэффициент межфазного теплообмена; у5 — доля теплоты газофазной реакции окисления газообразных продуктов пиролиза, усвоенная конденсированной фазой; у5 < 1; $ — удельная площадь поверхности фитомассы; сз — эмпирический коэффициент сопротивления; g — ускорение свободного падения.

Индексы

— невозмущенный поток; е — окружающая среда;

g — поверхность;

— параметры на верхней границе полога леса; + — параметры на нижней границе полога леса.

Характеристики очага горения

х0, у0 — координаты центра очага горения; А Гх, А Гу — ширина очага загорания по координатам х, у соответственно.

СПИСОК ЛИТЕРАТУРЫ

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

2. Катаева Л. Ю. Особенности дискретизации многомерных нелинейных задач // Наука и техни-катранспорта. — 2008. — №4. — С. 13-16.

3. Катаева Л. Ю., Романов А. В. Метод Патанкара и возможности его оптимизации // Наука и техника транспорта. — 2009. — №3.— С. 88-97.

4. Butler В. W., Anderson W. Я., Catchpole Е. A. Influence of Slope on Fire Spread Rate / The Fire Environment-Innovations. Management and Policy. — P. 75-85.

5. Доррер Г. А. Модель распространения фронта лесного пожара // Теплофизика лесных пожаров. — Новосибирск : ИТФ СО АН СССР, 1984. — С. 86-99.

6. Яо^егте1 Я. C. A Mathematical Model for Predicting Fire Spread in Wildland Fuels // Intermountain Forest & Range Experiment Station, Forest Service, US Dept. of Agriculture, 1972.

Материал поступил в редакцию 21 июня 2010 г.

Электронный адрес авторов: dmitrymaslennikov@rambler.ru.

3о| ISSN 0869-7493 ООЖАРООЗРЬЮОБЕЗООАСООСТЬ 2010 ТОМ 19 №11