РЕШЕНИЕ ИНТЕГРОДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ТЕПЛОПЕРЕНОСА В ЗАДАЧЕ МОДЕЛИРОВАНИЯ ПРОЦЕССА РАСПРОСТРАНЕНИЯ ЛЕСНОГО ПОЖАРА
ТЕПЛОВЫЕ РЕЖИМЫ ПРИБОРОВ И СИСТЕМ
УДК 519.8
С. А. АСТАФЬЕВ
РЕШЕНИЕ ИНТЕГРОДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ТЕПЛОПЕРЕНОСА
В ЗАДАЧЕ МОДЕЛИРОВАНИЯ ПРОЦЕССА РАСПРОСТРАНЕНИЯ ЛЕСНОГО ПОЖАРА
Изложен ход решения интегродифференциального уравнения теплопереноса, используемого как базовое для описания изменения параметров теплового потока в модели распространения лесного пожара. Ключевые слова: моделирование процесса распространения пожара, тепловой поток, интегродифференциальное уравнение теплопереноса, напряжение трения теплового потока.
Введение. На активно охраняемой территории лесного фонда России ежегодно регистрируется от 10 до 35 тыс. лесных пожаров, охватывающих площади 0,5—2,5 млн га. Уничтожение бесценных запасов древесины в результате лесных пожаров сопровождается выгоранием большого объема кислорода и выбросом в атмосферу сажи, копоти, двуокиси углерода, а также радиоактивных частиц [1].
Физика распространения лесного пожара состоит в следующем: тепловой поток, создаваемый низовым пожаром, поднимаясь наклонно по направлению ветра, подогревает кроны деревьев на значительном расстоянии впереди фронта огня. При воспламенении хотя бы одной из крон почти мгновенно воспламеняются и другие, и огонь „скачет“ по подогретым кронам, но затем вне сферы действия подогрева затухает. Когда низовой огонь приближается к фронту пожара, процесс подогрева полога повторяется и опять происходит „скачок огня“. Верховые пожары, выделяя большое количество теплоты, вызывают восходящие потоки продуктов горения и нагретого воздуха и образуют конвективные колонки диаметром в несколько сотен метров. Их поступательное движение совпадает с направлением продвижения фронта пожара. Пламя в середине колонки может подниматься на высоту до 100—120 м. Конвективная колонка увеличивает приток воздуха в зону пожара и порождает ветер, который усиливает горение. Таким образом, низовой пожар стимулирует развитие верхового и наоборот [1].
Основными задачами охраны лесов от пожаров являются их предупреждение, обнаружение, ограничение распространения и тушение. Для эффективной борьбы с лесными пожарами необходим прогноз возможного положения кромки пожара и силы горения.
При математическом моделировании процесса распространения лесного пожара можно выделить четыре группы моделей:
1) модели прогноза скорости распространения пожара; 2) модели прогноза контура и площади пожара;
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
Решение интегродифференциального уравнения теплопереноса
51
3) модели прогноза характеристик теплового течения и массопереноса во фронте и в зоне пожара;
4) общие математические модели, в рамках которых могут быть предсказаны все характеристики во фронте и в зоне пожара (в том числе, поля температур, концентрации и скорости распространения газовых и мелкодисперсных компонентов окружающей среды) [2].
Однако большинство разрабатываемых в настоящее время моделей базируются либо на статистических данных о скорости распространения огня в зависимости от погодных условий, либо на рассмотрении лесного пожара с позиций аэротермохимии. Исходя из этого более перспективным в решении задачи моделирования процесса распространения лесного пожара представляется использование одновременно двух подходов.
При первом — вероятностном — подходе [3, 4] модель прогноза контура и площади пожара к заданному моменту времени строится на основе определения вероятностей распространения огня через дискретные области (по результатам обработки аэрофотоснимков местности) и среднестатистических скоростей распространения пожара в разных направлениях в зависимости от скорости и направления ветра.
При втором подходе — с использованием интегродифференциального уравнения теплопереноса — для расчета параметров теплового потока в зависимости от местоположения относительно фронта пожара используются физические характеристики окружающей среды; кроме того, такой подход позволяет учитывать данные о рельефе местности.
В настоящей статье предложено решение интегродифференциального уравнения теплопереноса для применения его при моделировании процесса распространения лесного пожара.
Интегродифференциальное уравнение. Для моделирования за основу было взято интегродифференциальное уравнение, приведенное в работе [5]. Автор данного уравнения (профессор Томского государственного университета А. М. Гришин) предложил его в качестве математической модели стационарных двумерных квазиравновесных тепловых течений, описывающих процесс распространения лесного пожара:
( )u
∂u ∂x
+
⎛ ⎜ ⎝
vw
−
y
∫
0
∂u ∂x
dy
⎞ ⎟ ⎠
∂u ∂y
=
∂ ∂y
µρ
∂u ∂y
−
sin α Fr
−
1 ρ
∂P ∂x
−
∂P ∂y
y
∫
0
1 ρ
∂
ln P ∂x
dy
,
где и — скорость теплового потока; vw — безразмерная скорость внешнего потока (вдува); µ=18,27·10–6 Па·с — динамическая вязкость воздуха; ρ=1,2 кг/м3 — плотность воздуха;
Fr — критерий Фруда (безразмерная величина) — один из критериев подобия движения жид-
костей и газов; α — угол между касательной к поверхности и горизонтальной плоскостью;
P — давление окружающей среды в зоне пожара.
Система координат xoy, используемая в уравнении, показана на рис. 1.
y
Поле ветра
ymax
Перенос энергии
Приземный слой атмосферы
Фронт пожара
Лес
x
0
Рис. 1
Метод решения. В 80-е годы прошлого века, когда прикладные компьютерные программы для ведения сложных математических расчетов еще не обладали достаточной производительностью и гибкостью и не были широко распространены, автором уравнения было осуществлено его аналитическое решение с использованием различных методов, допущений
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
52 С. А. Астафьев
и замен переменных (например, вывод уравнений методом Прандтля, приближение Буссинеска, переход к переменным Дородницына — Хоуарта, решение уравнений методом Швеца). Как следствие, аналитическое решение является достаточно громоздким, а большое количество преобразований повышает вероятность ошибки; кроме того, к недостаткам такого способа решения можно отнести большие затраты времени на пересчет результатов при изменении входных данных и отсутствие быстрой визуализации результатов вычислений.
На современном этапе степень развития вычислительных систем и приложений позволяет произвести численное решение сложного интегродифференциального уравнения с использованием среды математического моделирования MatLab. Данное программное средство содержит большой набор встроенных функций, предназначенных для решения обыкновенных дифференциальных уравнений. Так, функции ODE реализуют одношаговый явный метод Рунге — Кутты 2-го, 3-го порядка (ode23) и 4-го, 5-го порядков (ode45); неявный метод Рунге — Кутты, использующий формулы обратного дифференцирования 2-го порядка (ode23tb); адаптивный многошаговый метод Адамса — Башворта — Мултона переменного порядка (ode113); одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка (ode23s), и метод трапеций с интерполяцией (ode23t). Функция pdepe служит для решения систем параболических и эллиптических дифференциальных уравнений в частных производных по двум переменным при заданных начальных и граничных условиях. Данная функция преобразует уравнения в частных производных в набор обыкновенных дифференциальных уравнений с использованием дискретизации 2-го порядка точности на основе фиксированного набора заданных пользователем узлов. Интегрирование осуществляется с помощью функции ode15s, которая реализует адаптивный многошаговый метод переменного порядка (от 1 до 5), использующий формулы численного дифференцирования.
Таким образом, главной задачей при решении исходного уравнения в среде MatLab является его запись в форме, адекватной предусмотренной в программном пакете.
Исходные данные, начальные и граничные условия. В качестве параметра vw и координат x, y в процессе решения уравнения использовались безразмерные величины (исходя из предположения, что vwmax=30 м/с — максимальная скорость внешнего потока; ymax=30 м — характерная высота пограничного теплового слоя; xmax=120 м — характерное расстояние теплового воздействия фронта пожара).
В качестве исходных данных условно было задано давление среды в зависимости от координат, выраженное формулой
{ }P(x, y) = 0, 2 + exp −0, 48(x − 0,5)2 − 0,16( y − 0,8)2 .
Р
1,15 1,1 1,05 х 0,4 0,3 0,2 0,1
0
Такой вид зависимости был выбран
1,18 исходя из предположения, что в зоне про-
1,16 грева давление среды экспоненциально рас-
1,14 тет с увеличением высоты (у) над поверхно-
1,12 1,1
стью и расстояния (х) от задней границы фронта пожара, и после достижения максимального значения давление вновь снижает-
1,08 ся до нормального. Заданное изменение
1,06 давления в зависимости от координат гра-
1,04 фически представлено на рис. 2.
0,2 0,4 0,6 0,8 y Рис. 2
1,02 В качестве начального условия (при x=0) было принято экспоненциальное увеличение скорости теплового потока в зависимости от высоты:
u0 ( y) = 1, 0 −1, 0 exp (−2,5 y) .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
Решение интегродифференциального уравнения теплопереноса
53
В качестве нижнего граничного условия (при y=0), согласно данным из работы [5], было
принято нулевое значение скорости теплового потока: uн = 0 , а в качестве верхнего гранично-
го условия (при y=ymax) — экспоненциальное уменьшение скорости теплового потока с увеличением расстояния от фронта пожара:
uв (x) = 1, 0 exp (−17x) .
Файл-функция, задающая граничные условия, содержит два параметра, определяющие
отсутствие „потокопроводности“ на нижней границе пограничного слоя и свободное распро-
странение теплового потока на его верхней границе.
Ход решения. На первом этапе решения исходное интегродифференциальное уравне-
ние теплопереноса было упрощено путем исключения из него интегрального слагаемого
∫ymax
∂u ∂x
dy
.
0
Решение параболического дифференциального уравнения в частных производных про-
изводилось с использованием функции pdepe, имеющей следующий вид:
sol=pdepe(m, @teplofun, @teploinit, @teplobound, y, x),
где u=sol(:,:,1) — искомая скорость теплового потока; m — параметр, соответствующий типу дифференциального уравнения; teplofun — функция, определяющая компоненты решаемого дифференциального
уравнения; teploinit — функция, определяющая начальные условия; teplobound — функция, определяющая граничные условия. Для решения дифференциального уравнения в пакете MatLab оно должно иметь сле-
дующую форму записи:
( ) ( ) ( )c
y,
x,
u,
∂u ∂y
∂u ∂x
=
y−m
∂ ∂y
⎛⎝⎜
ym
f
y,
x,
u,
∂u ∂y
⎞⎟⎠ + s
y,
x,
u,
∂u ∂y
.
Для параболического дифференциального уравнения m=0. Таким образом, исходное уравнение должно быть записано в виде следующей системы:
c = u;
⎫
f
=
⎛ ⎜ ⎝
−vw
+
y
∫
0
∂u ∂x
dy
⎞ ⎟ ⎠
u
+
µρ
∂u ∂y
,
⎪ ⎪⎪ ⎬
s
=
−
sin α Fr
−
1 ρ
∂P ∂x
−
∂P ∂y
y
∫
0
1 ρ
∂
ln P ∂x
⎪ dy.⎪
⎭⎪
В
первую
очередь
были
рассчитаны
частная
производная
давления
по
координате
x
—
∂P ∂x
,
частная производная давления по координате y —
∂P ∂y
и
частная
производная
логарифма
дав-
∫ления по координате x —
∂ ln P ∂x
.
График
значений
интеграла
A
=
ymax
1 ρ
∂ ln P ∂x
dy
представлен
0
на рис. 3.
В результате решения дифференциального уравнения теплопереноса были рассчитаны
значения скорости теплового потока в зависимости от координат u(x,y), частная производная
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
54 С. А. Астафьев
скорости теплового потока по координате x —
∂u ∂x
и частная производная скорости теплового
потока по координате y —
∂u ∂y
(рис. 4, а—в).
А
0,35
0,3 0,25 0,2
0,15
0,1 0,05
0
0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х Рис. 3
а) и
1,2 1 0,8 0,6 0,4 0,2 х
0,3 0,2 0,1
0
0,2 0,4 0,6 0,8
y 1
б)
du dx 40 20 0
–20 х 0,3 0,2 0,1
0,40,60,81 y 0,2 0
в) du dy 2 1,5 1 0,5
х 0,3 0,2 0,1
0 0,2 0,4 0,6 0,8 1 y
Рис. 4
Затем по формулам
τw
=
µρ
∂u ∂y
y=0
и
τe
=
µρ
∂u ∂y
y= ymax
были вычислены напряжение тре-
ния теплового потока на обтекаемой поверхности и напряжение трения теплового потока на
внешней границе пограничного слоя (рис. 5).
Физический смысл параметра τ заключается в том, что при τw = 0 происходит „отрыв“
теплового потока от поверхности, что означает смену типа теплового течения: вместо течения, при котором линии потока вдали от фронта пожара почти параллельны обтекаемой поверхности, образуется конвективная колонка.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
Решение интегродифференциального уравнения теплопереноса
55
Вид зависимости τw(x) вполне согласуется с графиком, представленным в работе [5]:
с увеличением х значение τw уменьшается, причем его резкое уменьшение в пределах фронта
пожара далее переходит в плавное и стремится к нулю.
Вместе с тем следует заметить, что в операциях интегрирования и дифференцирования
(с использованием функций trapz и diff) результат вычислений зависит от шага интегри-
рования (дифференцирования), и поэтому необходимо использовать масштабные коэффи-
циенты.
Одна из особенностей решения состоит в том, что вычисления производились методом
последовательных приближений. Так, в первом цикле вычислений рассчитывалось значение
∫интеграла
B
=
ymax
∂u ∂x
dy
,
которое
затем
использовалось
в последующей
итерации. В резуль-
0
тате после нескольких расчетных циклов значение интеграла, вычисленное по окончании
программы, сходится к значению, заданному на входе. График значений интеграла В пред-
ставлен на рис. 6.
τw, τe, о.е.
0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х
у у=ymax
2 –1
1,5 –2 –3
1 –4
0,5 –5
0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х
В
Рис. 5
Рис. 6
Заключение. Представленное решение интегродифференциального уравнения теплопе-
реноса осуществлено при условно принятых начальных и граничных условиях, подобранных
таким образом, что обеспечивается устойчивое решение уравнения. Ход вычислений отдель-
ных параметров имеет графическую визуализацию. Впоследствии необходимо произвести
согласование параметров модели с реальными условиями распространения огня путем прове-
дения натурных экспериментов и исследований параметров реальных пожаров.
Таким образом, использование уравнения теплопереноса в основе модели распростра-
нения лесного пожара позволяет в перспективе обеспечить возможность прогнозирования
эволюции пожара путем учета конкретных физических параметров окружающей среды. Сре-
ди параметров модели наиболее значимым, с точки зрения предсказания развития лесного
пожара во времени, представляется напряжение трения теплового потока. Дальнейшие иссле-
дования целесообразно направить на выявление характерных изменений данного параметра в
зависимости от внешних условий.
СПИСОК ЛИТЕРАТУРЫ
1. Воробьев Ю. Л., Акимов В. А., Соколов Ю. И. Лесные пожары на территории России: Состояние и проблемы / Под общ. ред. Ю. Л. Воробьева. М.: ДЭКС-ПРЕСС, 2004. 312 с.
2. Гришин А. М. О математическом моделировании природных пожаров и катастроф // Вестн. Томск. гос. ун-та. Математика и механика. 2008. № 2(3).
3. Астафьев С. А., Лысенко Д. Ю., Широков А. С. Моделирование процесса распространения лесного пожара с применением теории перколяции // Изв. вузов. Приборостроение. 2012. Т. 55, № 6. С. 70—74.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
56 Ю. В. Баёва, Е. В. Лаповок, С. И. Ханков
4. Астафьев С. А. Применение вероятностного подхода в задаче моделирования распространения лесного пожара // IX Всерос. конф. молодых ученых; V сессия науч. школы „Проблемы механики и точности в приборостроении“: Сб. докл. СПб: НИУ ИТМО, 2012. Вып. 1. 152 с.
5. Гришин А. М. Математические модели лесных пожаров. Томск: Изд-во Томск. ун-та, 1981. 277 с.
Сергей Алексеевич Астафьев
Сведения об авторе — аспирант; Санкт-Петербургский национальный исследовательский
университет информационных технологий, механики и оптики, кафедра мехатроники; E-mail: Rokkolo287@yandex.ru
Рекомендована кафедрой мехатроники
Поступила в редакцию 14.12.12 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
УДК 519.8
С. А. АСТАФЬЕВ
РЕШЕНИЕ ИНТЕГРОДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ТЕПЛОПЕРЕНОСА
В ЗАДАЧЕ МОДЕЛИРОВАНИЯ ПРОЦЕССА РАСПРОСТРАНЕНИЯ ЛЕСНОГО ПОЖАРА
Изложен ход решения интегродифференциального уравнения теплопереноса, используемого как базовое для описания изменения параметров теплового потока в модели распространения лесного пожара. Ключевые слова: моделирование процесса распространения пожара, тепловой поток, интегродифференциальное уравнение теплопереноса, напряжение трения теплового потока.
Введение. На активно охраняемой территории лесного фонда России ежегодно регистрируется от 10 до 35 тыс. лесных пожаров, охватывающих площади 0,5—2,5 млн га. Уничтожение бесценных запасов древесины в результате лесных пожаров сопровождается выгоранием большого объема кислорода и выбросом в атмосферу сажи, копоти, двуокиси углерода, а также радиоактивных частиц [1].
Физика распространения лесного пожара состоит в следующем: тепловой поток, создаваемый низовым пожаром, поднимаясь наклонно по направлению ветра, подогревает кроны деревьев на значительном расстоянии впереди фронта огня. При воспламенении хотя бы одной из крон почти мгновенно воспламеняются и другие, и огонь „скачет“ по подогретым кронам, но затем вне сферы действия подогрева затухает. Когда низовой огонь приближается к фронту пожара, процесс подогрева полога повторяется и опять происходит „скачок огня“. Верховые пожары, выделяя большое количество теплоты, вызывают восходящие потоки продуктов горения и нагретого воздуха и образуют конвективные колонки диаметром в несколько сотен метров. Их поступательное движение совпадает с направлением продвижения фронта пожара. Пламя в середине колонки может подниматься на высоту до 100—120 м. Конвективная колонка увеличивает приток воздуха в зону пожара и порождает ветер, который усиливает горение. Таким образом, низовой пожар стимулирует развитие верхового и наоборот [1].
Основными задачами охраны лесов от пожаров являются их предупреждение, обнаружение, ограничение распространения и тушение. Для эффективной борьбы с лесными пожарами необходим прогноз возможного положения кромки пожара и силы горения.
При математическом моделировании процесса распространения лесного пожара можно выделить четыре группы моделей:
1) модели прогноза скорости распространения пожара; 2) модели прогноза контура и площади пожара;
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
Решение интегродифференциального уравнения теплопереноса
51
3) модели прогноза характеристик теплового течения и массопереноса во фронте и в зоне пожара;
4) общие математические модели, в рамках которых могут быть предсказаны все характеристики во фронте и в зоне пожара (в том числе, поля температур, концентрации и скорости распространения газовых и мелкодисперсных компонентов окружающей среды) [2].
Однако большинство разрабатываемых в настоящее время моделей базируются либо на статистических данных о скорости распространения огня в зависимости от погодных условий, либо на рассмотрении лесного пожара с позиций аэротермохимии. Исходя из этого более перспективным в решении задачи моделирования процесса распространения лесного пожара представляется использование одновременно двух подходов.
При первом — вероятностном — подходе [3, 4] модель прогноза контура и площади пожара к заданному моменту времени строится на основе определения вероятностей распространения огня через дискретные области (по результатам обработки аэрофотоснимков местности) и среднестатистических скоростей распространения пожара в разных направлениях в зависимости от скорости и направления ветра.
При втором подходе — с использованием интегродифференциального уравнения теплопереноса — для расчета параметров теплового потока в зависимости от местоположения относительно фронта пожара используются физические характеристики окружающей среды; кроме того, такой подход позволяет учитывать данные о рельефе местности.
В настоящей статье предложено решение интегродифференциального уравнения теплопереноса для применения его при моделировании процесса распространения лесного пожара.
Интегродифференциальное уравнение. Для моделирования за основу было взято интегродифференциальное уравнение, приведенное в работе [5]. Автор данного уравнения (профессор Томского государственного университета А. М. Гришин) предложил его в качестве математической модели стационарных двумерных квазиравновесных тепловых течений, описывающих процесс распространения лесного пожара:
( )u
∂u ∂x
+
⎛ ⎜ ⎝
vw
−
y
∫
0
∂u ∂x
dy
⎞ ⎟ ⎠
∂u ∂y
=
∂ ∂y
µρ
∂u ∂y
−
sin α Fr
−
1 ρ
∂P ∂x
−
∂P ∂y
y
∫
0
1 ρ
∂
ln P ∂x
dy
,
где и — скорость теплового потока; vw — безразмерная скорость внешнего потока (вдува); µ=18,27·10–6 Па·с — динамическая вязкость воздуха; ρ=1,2 кг/м3 — плотность воздуха;
Fr — критерий Фруда (безразмерная величина) — один из критериев подобия движения жид-
костей и газов; α — угол между касательной к поверхности и горизонтальной плоскостью;
P — давление окружающей среды в зоне пожара.
Система координат xoy, используемая в уравнении, показана на рис. 1.
y
Поле ветра
ymax
Перенос энергии
Приземный слой атмосферы
Фронт пожара
Лес
x
0
Рис. 1
Метод решения. В 80-е годы прошлого века, когда прикладные компьютерные программы для ведения сложных математических расчетов еще не обладали достаточной производительностью и гибкостью и не были широко распространены, автором уравнения было осуществлено его аналитическое решение с использованием различных методов, допущений
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
52 С. А. Астафьев
и замен переменных (например, вывод уравнений методом Прандтля, приближение Буссинеска, переход к переменным Дородницына — Хоуарта, решение уравнений методом Швеца). Как следствие, аналитическое решение является достаточно громоздким, а большое количество преобразований повышает вероятность ошибки; кроме того, к недостаткам такого способа решения можно отнести большие затраты времени на пересчет результатов при изменении входных данных и отсутствие быстрой визуализации результатов вычислений.
На современном этапе степень развития вычислительных систем и приложений позволяет произвести численное решение сложного интегродифференциального уравнения с использованием среды математического моделирования MatLab. Данное программное средство содержит большой набор встроенных функций, предназначенных для решения обыкновенных дифференциальных уравнений. Так, функции ODE реализуют одношаговый явный метод Рунге — Кутты 2-го, 3-го порядка (ode23) и 4-го, 5-го порядков (ode45); неявный метод Рунге — Кутты, использующий формулы обратного дифференцирования 2-го порядка (ode23tb); адаптивный многошаговый метод Адамса — Башворта — Мултона переменного порядка (ode113); одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка (ode23s), и метод трапеций с интерполяцией (ode23t). Функция pdepe служит для решения систем параболических и эллиптических дифференциальных уравнений в частных производных по двум переменным при заданных начальных и граничных условиях. Данная функция преобразует уравнения в частных производных в набор обыкновенных дифференциальных уравнений с использованием дискретизации 2-го порядка точности на основе фиксированного набора заданных пользователем узлов. Интегрирование осуществляется с помощью функции ode15s, которая реализует адаптивный многошаговый метод переменного порядка (от 1 до 5), использующий формулы численного дифференцирования.
Таким образом, главной задачей при решении исходного уравнения в среде MatLab является его запись в форме, адекватной предусмотренной в программном пакете.
Исходные данные, начальные и граничные условия. В качестве параметра vw и координат x, y в процессе решения уравнения использовались безразмерные величины (исходя из предположения, что vwmax=30 м/с — максимальная скорость внешнего потока; ymax=30 м — характерная высота пограничного теплового слоя; xmax=120 м — характерное расстояние теплового воздействия фронта пожара).
В качестве исходных данных условно было задано давление среды в зависимости от координат, выраженное формулой
{ }P(x, y) = 0, 2 + exp −0, 48(x − 0,5)2 − 0,16( y − 0,8)2 .
Р
1,15 1,1 1,05 х 0,4 0,3 0,2 0,1
0
Такой вид зависимости был выбран
1,18 исходя из предположения, что в зоне про-
1,16 грева давление среды экспоненциально рас-
1,14 тет с увеличением высоты (у) над поверхно-
1,12 1,1
стью и расстояния (х) от задней границы фронта пожара, и после достижения максимального значения давление вновь снижает-
1,08 ся до нормального. Заданное изменение
1,06 давления в зависимости от координат гра-
1,04 фически представлено на рис. 2.
0,2 0,4 0,6 0,8 y Рис. 2
1,02 В качестве начального условия (при x=0) было принято экспоненциальное увеличение скорости теплового потока в зависимости от высоты:
u0 ( y) = 1, 0 −1, 0 exp (−2,5 y) .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
Решение интегродифференциального уравнения теплопереноса
53
В качестве нижнего граничного условия (при y=0), согласно данным из работы [5], было
принято нулевое значение скорости теплового потока: uн = 0 , а в качестве верхнего гранично-
го условия (при y=ymax) — экспоненциальное уменьшение скорости теплового потока с увеличением расстояния от фронта пожара:
uв (x) = 1, 0 exp (−17x) .
Файл-функция, задающая граничные условия, содержит два параметра, определяющие
отсутствие „потокопроводности“ на нижней границе пограничного слоя и свободное распро-
странение теплового потока на его верхней границе.
Ход решения. На первом этапе решения исходное интегродифференциальное уравне-
ние теплопереноса было упрощено путем исключения из него интегрального слагаемого
∫ymax
∂u ∂x
dy
.
0
Решение параболического дифференциального уравнения в частных производных про-
изводилось с использованием функции pdepe, имеющей следующий вид:
sol=pdepe(m, @teplofun, @teploinit, @teplobound, y, x),
где u=sol(:,:,1) — искомая скорость теплового потока; m — параметр, соответствующий типу дифференциального уравнения; teplofun — функция, определяющая компоненты решаемого дифференциального
уравнения; teploinit — функция, определяющая начальные условия; teplobound — функция, определяющая граничные условия. Для решения дифференциального уравнения в пакете MatLab оно должно иметь сле-
дующую форму записи:
( ) ( ) ( )c
y,
x,
u,
∂u ∂y
∂u ∂x
=
y−m
∂ ∂y
⎛⎝⎜
ym
f
y,
x,
u,
∂u ∂y
⎞⎟⎠ + s
y,
x,
u,
∂u ∂y
.
Для параболического дифференциального уравнения m=0. Таким образом, исходное уравнение должно быть записано в виде следующей системы:
c = u;
⎫
f
=
⎛ ⎜ ⎝
−vw
+
y
∫
0
∂u ∂x
dy
⎞ ⎟ ⎠
u
+
µρ
∂u ∂y
,
⎪ ⎪⎪ ⎬
s
=
−
sin α Fr
−
1 ρ
∂P ∂x
−
∂P ∂y
y
∫
0
1 ρ
∂
ln P ∂x
⎪ dy.⎪
⎭⎪
В
первую
очередь
были
рассчитаны
частная
производная
давления
по
координате
x
—
∂P ∂x
,
частная производная давления по координате y —
∂P ∂y
и
частная
производная
логарифма
дав-
∫ления по координате x —
∂ ln P ∂x
.
График
значений
интеграла
A
=
ymax
1 ρ
∂ ln P ∂x
dy
представлен
0
на рис. 3.
В результате решения дифференциального уравнения теплопереноса были рассчитаны
значения скорости теплового потока в зависимости от координат u(x,y), частная производная
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
54 С. А. Астафьев
скорости теплового потока по координате x —
∂u ∂x
и частная производная скорости теплового
потока по координате y —
∂u ∂y
(рис. 4, а—в).
А
0,35
0,3 0,25 0,2
0,15
0,1 0,05
0
0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х Рис. 3
а) и
1,2 1 0,8 0,6 0,4 0,2 х
0,3 0,2 0,1
0
0,2 0,4 0,6 0,8
y 1
б)
du dx 40 20 0
–20 х 0,3 0,2 0,1
0,40,60,81 y 0,2 0
в) du dy 2 1,5 1 0,5
х 0,3 0,2 0,1
0 0,2 0,4 0,6 0,8 1 y
Рис. 4
Затем по формулам
τw
=
µρ
∂u ∂y
y=0
и
τe
=
µρ
∂u ∂y
y= ymax
были вычислены напряжение тре-
ния теплового потока на обтекаемой поверхности и напряжение трения теплового потока на
внешней границе пограничного слоя (рис. 5).
Физический смысл параметра τ заключается в том, что при τw = 0 происходит „отрыв“
теплового потока от поверхности, что означает смену типа теплового течения: вместо течения, при котором линии потока вдали от фронта пожара почти параллельны обтекаемой поверхности, образуется конвективная колонка.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
Решение интегродифференциального уравнения теплопереноса
55
Вид зависимости τw(x) вполне согласуется с графиком, представленным в работе [5]:
с увеличением х значение τw уменьшается, причем его резкое уменьшение в пределах фронта
пожара далее переходит в плавное и стремится к нулю.
Вместе с тем следует заметить, что в операциях интегрирования и дифференцирования
(с использованием функций trapz и diff) результат вычислений зависит от шага интегри-
рования (дифференцирования), и поэтому необходимо использовать масштабные коэффи-
циенты.
Одна из особенностей решения состоит в том, что вычисления производились методом
последовательных приближений. Так, в первом цикле вычислений рассчитывалось значение
∫интеграла
B
=
ymax
∂u ∂x
dy
,
которое
затем
использовалось
в последующей
итерации. В резуль-
0
тате после нескольких расчетных циклов значение интеграла, вычисленное по окончании
программы, сходится к значению, заданному на входе. График значений интеграла В пред-
ставлен на рис. 6.
τw, τe, о.е.
0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х
у у=ymax
2 –1
1,5 –2 –3
1 –4
0,5 –5
0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х
В
Рис. 5
Рис. 6
Заключение. Представленное решение интегродифференциального уравнения теплопе-
реноса осуществлено при условно принятых начальных и граничных условиях, подобранных
таким образом, что обеспечивается устойчивое решение уравнения. Ход вычислений отдель-
ных параметров имеет графическую визуализацию. Впоследствии необходимо произвести
согласование параметров модели с реальными условиями распространения огня путем прове-
дения натурных экспериментов и исследований параметров реальных пожаров.
Таким образом, использование уравнения теплопереноса в основе модели распростра-
нения лесного пожара позволяет в перспективе обеспечить возможность прогнозирования
эволюции пожара путем учета конкретных физических параметров окружающей среды. Сре-
ди параметров модели наиболее значимым, с точки зрения предсказания развития лесного
пожара во времени, представляется напряжение трения теплового потока. Дальнейшие иссле-
дования целесообразно направить на выявление характерных изменений данного параметра в
зависимости от внешних условий.
СПИСОК ЛИТЕРАТУРЫ
1. Воробьев Ю. Л., Акимов В. А., Соколов Ю. И. Лесные пожары на территории России: Состояние и проблемы / Под общ. ред. Ю. Л. Воробьева. М.: ДЭКС-ПРЕСС, 2004. 312 с.
2. Гришин А. М. О математическом моделировании природных пожаров и катастроф // Вестн. Томск. гос. ун-та. Математика и механика. 2008. № 2(3).
3. Астафьев С. А., Лысенко Д. Ю., Широков А. С. Моделирование процесса распространения лесного пожара с применением теории перколяции // Изв. вузов. Приборостроение. 2012. Т. 55, № 6. С. 70—74.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7
56 Ю. В. Баёва, Е. В. Лаповок, С. И. Ханков
4. Астафьев С. А. Применение вероятностного подхода в задаче моделирования распространения лесного пожара // IX Всерос. конф. молодых ученых; V сессия науч. школы „Проблемы механики и точности в приборостроении“: Сб. докл. СПб: НИУ ИТМО, 2012. Вып. 1. 152 с.
5. Гришин А. М. Математические модели лесных пожаров. Томск: Изд-во Томск. ун-та, 1981. 277 с.
Сергей Алексеевич Астафьев
Сведения об авторе — аспирант; Санкт-Петербургский национальный исследовательский
университет информационных технологий, механики и оптики, кафедра мехатроники; E-mail: Rokkolo287@yandex.ru
Рекомендована кафедрой мехатроники
Поступила в редакцию 14.12.12 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2013. Т. 56, № 7