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

Статья в формате PDF

"Техническое Руководство к FDS" в формате PDF (англ.)

Выходные величины

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

• Температуру газа
• Скорость газа - Концентрацию частиц газа (водяной пар, $CO_{2},CO,N_{2}$)
• Концентрация дыма и оценки выдимости
• Давление - Скорость выделения теплоты на единицу объема
• Доля в смеси (или отношение воздух/топливо)
• Плотность газа
• Масса водяных капель на единицу объема

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

• Поверхностная и внутренняя температура
• Тепловой поток - и радиационный и конвективный
• Скорость горения
• Масса водяных капель на единицу площади

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

• Общую скорость выделения теплоты (HRR)
• Время активации спринклеров и детекторов
• Потоки массы и энергии через проемы и твердые тела

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

Основные уравнения и процедура решения

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

Основные уравнения

Эта секция знакомит с основными уравнениями сохранения массы, импульса и энергии для Ньютоновской жидкости. Это те же самые уравнения, которые могут быть найдены почти в любой книге по динамике жидкостей или CFD. Хорошая ссылка с описанием уравнений, указаниями и различными приближениями - это Anderson[1]. Обратите внимание, что это набор дифф.уравнений в частных производных, содержащий шесть уравнений для шести неизвестных, все функции в трех пространственных координатах и времени: плотность $\rho$, три компоненты скорости $\mathbf{\vec{u}}=[u,v,w]^{T}$, температура $T$ и давление $p$.
 

Перенос массы и веществ

Сохранение массы может быть выражено либо в терминах плотности $\rho$:
\begin{eqnarray} \frac{\partial\rho}{\partial t}+\nabla\cdot\rho\mathbf{\vec{u}} & = & \dot{m}_{b}^{'''} \end{eqnarray} либо, в терминах индивидуальных газовых компонент, $Y_{\alpha}$:
\begin{eqnarray} \frac{\partial}{\partial t}(\rho Y_{\alpha})+\nabla\cdot\rho Y_{\alpha}\mathbf{u=}\nabla\cdot\rho D_{\alpha}\nabla Y_{\alpha}+\dot{m}_{\alpha}^{'''}+\dot{m}_{b,\alpha}^{'''} \end{eqnarray} где $\dot{m}_{b}^{'''}=\sum_{\alpha}\dot{m}_{b,\alpha}^{'''}$ - скорость образования газовых компонент из-за испарения капель или частиц. Суммируя эти уравнения по всем компонентам, получим оригинальное уравнение сохранения массы, потому что $\sum Y_{\alpha}=1\textrm{ и }\sum\dot{m}_{\alpha}^{'''}=0\textrm{ и }\sum\dot{m}_{b,\alpha}^{'''}=\dot{m}_{b}^{'''}$, по определению, и потому что предполагается, что $\sum\rho D_{\alpha}\nabla Y_{\alpha}=0$. Последнее утверждение, вообще, не верно. Однако, уравнения переноса решаются для общей массы за исключением одного из компонент, предполагая что диффузионный коэффициент подразумеваемых компонент выбран так, что сумма всех диффузионных потоков равна нулю.
 

Перенос импульса

Уравнение сохранениия импульса записывается как:
\begin{eqnarray} \frac{\partial}{\partial t}(\rho\mathbf{u})+\nabla\cdot\rho\mathbf{uu}+\nabla p=\rho\mathbf{g}+\mathbf{f}_{b}+\nabla\cdot\tau_{ij} \end{eqnarray} Здесь $\mathbf{\boldsymbol{\mathbf{uu}}}$ - двухэлементный тензор. В матричной записи с $\mathbf{\vec{u}}=[u,v,w]^{T}$, двухэлементный тензор представлен тензорным произведением векторов $\mathbf{u}$ и $\mathbf{u}^{T}$. Термин $\nabla\cdot\rho\mathbf{uu}$- соответственно вектор, полученный применением векторного оператора $\nabla=(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z})$ к тензору. Сила $\mathbf{f}_{b}$ в уравнении импульса представляет внешние силы, такие как сопротивление движению через капли жидкости. Тензор напряжения $\tau_{ij}$ определяется как:
\begin{eqnarray*} \tau_{ij} & = & \mu\,\Big(\,2\mathbf{S}_{ij}-\frac{2}{3}\delta_{ij}(\nabla\cdot\mathbf{u})\,\Big);\qquad\delta_{ij}=\begin{cases} \begin{array}{cc} 1\quad & i=j\\ 0\quad & i\neq j \end{array}\end{cases}; \end{eqnarray*} \begin{equation} S_{ij}=\frac{1}{2}\Big(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\Big)\qquad i,j=1,2,3 \end{equation} Элемент $\mathbf{S}_{ij}$ - симметричный тензор скорости деформации, выписанный с использованием традиционной тензорной записи. Символ $\mu$ - динамическая вязкость газа.
Полное вычисление может быть выполнено либо с применением Прямого Численного Моделирования (DNS), в котором рассеивающие (диссипативные) члены вычисляются напрямую, либо как Масштабное Моделирование Вихрей (LES), в котором масштабные вихри вычисляются напрямую и моделируются процессы рассеивания в подсеточных масштабах. Численные алгоритмы разработаны так что LES становится DNS-ом когда решетка оптимизирована. Наиболее частое использование FDS - модель LES. Например, при моделировании течении дыма через большие многокомнатные помещения невозможно решить горение и процессы переноса напрямую. Однако, для мелкомасштабных экспериментов с горением, возможно вычислить перенос и процессы горения напрямую.
Глава 5 содержит детальное описание численного решения уравнений импульса и давления. Для описания дальнейшей процедуры решения достаточно рассмотреть уравнение импульса, записанное как:
\begin{eqnarray} \frac{\partial\mathbf{u}}{\partial t}+\mathbf{F}+\nabla\mathcal{H}=0 \end{eqnarray} и уравнение давления как: \begin{eqnarray} \nabla^{2}\mathcal{H}=-\frac{\partial}{\partial t}(\nabla\cdot\mathbf{u})-\nabla\cdot\mathbf{F} \end{eqnarray} которое получено применением дивергенции к уравнению импульса.
 

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

Уравнение сохраниния энергии записано в терминах ЯВНОЙ ЭНТАЛЬПИИ, $h_{s}$ :
\begin{eqnarray} \frac{\partial}{\partial t}(\rho h_{s})+\nabla\cdot\rho h_{s}\mathbf{u}=\frac{\mathrm{D}p}{\mathrm{D}t}+\dot{q}^{'''}-\dot{q}_{b}^{'''}-\nabla\cdot\dot{\mathbf{q}}^{''}+\varepsilon \end{eqnarray} Явная энтальпия является функцией температуры:
\begin{eqnarray} h_{s}=\sum_{\alpha}Y_{\alpha}h_{s,\alpha};\qquad h_{s,\alpha}(T)=\int_{T_{0}}^{T}c_{p,\alpha}(T^{'})dT^{'} \end{eqnarray} Обратите внимание на использование индивидуальной производной, $D()/Dt=\partial()/\partial t+\mathbf{u}\cdot\nabla()$. Элемент $\dot{q}^{'''}$ - скорость выделения теплоты на единицу объема из химической реакции. Элемент $\dot{q}_{b}^{'''}$ - затраченная энергия на испарение капель. Элемент $\mathbf{\dot{q}}^{''}$ представляет проводящие и радиационные тепловые течения:
\begin{eqnarray} \mathbf{\dot{q}}^{''}=-k\nabla T-\sum_{\alpha}h_{s,\alpha}\rho D_{\alpha}\nabla Y_{\alpha}+\mathbf{\dot{q}}_{r}^{''} \end{eqnarray} где $k$ - теплопроводность.
 

Уравнение состояния

\begin{eqnarray} p=\frac{\rho\Re T}{\overline{W}} \end{eqnarray} Приблизительная форма уравнений Навье-Стокса подходит для низкоскоростных (малых чисел Маха, отношение скорости течения газового потока к скорости распространения звука) задач, используемых в модели. Приближение касается отфильтрованных звуковых волн с большим разбросом по температуре и плотности. Это дает уравнения элиптического вида, согласующиеся с низкоскоростными, температурно-конвективными процессами. На практике это значит, что давление по пространству $p(x,y,z)$ заменено средним или фоновым давлением $\overline{p}_{m}(z,t)$, которое только фунция времени и высоты над землей.
\begin{eqnarray} \overline{p}_{m}(z,t) & = & \rho T\Re\sum_{\alpha}Y_{\alpha}/W_{\alpha} \end{eqnarray} Беря индивидуальную производную от фонового давления и заменяя результат в уравнении сохранения энергии, получим выражение для дивергенции скорости, $\nabla\dot{\mathbf{u}}$, которая является важной величиной в численном алгоритме, потому что она эффективно убирает необходимость решать уравнение переноса для удельной энтальпии. Исходные члены уравнения сохранения энергии объединяются в дивергенции, которая появляется в уравнениях переноса массы. Температура находится из плотности и фонового давления через уравнение состояния.
 

Процедура решения

FDS использует конечно-разностное приближение второго порядка точности для главных уравнений на группе связанных прямоугольных сеток. Переменные потока обновляются во времени с использованием явной схемы Рунге-Кутта второго порядка. Эта глава описывает как этот алгоритм используется для вычисления во времени плотности, массовых долей веществ, скоростей компонентов, и фонового и возмущенного давления. Положим, что $\rho^{n},\, Y_{\alpha}^{n},\,\mathbf{u}^{n},\,\overline{p}_{m}^{n}\text{и }\mathcal{H}^{n}$ обозначают эти переменные в $n$-ый шаг времени.
1. Вычислим локально-среднее поле скорости $\mathbf{\bar{u}}^{n}$ (см. раздел 5.6.3 технического руководства к FDS)
2. Оценим $\rho,Y_{\alpha}\text{и}\overline{p}_{m}$ в следующий шаг времени с явным шагом Эйлера. Например, плотность оценивается как:
\begin{eqnarray} \frac{\rho^{*}-\rho^{n}}{\delta t}+\nabla\cdot\rho^{n}\bar{\mathbf{u}}^{n}=0 \end{eqnarray}

3. Обменяемся величинами $\rho^{*}\text{и}\, Y_{\alpha}^{*}$ на границах сетки.
4. Применим граничные условия для $\rho^{*}\text{и}\, Y_{\alpha}^{*}$.
5. Вычислим дивергенцию $\nabla\cdot\mathbf{\bar{u}}^{*}$, используя оцененные термодинамические величины. Важно, что на этом этапе поле скорости еще не оценено в следующий шаг времени, а только его дивергенция.
6. Решим уравнение Пуассона для флуктуации давления прямым решением на каждой индивидуальной сетке:
\begin{eqnarray} \nabla^{2}\mathcal{H}^{n}=-\left[\frac{\nabla\cdot\mathbf{u}^{*}-\nabla\cdot\bar{\mathbf{u}}^{n}}{\delta t}\right]-\nabla\cdot\bar{\mathbf{F}}^{n} \end{eqnarray} важно, что вектор $\mathbf{\bar{F}}^{n}=\mathbf{F}(\rho^{n},\,\bar{\mathbf{u}}^{n})$ вычислен с использованием локально-средних скоростей и эта дивергенция локально-среднего поля вычислена явно.
7. Оценим скорость в следующий шаг времени
\begin{eqnarray} \frac{\mathbf{u}^{*}-\bar{\mathbf{u}}^{n}}{\delta t}+\bar{\mathbf{F}}^{n}+\nabla\mathcal{H}^{n}=0 \end{eqnarray} важно, что дивергенция оцененного поля скоростей идентично равна оцененной дивергенции $\nabla\cdot\mathbf{u}^{*}$, которая была получена из оцененных термодинамических величин.
8. Проверим временной шаг в этой точке чтобы убедиться, что:
\begin{eqnarray} \delta t\, max\,\left(\frac{|u|}{\delta x},\frac{|v|}{\delta y},\frac{|w|}{\delta z}\right)<1;\qquad2\,\delta t\, v\,\left(\frac{1}{\delta x^{2}}+\frac{1}{\delta y^{2}}+\frac{1}{\delta z^{2}}\right)<1 \end{eqnarray} если временной шаг слишком большой, он уменьшается так, чтобы он удовлетворял обоим условиям, и процедура возвращается к началу временного шага. Если шаг по времени удовлетворяет критерию стабильности, процедура продолжает шаг корректора. (см. главу 5.5 ``Технического руководства к FDS'' для подробностей по стабильности).
Это завершает этап Предиктора на временном шаге. В этой точке величины $\mathcal{H}^{n}$ и компоненты $\mathbf{u}^{*}$ обмениваются на границах сетки через MPI-вызовы.
1. Вычислим локально-среднее поле скоростей $\bar{\mathbf{u}}^{*}$.
2. Применим вторую часть обновлений(вычислений) Рунге-Кутта для переменных массы. Например, плотность корректируется:
\begin{eqnarray} \frac{\rho^{n+1}-\frac{1}{2}(\rho^{n}+\rho^{*})}{\delta t/2}+\nabla\cdot\rho^{*}\mathbf{\bar{u}}^{*}=0 \end{eqnarray}

3. Обменяемся величинами $\rho^{n}$ и $Y_{\alpha}^{n}$ на границах сетки.
4. Применим граничные условия для $\rho^{n}$ и $Y_{\alpha}^{n}$.
5. Вычислим дивергенцию $\nabla\cdot\mathbf{u}^{n+1}$ из скорректированных термодинамических величин. Важно снова, что поле скорости еще не скорректировано на этом этапе.
6. Вычислим флуктуацию давления используя оцененные величины:
\begin{eqnarray} \nabla^{2}\mathcal{H}^{*}=-\left[\frac{\nabla\cdot\mathbf{u}^{n+1}-\frac{1}{2}(\nabla\cdot\bar{\mathbf{u}}^{*}+\nabla\cdot\bar{\mathbf{u}}^{n})}{\delta t/2}\right]-\nabla\cdot\bar{\mathbf{F}}^{*} \end{eqnarray}

7. Обновим скорость через вторую часть схемы Рунге-Кутта
\begin{eqnarray} \frac{u^{n+1}-\frac{1}{2}(\bar{\mathbf{u}}^{*}+\bar{\mathbf{u}}^{n})}{\delta t/2}+\bar{\mathbf{F}}^{*}+\nabla\mathcal{H}^{*}=0 \end{eqnarray} важно снова что дивергенция скорректированного поля скорости идентично равна дивергенции, вычисленной раньше.
8. В заключении временного шага величины $\mathcal{H}^{*}$ и $\mathbf{u}^{n+1}$ компоненты обмениваются на граниицах сетки через MPI-вызовы.
 

Пространственная дискретизация

Пространственные производные в главных уравнениях записаны как конечные разности второго порядка точности на прямоугольной сетке. Общий домен - это прямоугольный ящик, разделенный на прямоугольные ячейки. Каждой ячейке присваиваются индексы $i,\, j,\, k$, представляющие позицию ячейки в $x,\, y,\, z$ направлениях соответственно. Скалярные величины присваиваются к центру каждой ячейки; так $\rho_{ijk}^{n}$ - плотность в $n$-ый временной шаг в центре ячейки с индексами $i,\, j,\, k$. Векторные величины такие как скорость, присваиваются к соотвветствующим граням ячейки. Например, $u_{ijk}^{n}$ - $x$-компонента скорости в положительно-ориентированной грани $i,j,k$-ячейки; $u_{i-1,jk}^{n}$ - определена для отрицательно-ориентированной грани этой же ячейки.
 

Список литературы

[1] D.A. Anderson, J.C. Tannehill, and R.H. Pletcher. Computational Fluid Mechanics and Heat Transfer. Hemisphere Publishing Corporation, Philadelphia, Pennsylvania, 1984.