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

ТЕПЛОВЫЕ РЕЖИМЫ ПРИБОРОВ И СИСТЕМ

УДК 519.8

С. А. Астафьев РЕШЕНИЕ

ИНТЕГРОДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ ТЕПЛОПЕРЕНОСА В ЗАДАЧЕ МОДЕЛИРОВАНИЯ ПРОЦЕССА РАСПРОСТРАНЕНИЯ ЛЕСНОГО ПОЖАРА

Изложен ход решения интегродифференциального уравнения теплопереноса, используемого как базовое для описания изменения параметров теплового потока в модели распространения лесного пожара.

Ключевые слова: моделирование процесса распространения пожара, тепловой поток, интегродифференциальное уравнение теплопереноса, напряжение трения теплового потока.

Введение. На активно охраняемой территории лесного фонда России ежегодно регистрируется от 10 до 35 тыс. лесных пожаров, охватывающих площади 0,5—2,5 млн га. Уничтожение бесценных запасов древесины в результате лесных пожаров сопровождается выгоранием большого объема кислорода и выбросом в атмосферу сажи, копоти, двуокиси углерода, а также радиоактивных частиц [1].

Физика распространения лесного пожара состоит в следующем: тепловой поток, создаваемый низовым пожаром, поднимаясь наклонно по направлению ветра, подогревает кроны деревьев на значительном расстоянии впереди фронта огня. При воспламенении хотя бы одной из крон почти мгновенно воспламеняются и другие, и огонь „скачет" по подогретым кронам, но затем вне сферы действия подогрева затухает. Когда низовой огонь приближается к фронту пожара, процесс подогрева полога повторяется и опять происходит „скачок огня". Верховые пожары, выделяя большое количество теплоты, вызывают восходящие потоки продуктов горения и нагретого воздуха и образуют конвективные колонки диаметром в несколько сотен метров. Их поступательное движение совпадает с направлением продвижения фронта пожара. Пламя в середине колонки может подниматься на высоту до 100—120 м. Конвективная колонка увеличивает приток воздуха в зону пожара и порождает ветер, который усиливает горение. Таким образом, низовой пожар стимулирует развитие верхового и наоборот [1].

Основными задачами охраны лесов от пожаров являются их предупреждение, обнаружение, ограничение распространения и тушение. Для эффективной борьбы с лесными пожарами необходим прогноз возможного положения кромки пожара и силы горения.

При математическом моделировании процесса распространения лесного пожара можно выделить четыре группы моделей:

1) модели прогноза скорости распространения пожара;

2) модели прогноза контура и площади пожара;

3) модели прогноза характеристик теплового течения и массопереноса во фронте и в зоне пожара;

4) общие математические модели, в рамках которых могут быть предсказаны все характеристики во фронте и в зоне пожара (в том числе, поля температур, концентрации и скорости распространения газовых и мелкодисперсных компонентов окружающей среды) [2].

Однако большинство разрабатываемых в настоящее время моделей базируются либо на статистических данных о скорости распространения огня в зависимости от погодных условий, либо на рассмотрении лесного пожара с позиций аэротермохимии. Исходя из этого более перспективным в решении задачи моделирования процесса распространения лесного пожара представляется использование одновременно двух подходов.

При первом — вероятностном — подходе [3, 4] модель прогноза контура и площади пожара к заданному моменту времени строится на основе определения вероятностей распространения огня через дискретные области (по результатам обработки аэрофотоснимков местности) и среднестатистических скоростей распространения пожара в разных направлениях в зависимости от скорости и направления ветра.

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

В настоящей статье предложено решение интегродифференциального уравнения тепло-переноса для применения его при моделировании процесса распространения лесного пожара.

Интегродифференциальное уравнение. Для моделирования за основу было взято ин-тегродифференциальное уравнение, приведенное в работе [5]. Автор данного уравнения (профессор Томского государственного университета А. М. Гришин) предложил его в качестве математической модели стационарных двумерных квазиравновесных тепловых течений, описывающих процесс распространения лесного пожара:

и % +

дх

(

\

и

-УФУ 0 их

ди = д / ди ду ду ^ ду

^п а _ 1 дР _дРу1 д 1п Р р дх ду 0 р дх

Бг

Ф,

где и — скорость теплового потока; у№ — безразмерная скорость внешнего потока (вдува); ц=18,27-10-6 Па-с — динамическая вязкость воздуха; р=1,2 кг/м3 — плотность воздуха; Бг — критерий Фруда (безразмерная величина) — один из критериев подобия движения жидкостей и газов; а — угол между касательной к поверхности и горизонтальной плоскостью; Р — давление окружающей среды в зоне пожара.

Система координат хоу, используемая в уравнении, показана на рис. 1.

Перенос энергии

Поле ветра

Приземный слой атмосферы

0 77// П~ГГП Г-ГТТТТ1—ГГТТ~П~ГГГТ1 гт

Рис. 1

Метод решения. В 80-е годы прошлого века, когда прикладные компьютерные программы для ведения сложных математических расчетов еще не обладали достаточной производительностью и гибкостью и не были широко распространены, автором уравнения было осуществлено его аналитическое решение с использованием различных методов, допущений

и замен переменных (например, вывод уравнений методом Прандтля, приближение Буссине-ска, переход к переменным Дородницына — Хоуарта, решение уравнений методом Швеца). Как следствие, аналитическое решение является достаточно громоздким, а большое количество преобразований повышает вероятность ошибки; кроме того, к недостаткам такого способа решения можно отнести большие затраты времени на пересчет результатов при изменении входных данных и отсутствие быстрой визуализации результатов вычислений.

На современном этапе степень развития вычислительных систем и приложений позволяет произвести численное решение сложного интегродифференциального уравнения с использованием среды математического моделирования 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Д8 1,16 1,14 1,12 1,1 1,08 1,06 1,04 1,02

1,15

0,4 0,6 Рис. 2

и.

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

В качестве начального условия (при х=0) было принято экспоненциальное увеличение скорости теплового потока в зависимости от высоты: 0(У) = 1,0 -1,0 ехр (-2,5у) .

В качестве нижнего граничного условия (при у=0), согласно данным из работы [5], было принято нулевое значение скорости теплового потока: ин = 0, а в качестве верхнего граничного условия (при у=утах) — экспоненциальное уменьшение скорости теплового потока с увеличением расстояния от фронта пожара:

ив (х) = 1,0 ехр (-17х) .

Файл-функция, задающая граничные условия, содержит два параметра, определяющие отсутствие „потокопроводности" на нижней границе пограничного слоя и свободное распространение теплового потока на его верхней границе.

Ход решения. На первом этапе решения исходное интегродифференциальное уравнение теплопереноса было упрощено путем исключения из него интегрального слагаемого

утах

Л"*.

дх 0

Решение параболического дифференциального уравнения в частных производных производилось с использованием функции pdepe, имеющей следующий вид:

sol=pdepe(m, @teplofun, @teploinit, @teplobound, у, x),

где

u=sol(:,:,1) — искомая скорость теплового потока; т — параметр, соответствующий типу дифференциального уравнения; teplofun — функция, определяющая компоненты решаемого дифференциального уравнения;

teploinit — функция, определяющая начальные условия; teplobound — функция, определяющая граничные условия.

Для решения дифференциального уравнения в пакете Ма1ЬаЬ оно должно иметь следующую форму записи:

(x, %)=y"" £ (y"/(y, x «■ |)>* (y, x, «•

du

dW

Для параболического дифференциального уравнения т=0.

Таким образом, исходное уравнение должно быть записано в виде следующей системы:

c = u;

í y ~ \

f = _ v +\dUdy w 1 dx

\ o°x J

u,

5 = sin a _1 dP -dP y 1 d ln P dy Fr p dx dy 0 p dx

dP

В первую очередь были рассчитаны частная производная давления по координате x — д-,

dP

частная производная давления по координате y — dy и частная производная логарифма дав- * ~ ymax , _

dln P Г 1 dln P

ления по координате x — dx ■ График значений интеграла A = J p dx dy представлен

o p

на рис. 3.

В результате решения дифференциального уравнения теплопереноса были рассчитаны значения скорости теплового потока в зависимости от координат u(x,y), частная производная

скорости теплового потока по координате х — -д— и частная производная скорости теплового потока по координате у — -ду (рис. 4, а—в).

а) и 1,2

ния теплового потока на обтекаемой поверхности и напряжение трения теплового потока на внешней границе пограничного слоя (рис. 5).

Физический смысл параметра т заключается в том, что при тw = 0 происходит „отрыв"

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

Вид зависимости тw(х) вполне согласуется с графиком, представленным в работе [5]:

с увеличением х значение тw уменьшается, причем его резкое уменьшение в пределах фронта

пожара далее переходит в плавное и стремится к нулю.

Вместе с тем следует заметить, что в операциях интегрирования и дифференцирования (с использованием функций 1;гарЕ и ^ГГ) результат вычислений зависит от шага интегрирования (дифференцирования), и поэтому необходимо использовать масштабные коэффициенты.

Одна из особенностей решения состоит в том, что вычисления производились методом последовательных приближений. Так, в первом цикле вычислений рассчитывалось значение

утах

интеграла В = | дхс^У, которое затем использовалось в последующей итерации. В резуль-

0

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

о.е

0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х

0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 х Рис. 5

-1 -2 -3 -4 -5 В

Рис. 6

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

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

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

1. Воробьев Ю. Л., Акимов В. А., Соколов Ю. И. Лесные пожары на территории России: Состояние и проблемы / Под общ. ред. Ю. Л. Воробьева. М.: ДЭКС-ПРЕСС, 2004. 312 с.

2. Гришин А. М. О математическом моделировании природных пожаров и катастроф // Вестн. Томск. гос. ун-та. Математика и механика. 2008. № 2(3).

3. Астафьев С. А., Лысенко Д. Ю., Широков А. С. Моделирование процесса распространения лесного пожара с применением теории перколяции // Изв. вузов. Приборостроение. 2012. Т. 55, № 6. С. 70—74.

0

56

Ю. В. Баёва, Е. В. Лаповок, С. И. Ханков

4. Астафьев С. А. Применение вероятностного подхода в задаче моделирования распространения лесного пожара // IX Всерос. конф. молодых ученых; V сессия науч. школы „Проблемы механики и точности в приборостроении": Сб. докл. СПб: НИУ ИТМО, 2012. Вып. 1. 152 с.

5. Гришин А. М. Математические модели лесных пожаров. Томск: Изд-во Томск. ун-та, 1981. 277 с.

Сведения об авторе

Сергей Алексеевич Астафьев — аспирант; Санкт-Петербургский национальный исследовательский

университет информационных технологий, механики и оптики, кафедра мехатроники; E-mail: [email protected]

Рекомендована кафедрой Поступила в редакцию

мехатроники 14.12.12 г.

УДК 520.224.2. 224.4

Ю. В. Баёва, Е. В. Лаповок, С. И. Ханков

МЕТОД ПОДДЕРЖАНИЯ ЗАДАННОГО ТЕМПЕРАТУРНОГО ДИАПАЗОНА

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

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

Ключевые слова: космический аппарат, степень черноты, коэффициент поглощения солнечного излучения, коэффициент облученности, тепловой режим, теплообмен излучением.

При движении космического аппарата (КА, далее — объект) по круговой орбите вокруг Земли с периодическим заходом в тень и выходом на освещенный Солнцем участок траектории его средняя температура может колебаться в больших пределах. Во многих случаях при наружном размещении на борту космического аппарата оптических систем (ОС) требуется обеспечить минимальные колебания температур их элементов относительно некоторого среднего значения Т0.

В настоящей статье представлена аналитическая методика расчета параметров КА (с выявлением наиболее значимых), обеспечивающих минимальное колебание его температуры. Суть методики заключается в расчете температур КА на теневом и солнечном участках траектории в стационарном тепловом режиме при условии, что на теневом участке роль Солнца как источника тепловыделений выполняет собственный поверхностный источник тепловыделений объекта.

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

— расчет проводится для экваториальной круговой орбиты, которая характеризуется максимальным изменением температур при двух вариантах положения объекта на линии центр Солнца — центр Земли: между Солнцем и Землей и в тени Земли;

— рассматриваются объекты сферической и цилиндрической конфигурации, что соответствует типовым конфигурациям КА и ОС;