Содержание к диссертации
Введение
Глава 1. Обзор исследований по моделированию излучения при лесных пожарах 13
1.1. Современные подходы к моделированию пожаров 13
1.2. Существующие подходы к моделированию процесса переноса лучистой энергии 27
1.3. Методы решения уравнений химической кинетики 37
Глава 2. Математическое моделирование динамики пожара с учётом излучения 43
2.1. Физическая и математическая постановка задачи 43
2.2. Теоретические основы дифференциальной модели излучения в диффузионном приближении 50
2.3. Теоретические основы дифференциальной модели излучения в диффузионно-волновом приближении 57
2.4. Результаты моделирования излучения по диффузионному и диффузионно-волновому приближениям дифференциальной модели и их анализ 61
2.5. Рекомендации по моделированию излучения при пожарах 63
Глава 3. Алгоритмы и численное моделирование динамики ландшафтных пожаров 64
3.1. Основные вычислительные сложности при моделировании динамики ландшафтных лесных пожаров 64
3.2. Особенности численного моделирования 65
3.3. Результаты моделирования динамики пожара под воздействием внешнего поля скоростей с различными моделями излучения и их анализ 76
Заключение 93
Список использованных источников
- Существующие подходы к моделированию процесса переноса лучистой энергии
- Методы решения уравнений химической кинетики
- Теоретические основы дифференциальной модели излучения в диффузионном приближении
- Особенности численного моделирования
Существующие подходы к моделированию процесса переноса лучистой энергии
Перенос тепла от фронта пожара к ещё не тронутым огнём лесным горючим материалам является одной из составляющих, необходимой для распространения пожара. Процессы конвекции, кондукции и излучения составляют теплоперенос. Как известно [38], радиационный теплоперенос преобладает над конвективным уже при 400 К, однако его моделирование существенно сложнее вследствие распространения лучистой энергии на несколько метров от нагретой поверхности, тогда как конвективный и диффузионный теплоперенос имеют место только при локальном перепаде температур. В ряде моделей, таких как [21, 25, 26, 39] основным фактором теплопереноса считается излучение.
Уравнение переноса как часть теории излучения было получено благодаря трудам Schuster [40] ещё на рубеже XIX и XX веков. В 50-х годах прошлого века, исследования уравнений переноса стали самостоятельным разделом математической физики, что позволило создать новые более эффективные методы решения интегральных и интегро-дифференциальных уравнений, описанные в монографиях [41-45]. В работах 60-80 годов делаются попытки вывести соотношения теории переноса на основе статистической теории волновых полей [46-52]. В работах [53-56] моделируются процессы излучения, возникающие при аварийных ситуациях. В [53] при построении полуэмпирической модели горения ракетного топлива принято равномерное тепловыделение в объеме огненного шара и нормальное распространение пламени на его поверхности. Существуют различные подходы при моделировании процесса переноса лучистой энергии, так например, в работе [54] используется максимально упрощенная схема, в основе которой лежат эмпирические зависимости излучающей способности огненного шара от массы аварийного выброса топлива. Кроме того, в данной работе предложена методика оценки безопасных расстояний при пожарах на опасных производствах в результате теплового излучения. В настоящее время существует множество моделей излучения. При этом основное уравнение радиационного переноса [57] можно записать в виде ?L = -(ka+ks)I + ka + %- J / (Q,Q )/(QVQ as n 4 я 4Я где / - интенсивность радиационного излучения в направлении Q; s -расстояние в направлении Q; Eg = аГё - энергия, излучаемая абсолютно черным газом при температуре газа Tg; ка и ks - коэффициенты поглощения и рассеяния; Р(ї, О. ) - вероятность того, что излучение в направлении Q после рассеяния попадет в телесный угол dQ. в окрестности направления Q. Это уравнение необходимо интегрировать по всем направлениям и длинам волн. Более того, в задачах, предполагающих изменение мощности и конфигурации источников излучения по времени необходимо решать такое уравнение на каждом шаге по времени. Для большинства практических задач точное решение невозможно [58], вместо него разработано несколько приближенных методов, которые и используются для моделирования динамики пожаров в помещениях.
Если разделить пространственное и угловое распределение интенсивности излучения, задачу можно существенно упростить. Этот подход используется в "потоковых методах" [59]. Если предположить, что спектральная интенсивность постоянна в пределах заданных интервалов телесного угла, то уравнение радиационного переноса сводится к нескольким связанным между собой обыкновенным линейным дифференциальным уравнениям относительно осредненных по пространству интенсивностей или потоков излучения. Однако этот метод имеет ряд недостатков, среди которых одним из главных, применительно к пожарам, является неточность метода при моделировании радиационного переноса под углом к декартовой сетке.
Модель, разработанная Локвудом и Шахом [60], преодолевает основной недостаток потоковых методов. Для данной модели характерны некоторые черты методов Монте-Карло, а именно прохождение "лучей" электромагнитного излучения через вычислительную область между границами. Однако в отличие от методов Монте-Карло, где направления лучей генерируются случайным образом, в этой модели они выбираются предварительно, таким же образом, как выбирается расположение гидродинамической сетки. Метод включает в себя решение уравнения радиационного переноса вдоль путей этих лучей, выбираемых обычно таким образом, чтобы они приходили в центры граничных поверхностей гидродинамических контрольных объемов. В работе [61] также используется модель излучения на основе метода Монте-Карло, основанный на моделировании распространения фотонов по случайным направлениям. При использовании методов Монте-Карло к решению распространение излучения представляется в виде случайной марковской цепи столкновений фотонов или дискретных порций энергии с веществом. Автором учитывается оптическая неоднородность среды. Данный метод может давать приемлемое время расчёта при моделировании очага пожара, размеры которого достаточно малы. В случае если излучающая поверхность пламени имеет большие размеры, применение данного метода затруднительно в силу того, что необходимо вычислить распространение излучения между большим количеством элементов как излучающей, так и поглощающей поверхности.
Методы решения уравнений химической кинетики
Зачастую физические явления включают в себя как быстро, так и медленно протекающие процессы. Например, в реакционно-способной системе, могут существовать несколько химических реакций с сильно различающимися скоростями, что требует наложения жёстких ограничений на шаг по времени. Свойство А-устойчивости [86, 87] позволяет снять такие ограничения и поэтому важно для численных схем, применяемых для решения жёстких систем. Иными словами, А-устойчивый метод гарантирует ограниченность погрешности для жёсткой системы при любом шаге по времени. Одним из наиболее применимых подходов к получению схем является расщепление [88, 89]. Моделирование химических реакций -является важной частью решения задач горения. Точность расчётов в большой степени зависит от адекватности применяемых моделей и численных схем. Самым простым способом учёта физико-химических процессов в неоднородной среде является использование предположения о мгновенности происходящих реакций в рамках одного шага по времени, при наличии реагирующих компонент в достаточном количестве в ячейке расчётной области. Для описания физико-химических процессов, как правило, используют систему обыкновенных дифференциальных уравнений следующего вида —t = fk(t,C),k = l,2,...,N. at Как было отмечено в [86], для данной системы уравнений химической кинетики из функции скорости изменения концентраций можно выделить скорость образования компоненты и поглощения, пропорциональную концентрации этой компоненты, что приведёт к системе уравнений в следующей форме = Vla(Ca)-CaV2a(Ca),k-l,2,...,N,Vla 0,V2a 0. at Как было показано Зельдовичем [90], для данной системы уравнений существует только одно стационарное решение в к - мерном пространстве. Нетрудно заметить, что правая часть состоит из положительного слагаемого, отвечающего за образование компонента и отрицательного слагаемого, отвечающего за его расход, причём предполагается его пропорциональность концентрации. Учитывая выше сказанное, можно воспользоваться подходом, аналогичным способу линеаризации источникового члена в методе Патанкара [59]. В силу того, что при решении жёстких систем обыкновенных дифференциальных уравнений собственные значения матрицы Якоби для правых частей системы имеют большой разброс, возникают трудности при их численном решении на основе классических методов, таких как Рунге Кутты и Адамса [91], связанные с жёсткими требованиями к величине шага по времени для обеспечения устойчивости. Иными словами, жесткость системы уравнений обуславливает высокие требования к устойчивости используемых численных схем.
Быстротекущие реакции в жёсткой задаче Коши даже на коротких интервалах проявляют некоторые свойства, характерные для бесконечных интервалов. Это делает возможным применение к ним методов аппроксимации на бесконечных интервалах. Ввиду сложности рассматриваемых задач, применяемые к ним схемы должны удовлетворять следующим требованиям: - численное решение должно быть А-устойчивым; - концентрации должны быть неотрицательными; - сумма всех расчётных значений концентраций должна быть равна единице; - метод должен быть экономичным. Требование экономичности диктует необходимость использовать одношаговые разностные схемы. В [86] представлена схема Кранка-Николсона, основанная на известном методе трапеции для аппроксимации неоднородной правой части системы обыкновенных дифференциальных уравнений. Схема Кранка-Николсона имеет следующий вид: С - «=у/«(0 + /««)] Данная схема является явно-неявной, и для её решения часто используется метод простой итерации (см., например, [29, 86, 93-99]). Итерационный процесс представляется с помощью формулы вида: /7+1 П I у / П+\ \ , Г / Л\ uak Ua =—[fa(Uak-0 + fa(Ua)\ где &=1,2,3 — номер итерации, а и" = иа. Итерационный процесс, выполняемый на каждом шаге по времени, продолжается до достижения заданной точности или предельного числа итераций. Известно, что для сходимости таких методов требуется использовать критерий устойчивости.
В другом подходе [55, 56] скорость горения ограничивается процессом смешения компонентов горючей смеси, притом что химическая реакция считается бесконечно быстрой.
Третий подход основан на использовании теории тепломассообмена для построения инженерных моделей горения смеси газов [59]. Используя метод Патанкара [59] можно получить неявную схему.
Теоретические основы дифференциальной модели излучения в диффузионном приближении
Диффузионное приближение [72, 92] основано на модификации закона распространения лучистой энергии (1.3) путём замены мгновенного характера распространения лучистой энергии на диффузионный с фиктивным коэффициентом диффузии vU37_1. При этом в уравнение (1.2) вносить изменения не требуется.
Для упрощения расчёта распространения излучения около границ целесообразно уменьшить мощность излучения от среды на величину, равную мощности излучения при начальной температуре. Такой приём позволяет уменьшить погрешность для граничных условий свободной границы за счёт того, что пожар рассматривается на удалении от неё, а излучение от нагретой среды ослабевает с расстоянием. С учётом вышесказанного, диффузионное приближение будет иметь следующий вид: (( з \ Т id(p5cp5UT)id(p5cp5WT) = dt дх dz (2.15) TP. p.+PsCpS Vi=l JJ_ дх\ дх J dz\ dz J q2R2 + q3R3 + q5R5 + k,(cUR - 4 j(T4 - Te4)\ R dUR _ д dt dx с dUR \Ъкъ дх + dz с dU \Ъкъ dz k5{cUR-4a(T4e4)),(2.\6) VUSJ, - параметр, характеризующий искусственную диффузию в диффузионной и диффузионно-волновой моделях излучения; сшл - параметр характеризующий процесс волнового распространения энергии в диффузионно-волновой модели.
В обеих постановках предполагается наличие однородных граничных условий первого, второго или третьего рода. Для анализа распространения энергии предположим, что излучение имеет место на отрезке времени t e[t],t2],0 tl t2. Интегрируя уравнения (2.17), (2.18) по времени t є[0,оо), с учётом &R = \и Rdt получим: f д_ dx d$K \Ъкъ дх VU3J$R со Q о дх г д_ dz ( дЗ„ \Ъкъ dz d$D \Ъкг dx + dz k{c3R) = )4ksaT dt + ks{c$R) d& + ks{c$R) yikz dz , (2.19) (2.20) \4kscrT dt + ks{clR)
Таким образом, различие между уравнениями сводится к члену vU3nG-R Так как в начальный момент излучения ещё нет, а на бесконечный момент оно уже рассеялось (что будет показано ниже в данном параграфе), следовательно, величина $Лне зависит от значения vU3Jl, иными словами, распределение энергии по пространству не зависит от фиктивного параметра v изл Вследствие линейности уравнений (2.17) и (2.18), излучение от источника любой конфигурации можно представить в виде суперпозиции излучений точечных источников. Используя это можно исследовать влияние диффузионного приближения на примере излучения мгновенного точечного источника при нулевых начальных условиях.
Отметим, для справедливости результатов, полученных на основе подхода, описанного в формулах (2.28) - (2.30), необходимо, чтобы численная схема не допускала появления отрицательных значений искомых величин. В противном случае, полученная оценка (2.30) не будет справедливой. В отличие от большинства упрощений, огрубляющих модели, диффузионное приближение не создаёт принципиально непреодолимой погрешности, давая возможность повышать точность за счёт роста количества вычислений, что характерно для численных схем. Коэффициент vmn в уравнении (2.17) не обусловлен физическими параметрами рассматриваемой предметной области, поэтому его следует выбирать исходя из целесообразности для численного расчёта. С одной стороны, как было показано выше, коэффициент \изл должен быть мал настолько, насколько это возможно. С другой стороны, его уменьшение приводит к росту скорости диффузии, и уменьшению размера требуемого шага по времени при использовании явных схем. Увеличение коэффициента диффузии при численных расчётах приводит к эффектам, аналогичным эффектам от увеличения шага по времени.
Особенности численного моделирования
Для учёта влияния геометрии местности на невозмущённое поле скоростей решается стационарная задача обтекания рельефа. При этом температура и состав газовой фазы считаются постоянными, а химические реакции не учитываются. После достижения заданной точности, получившиеся поля скоростей и давлений используются как начальные условия для задачи горения. В случае однородного рельефа (равнина) скорость газовой фазы пропорциональна высоте.
Расщепление по физическим процессам является частью метода крупных частиц [101], согласно которому на каждом шаге по времени выполняются поочерёдно Эйлеров и Лагранжев этапы. Этим этапам, соответственно сопоставлены системы уравнений (3.3) и (3.4). Особенностью алгоритма, представленного в настоящей диссертации, является выделение процесса распространения лучистой энергии при расщеплении, что соответствует системе уравнений (3.5).
Система уравнений (3.3) получается из уравнений (2.1) - (2.4), если в них опустить дивергентные слагаемые. С физической точки зрения такая система описывает процесс изменения параметров системы за счет работы сил давления, действующих на границе области, эффекты химико-физических процессов и диффузии. В свою очередь, система (3.4), содержащая дивергентные слагаемые, отвечает процессу конвективного переноса газодинамических величин.
Для расчётов использовались конфигурации, показанные на рис. 2.1 — 2.3. На начальном этапе программы выполняется дискретизация рельефа и области РГМ [103-106]. Трапеция, формирующая рельеф, описывается как кусочно-линейная функция высоты от горизонтальной координаты. Каждый элемент расчётной области, имеющий высоту ниже данной трапеции, считается элементом препятствия. Если трапеция имеет отрицательную высоту (для случая оврага), то для корректного моделирования рельефа за нулевой уровень в численной реализации принимается нижняя точка оврага. После формирования рельефа вычисляется высота леса в шагах расчётной сетки и заполняется соответствующее количество ячеек. Размер шага по координате выбирается таким образом, чтобы высота леса была кратна ему, в противном случае она будет фактически округлена.
Для реализации границы, связанной с неоднородным рельефом, каждой ячейке расчётной области ставится в соответствие В,у, значение которой равно нулю, если (i,j) содержит преграду и единицу в противном случае. В ячейках, содержащих препятствие, не требуется проводить вычисления. Предполагается, что соотношения (2.14) выполняются на каждой границе между свободной ячейкой и ячейкой, содержащей преграду. Реализация граничных условий выполняется за счёт комбинирования их дискретного аналога с общей численной схемой, что позволяет упростить обсчёт угловых точек.
В рамках рассматриваемой системы уравнений и методов её решения, на каждом шаге по времени численную схему можно записать в виде шаблона (3.6). Неизвестные значения на половинных узлах вычисляются как среднеарифметическое известных значений вокруг них. Уравнение (3.6) интерпретируется как изменение величины % за счёт потоков через границы с соседними ячейками и источникового члена. С=&+ Де„2„ -f,: ,few, -с)+ + №., -$:,-m)+ Mj - ,./ )+ , где а , (А 0...4) - коэффициенты, - обобщённая переменная, соответствующая значениям искомых величин, " - её дискретное представление, i,j - пространственные дискретные координаты, т- номер шага по времени. В рассматриваемой задаче существует два вида граничных условий на твёрдых границах.
На первом этапе алгоритма выполнения расчётного шага программа вычисляет скорости химико-физических процессов R], R2, R3, а также массовую скорость выделения компонентов газовой фазы R5!, R52.
Расчёты проводятся только в точках, где температура выше некоторого критического значения, выбираемого исходя из допустимости пренебрежения эффектами химико-физических процессов при более низких температурах. В программной реализации эта температура составляет 400 К. Данный подход позволяет уменьшить объём таких трудоёмких операций, как вычисление экспоненты.
Особенностью выбранной модели лесного пожара является расчёт скоростей химико-физических процессов согласно закону Аррениуса, вследствие чего при высоких температурах возможно возникновение скоростей реакций, при которых расход вещества за шаг времени превышает его запас в ячейке. При использовании явной схемы такая ситуация приводит к отрицательным концентрациям, что не только является нефизичным результатом, но и может привести к полному разрушению расчётной схемы.
