Теплофизика высоких температур, 2022, T. 60, № 5, стр. 752-762

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

И. С. Григорьев^1, 2, *

¹ АО “НПО ГИПО”
Казань, Россия

² КНИТУ–КАИ
Казань, Россия

^* E-mail: gipo@telebit.ru

Поступила в редакцию 24.05.2021
После доработки 07.04.2022
Принята к публикации 07.06.2022

EDN: MEDUNO
DOI: 10.31857/S0040364422050052

Полный текст (PDF)

Аннотация

Представлено описание статистических алгоритмов расчета рассеянного излучения по методу Монте-Карло. Акцент сделан на расчете светорассеяния в среде с распределенными объемными источниками излучения, для чего был выбран объект в виде осесимметричного объема газа, имитирующего струю самолетного двигателя с неравномерными распределениями температуры и концентрации рассеивающих частиц. Продемонстрированы теория и реализация на графическом процессоре трех алгоритмов Монте-Карло: прямого метода прослеживания траекторий фотонов и методов, основанных на разложении в ряд Неймана решений прямого и сопряженного уравнений переноса излучения (локальные оценки), с учетом весовых модификаций и различных индикатрис рассеяния. Наглядному сравнению работы алгоритмов способствует представление данных в виде функционала от интенсивности излучения и площади излучающей поверхности (силы излучения).

ВВЕДЕНИЕ

Процесс распространения света можно рассматривать как случайную марковскую цепь столкновений фотонов с веществом, которые приводят либо к рассеянию, либо к поглощению фотонов. Цепь Маркова характеризуется тем, что вероятность наступления каждого события зависит только от состояния, достигнутого в предыдущем событии. Метод Монте-Карло заключается в моделировании траекторий этой цепи на ЭВМ и вычислении статистической оценки для искомых функционалов. Под прямым методом расчета Монте-Карло переноса излучения понимается прослеживание траекторий частиц, испущенных из источника, и регистрация потока частиц на исследуемой поверхности. С использованием данного метода удобно определять полусферический поток излучения через регистрирующую площадку, а также угловое распределение потока, если разыгрывается достаточное число траекторий и размер регистрирующей площадки достаточно велик. Однако при переходе к точечному детектору данный метод становится непригодным, так как вероятность пересечения частицей соответствующей площадки стремится к нулю. Основа применения локальных оценок состоит в представлении некоторого функционала от плотности потока частиц и функции детектора в виде ряда Неймана, который сходится при наличии поглощения в среде или возможности вылета частицы за пределы расчетной области. При реализации этих расчетов также моделируется марковская траектория частицы, испытывающей поглощение и рассеяние, однако теперь на каждом акте рассеяния высчитывается вклад (локальная оценка), который частица вносит в регистрируемую интенсивность на детекторе, что позволяет проводить расчет с произвольным угловым разрешением. Методы на основе локальных оценок применялись в задачах расчета ядерных реакторов, а в задачах атмосферной оптики они впервые получили сильное развитие в работе [1]. С использованием данных методов реализовано большое количество кодов для расчета яркости атмосферы (с наименованием которых можно ознакомиться в проекте I3RC [2]). Также они широко применяются в задачах лазерной локации при решении уравнения переноса излучения с учетом поляризации [3–6].

Применение метода локальных оценок на основе решения сопряженного уравнения переноса излучения не так широко представлено в печати, как методов, основанных на решении прямого уравнения. Обстоятельством, упрощающим расчет, является рассмотрение вклада только точечных источников, что продемонстрировано в работах [7, 8], хотя в [7] представлены выражения, учитывающие излучение газов. Способ расчета объемных источников при описании излучения высокотемпературных струй с использованием локальных оценок представлен в работах [9, 10]. В [11] продемонстрировано применение метода прямого расчета Монте-Карло переноса излучения в осесимметричной струе. Особенностью данной работы является применение комбинированной локальной оценки при решении прямого уравнения переноса, а также использование весовых модификаций, заменяющих розыгрыш поглощения и вылет частицы из расчетной области.

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

ОПИСАНИЕ ОБЪЕКТА

Граница струи аппроксимирована параболой $\rho (x)$:

$\begin{gathered} \rho \left( x \right) = a{{x}^{2}} + bx + c, \\ y = \rho {\kern 1pt} \cos (\varphi ), \\ z = \rho {\kern 1pt} \sin (\varphi ). \\ \end{gathered} $

Используется сетка следующего типа (рис. 1): вдоль оси X разбиение равномерное, и оно также равномерно в каждом сечении $x = {\text{const}}$ с уровня $\rho = 0$ до границы $\rho (x)$.

Рис. 1.

Расчетная область.

Задается разбиение по углу $\varphi $, в результате чего поверхность струи оказывается покрыта площадками с площадью $dS = \rho d\rho d\varphi $ для плоских сторон и $d{{S}_{i}} = d\varphi \left( {a{\text{/}}3\left( {x_{{i + 1}}^{3} - x_{i}^{3}} \right) + } \right.b{\text{/2}}\left( {x_{{i + 1}}^{2} - x_{i}^{2}} \right) + $ $ + \,\,\left. {c\left( {{{x}_{{i + 1}}} - {{x}_{i}}} \right)} \right)$ в интервалах $[{{x}_{i}},{{x}_{{i + 1}}}]$ для боковой поверхности, образованной параболоидом вращения. Для регистрации покидающего струю излучения над каждой площадкой $d{{S}_{i}}$ строится локальная сферическая система координат, в которой зенитный угол изменяется в диапазоне $[0,\pi {\text{/}}2]$, а азимутальный в $[0,\pi ]$ (из-за наличия осевой симметрии). В расчетах использовалась сетка 120 × 120 для дискретизации сферических систем координат и 19 × 19 ячеек разбиения по осям X и Y для дискретизации расчетного объема. Ввиду осевой симметрии распределения температур и концентраций примеси излучение от струи также принято осесимметричным. Поэтому запись результатов происходит только на полосе $\varphi = 0$. При использовании методов локальных оценок подразумевается наличие точечного датчика в середине каждой площадки $d{{S}_{i}}$, а в методе прямого расчета запись в локальную систему координат над $d{{S}_{i}}$ осуществляется со всей осесимметричной поверхности (в интервале $[{{x}_{i}},{{x}_{{i + 1}}}]$) с весом ${{d\varphi } \mathord{\left/ {\vphantom {{d\varphi } {2\pi }}} \right. \kern-0em} {2\pi }}$. В качестве рассеивающих частиц использовались псевдочастицы с заданной индикатрисой рассеяния, концентрация которых распределена по объему пропорционально распределению температуры при значении на срезе сопла ${{10}^{8}}$ м^–3. Предполагается, что дополнительного поглощения данные частицы не вносят. Изменение коэффициента рассеяния производилось путем изменения сечения рассеяния частиц в диапазоне от $0$ до $80 \times {{10}^{{ - 10}}}$ м².

ИСПОЛЬЗУЕМЫЕ МОДЕЛИ И ПРИБЛИЖЕНИЯ

Рассмотрим выражение для интенсивности излучения $I({\mathbf{r}},{\mathbf{\omega }})$ в точке ${\mathbf{r}}$ в направлении ${\mathbf{\omega }}$ в поглощающей и рассеивающей среде $\mathbb{R}$ с заданными оптическими характеристиками: коэффициентом ослабления ${{\sigma }_{e}}({\mathbf{r}})$, альбедо однократного рассеяния $\varpi ({\mathbf{r}})$, нормированной фазовой функцией (индикатрисы) рассеяния $P({\mathbf{r}},{\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }})$ и источниками теплового излучения $B[T({\mathbf{r}})]$ (описываемого функцией Планка) с температурой $T({\mathbf{r}})$ в состоянии локального термодинамического равновесия:

$\omega \cdot \nabla I\left( {r,\omega } \right) = {{\sigma }_{e}}\left( r \right)\left( { - I\left( {r,\omega } \right) + S\left( {r,\omega } \right)} \right),\,\,\,\,{\mathbf{r}} \in \mathbb{R},$

где функция источников излучения представлена в виде

$\begin{gathered} S({\mathbf{r}},{\mathbf{\omega }}) = \varpi ({\mathbf{r}})\int\limits_{4\pi } {P(} {\mathbf{r}},\omega {\kern 1pt} ' \to {\mathbf{\omega }})I({\mathbf{r}},\omega {\kern 1pt} ')d\omega {\kern 1pt} '\,\, + \\ + \,\,\left( {1 - \varpi ({\mathbf{r}})} \right)B[T({\mathbf{r}})]. \\ \end{gathered} $

На границе расчетной области $\partial \mathbb{R}$ в точке ${{r}_{b}} \in \partial \mathbb{R}$ могут присутствовать источник внешнего излучения (например, солнечное излучение) ${{S}_{0}}({{{\mathbf{r}}}_{{\mathbf{b}}}},{\mathbf{\omega }})$ или твердая граница, характеризуемая направляющим вектором ${\mathbf{n}}({{{\mathbf{r}}}_{{\mathbf{b}}}})$, альбедо отражения плоского слоя $a({{{\mathbf{r}}}_{{\mathbf{b}}}},{\mathbf{\omega }})$, нормированной фазовой функцией отражения ${{P}_{s}}({{{\mathbf{r}}}_{{\mathbf{b}}}},\omega {\kern 1pt} ' \to {\mathbf{\omega }})$, источником теплового излучения $B[T({{{\mathbf{r}}}_{{\mathbf{b}}}})]$ с температурой $T({{{\mathbf{r}}}_{{\mathbf{b}}}})$ и направленной излучательной способностью $\varepsilon ({{{\mathbf{r}}}_{{\mathbf{b}}}},{\mathbf{\omega }})$. Интенсивность излучения на границе $\partial \mathbb{R}$ может быть записана в виде

$\begin{gathered} I\left( {{{r}_{b}},\omega } \right) = \\ = \mathop \smallint \limits_{\left\langle {n\left( {{{r}_{b}}} \right),\omega {\kern 1pt} '} \right\rangle \, < \,0} a({{r}_{b}},\omega {\kern 1pt} '){{P}_{s}}\left( {{{r}_{b}},\omega {\kern 1pt} ' \to \omega } \right)I\left( {{{r}_{b}},\omega {\kern 1pt} '} \right)d{{\omega }}{\kern 1pt} {\text{'}} + \\ + \,\,\varepsilon \left( {{{r}_{b}},\omega } \right)B\left[ {T\left( {{{r}_{b}}} \right)} \right] + {{S}_{0}}\left( {{{r}_{b}},\omega } \right), \\ {{{\mathbf{r}}}_{b}} \in \partial \mathbb{R},\,\,\,\,\left\langle {n\left( {{{r}_{b}}} \right),\omega } \right\rangle > 0. \\ \end{gathered} $

Силу излучения $F\left( {{{r}_{{{\text{obs}}}}}} \right)$, регистрируемую наблюдателем из положения ${{r}_{{{\text{obs}}}}}$ от излучающего объекта площадью $\Sigma $, можно записать в виде

(1)

$\begin{gathered} F\left( {{{r}_{{{\text{obs}}}}}} \right) = \mathop \sum \limits_{\Sigma ,\left\langle {{{n}_{0}},{{n}_{b}}} \right\rangle \, > \,0} I\left( {{{r}_{b}}\left( {d{{S}_{i}}} \right),\,\,{{n}_{b}}\left( {d{{S}_{i}}} \right)} \right) \times \\ \times \,\,\left\langle {{{n}_{0}},{{n}_{b}}} \right\rangle d{{S}_{i}},~\,\,\,\,{{n}_{0}} = \frac{{{{r}_{{{\text{obs}}}}} - {{r}_{b}}}}{{\left| {{{r}_{{{\text{obs}}}}} - {{r}_{b}}} \right|}}, \\ \end{gathered} $

где ${{{\mathbf{n}}}_{{\mathbf{b}}}}$ – вектор внешней нормали к границе струи. На практике данная величина используется для оценки облученности зрачка оптико-электронной системы. Выражение для интенсивности излучения в операторной форме имеет вид

(2)

$I\left( {r,\omega } \right) = \mathcal{K}I\left( {r,\omega } \right) + {{I}_{0}}\left( {r,\omega } \right),$

где интегральный оператор переноса излучения $\mathcal{K}$ определен как

$\tau ({\mathbf{r}}{\kern 1pt} ',{\mathbf{r}})$ – оптическая толщина среды между точками ${\mathbf{r}}$ и ${\mathbf{r}}{\kern 1pt} '$. Вклад нерассеянного (прямого) излучения от объемных источников и внешней подсветки ${{I}_{0}}({\mathbf{r}},{\mathbf{\omega }})$ имеет вид

(3)

$\begin{gathered} {{I}_{0}}\left( {r,\omega } \right) = \left( {\varepsilon \left( {{{r}_{b}},\omega } \right)B\left[ {T\left( {{{r}_{b}}} \right)} \right] + {{S}_{0}}\left( {{{r}_{b}},\omega } \right)} \right){{{\text{e}}}^{{ - \tau \left( {{{r}_{b}},r} \right)}}} + \\ \underbrace { + \,\,\mathop \smallint \limits_0^{{{r}_{b}}} {{\sigma }_{e}}\left( {r{\kern 1pt} '} \right){{{\text{e}}}^{{ - \tau \left( {r{\kern 1pt} ',r} \right)}}}\left( {1 - \varpi \left( {r{\kern 1pt} '} \right)} \right)B\left[ {T\left( {r{\kern 1pt} '} \right)} \right]dr{\kern 1pt} '}_{I_{0}^{V}\left( {r,\omega } \right)~\,\, - ~\,\,{\text{объемные\;источники}}}. \\ \end{gathered} $

Стоит отметить, что на однородном (${{\sigma }_{e}}({\mathbf{r}}),T({\mathbf{r}}) = $ = const) участке пути длиной $s$ выражение для нерассеянного излучения объемных источников в (3) имеет аналитическое решение

$I_{0}^{V}({\mathbf{r}},s) = (1 - \varpi ({\mathbf{r}}))B[T({\mathbf{r}})]\left( {1 - \exp ( - {{\sigma }_{e}}({\mathbf{r}})s)} \right).$

Решение уравнения (2) представляется в виде ряда Неймана

$I = \mathop \sum \limits_{n = 0}^\infty {{I}_{n}},\,\,\,\,{{I}_{n}} = {{\mathcal{K}}^{n}}{{I}_{0}} = \mathcal{K}{{I}_{{n - 1}}}.$

Далее индекс $n$ обозначает кратность рассеивания. Функция плотности столкновений частиц ${{\Phi }_{n}}({\mathbf{r}},{\mathbf{\omega }})$ представляет собой произведение коэффициента ослабления ${{\sigma }_{e}}({\mathbf{r}})$ и интенсивности ${{I}_{n}}({\mathbf{r}},{\mathbf{\omega }})$, а функция плотности источников излучения $\psi ({\mathbf{r}},{\mathbf{\omega }})$ имеет вид

$\psi \left( {r,\omega } \right) = {{\sigma }_{e}}\left( r \right)\left( {1 - \varpi \left( r \right)} \right)B\left[ {T\left( r \right)} \right],\,\,\,\,r \in \mathbb{R}.{\text{\;}}$

Функция плотности столкновений $n$-го порядка записывается как

$\begin{gathered} {{\Phi }_{n}}({\mathbf{r}},{\mathbf{\omega }}) = TC{{\Phi }_{{n - 1}}}({\mathbf{r}},{\mathbf{\omega }}), \\ {{\Phi }_{1}}({\mathbf{r}},{\mathbf{\omega }}) = T\psi ({\mathbf{r}},{\mathbf{\omega }}), \\ \end{gathered} $

где оператор рассеяния частицы

$Cf({\mathbf{r}},{\mathbf{\omega }}) = \int\limits_{4\pi } {\varpi ({\mathbf{r}})P({\mathbf{r}},{\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }})f} ({\mathbf{r}},{\mathbf{\omega }})d\omega {\kern 1pt} '$

представляет собой произведение вероятности частицы выжить (не быть поглощенной средой) $\varpi ({\mathbf{r}})$ и вероятности рассеяться из направления ${\mathbf{\omega }}{\kern 1pt} '$ в направление ${\mathbf{\omega }}$, а оператор переноса частицы

$\begin{gathered} Tf\left( {r,\omega } \right) = \mathop \smallint \limits_\mathbb{R} {{\sigma }_{e}}\left( r \right)\frac{{\exp \left( { - \tau \left( {r{\kern 1pt} ',r} \right)} \right)}}{{{{{\left| {r - r{\kern 1pt} '} \right|}}^{2}}}}\delta \times \\ \times \,\,\left( {\omega - \frac{{\left( {r - r{\kern 1pt} '} \right)}}{{\left| {r - r{\kern 1pt} '} \right|}}} \right)f\left( {r{\kern 1pt} ',\omega } \right)dr{\kern 1pt} ', \\ \end{gathered} $

в котором интегрирование происходит по расчетной области, обозначает перенос частицы вдоль вектора ${\mathbf{\omega }}$ на длину свободного пробега, имеющую плотность распределения

(4)

$p({{{\mathbf{r}}}_{{\mathbf{n}}}}) = {{\sigma }_{e}}({{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}})\exp \left( { - \tau ({{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}},{{{\mathbf{r}}}_{{\mathbf{n}}}})} \right)$

из ($n - 1$)-й точки ${\mathbf{r}}{\kern 1pt} '$ в $n$-ю точку столкновения ${\mathbf{r}}$. Сопряженные операторы рассеяния и переноса имеют вид

$\begin{gathered} {{C}^{\dag }}f\left( {r,\omega } \right) = \mathop \smallint \limits_{4\pi } \varpi \left( r \right)P\left( {r,\omega \to \omega {\kern 1pt} '} \right)f\left( {r,\omega {\kern 1pt} '} \right)d\omega {\kern 1pt} ', \\ {{T}^{\dag }}f\left( {r,\omega } \right) = \mathop \smallint \limits_\mathbb{R} {{\sigma }_{e}}(r{\kern 1pt} ')\frac{{\exp \left( { - \tau \left( {r{\kern 1pt} ',r} \right)} \right)}}{{{{{\left| {r - r{\kern 1pt} '} \right|}}^{2}}}} \times \\ \times \,\,\delta \left( {\omega + \frac{{\left( {r - r{\kern 1pt} '} \right)}}{{\left| {r - r{\kern 1pt} '} \right|}}} \right)f\left( {r{\kern 1pt} ',\omega } \right)dr{\kern 1pt} '. \\ \end{gathered} $

ЛОКАЛЬНЫЕ ОЦЕНКИ ПРИ РЕШЕНИИ ПРЯМОГО УРАВНЕНИЯ ПЕРЕНОСА

Выражения для локальных оценок при решении прямого и сопряженного уравнений переноса излучения представлены в работе [7]. Интенсивность излучения, регистрируемая приемником с функцией отклика $\varphi ({\mathbf{r}},{\mathbf{\omega }})$, может быть выражена через сумму вкладов функционала от плотности столкновений

(5)

$\begin{gathered} {{I}_{\varphi }} = \sum\limits_{n = 1}^\infty {({{\Phi }_{{n + 1}}},{\varphi \mathord{\left/ {\vphantom {\varphi {{{\sigma }_{e}}}}} \right. \kern-0em} {{{\sigma }_{e}}}})} = \\ = \sum\limits_{n = 1}^\infty {\left( {{{{\left( {TC} \right)}}^{{n - 1}}}T\psi ,{{C}^{\dag }}{{T}^{\dag }}{\varphi \mathord{\left/ {\vphantom {\varphi {{{\sigma }_{e}}}}} \right. \kern-0em} {{{\sigma }_{e}}}}} \right)} . \\ \end{gathered} $

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

(6)

$\Psi _{n}^{F} = C_{{}}^{\dag }T_{{}}^{\dag }{\varphi \mathord{\left/ {\vphantom {\varphi {{{\sigma }_{e}}}}} \right. \kern-0em} {{{\sigma }_{e}}}}{\text{.}}$

Запись (6) наглядно демонстрирует принцип реализации прямого уравнения с помощью локальных оценок: в левой части функционала описана марковская цепь траектории частицы, стартующей из источника излучения, над которой последовательно разыгрывается сначала транспортное ядро, а затем ядро рассеяния. При этом перед розыгрышем ядра рассеяния вычисляется вклад (6), стоящий в правой части функционала (5). Рассматривается точечный детектор в точке ${\mathbf{r}}{\kern 1pt} *$, который фиксирует интенсивность излучения в направлении $\Delta \Omega $. Соответствующая функция отклика детектора имеет вид

(7)

$\begin{gathered} {{\varphi }_{{\Delta \Omega }}}\left( {r,\omega } \right) = \delta \left( {r - r{\kern 1pt} *} \right){{\chi }_{{\Delta \Omega }}}\left( \omega \right), \\ {{\chi }_{{\Delta \Omega }}}\left( \omega \right) = \left\{ {\begin{array}{*{20}{c}} {\frac{1}{{\Delta \Omega }}}&{\omega \in \Delta \Omega ,} \\ 0&{\omega \notin \Delta \Omega .} \end{array}} \right.{\text{\;\;}} \\ \end{gathered} $

Выпишем явный вид (6) с учетом (7):

$\begin{gathered} \Psi _{n}^{{1,F}}({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}}) = \varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})P({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}} \to {\mathbf{\omega }}{\kern 1pt} *) \times \\ \times \,\,{{\exp \left( { - \tau ({{{\mathbf{r}}}_{{\mathbf{n}}}},{\mathbf{r}}{\kern 1pt} *)} \right)} \mathord{\left/ {\vphantom {{\exp \left( { - \tau ({{{\mathbf{r}}}_{{\mathbf{n}}}},{\mathbf{r}}{\kern 1pt} *)} \right)} {{{{\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}^{2}}}}} \right. \kern-0em} {{{{\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}^{2}}}}{{\chi }_{{\Delta \Omega }}}\left[ {{{\left( {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right)} \mathord{\left/ {\vphantom {{\left( {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right)} {\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}} \right. \kern-0em} {\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}} \right]. \\ \end{gathered} $

Такая форма вклада в литературе носит название простой локальной оценки [1], согласно которой вклад в направление $\Delta \Omega $ вносится только тогда, когда положение частицы при рассеянии входит в данный телесный угол. При уменьшении угла $\Delta \Omega $ дисперсия такой оценки возрастает. Наряду с простой локальной оценкой существует двойная локальная оценка, строящаяся с использованием дополнительной промежуточной точки ${\mathbf{r}}{\kern 1pt} ' \in {\mathbf{\Delta \Omega }}$, распределенной в ${\mathbf{\Delta \Omega }}$ с плотностью (4), и имеющая вид

(8)

$\begin{gathered} \Psi _{n}^{{2,F}}({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}}) = \varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})\varpi ({\mathbf{r}}{\kern 1pt} ')\frac{{\exp \left( { - \tau ({{{\mathbf{r}}}_{{\mathbf{n}}}},{\mathbf{r}}{\kern 1pt} ')} \right)}}{{{{{\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} '} \right|}}^{2}}}} \times \\ \times \,\,P({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}} \to {\mathbf{\omega }}{\kern 1pt} ')P({\mathbf{r}}{\kern 1pt} ',{\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }}{\kern 1pt} *), \\ \end{gathered} $

где ${\mathbf{\omega }}{\kern 1pt} ' = {{({\mathbf{r}}{\kern 1pt} '\,\, - {{{\mathbf{r}}}_{{\mathbf{n}}}})} \mathord{\left/ {\vphantom {{({\mathbf{r}}{\kern 1pt} '\,\, - {{{\mathbf{r}}}_{{\mathbf{n}}}})} {\left| {{\mathbf{r}}{\kern 1pt} '\, - {{{\mathbf{r}}}_{{\mathbf{n}}}}} \right|,}}} \right. \kern-0em} {\left| {{\mathbf{r}}{\kern 1pt} '\, - {{{\mathbf{r}}}_{{\mathbf{n}}}}} \right|,}}\,\,\,\,{\mathbf{\omega }}{\kern 1pt} * = {\text{(}}{{{\mathbf{r}}{\kern 1pt} * - \,\,{\mathbf{r}}{\kern 1pt} ')} \mathord{\left/ {\vphantom {{{\mathbf{r}}{\kern 1pt} * - \,\,{\mathbf{r}}{\kern 1pt} ')} {\left| {{\mathbf{r}}{\kern 1pt} * - \,\,{\mathbf{r}}{\kern 1pt} '} \right|}}} \right. \kern-0em} {\left| {{\mathbf{r}}{\kern 1pt} * - \,\,{\mathbf{r}}{\kern 1pt} '} \right|}}$. Вклад двойной локальной оценки строится для каждого текущего положения $({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}})$. Простую и двойную локальные оценки можно комбинировать, как описано в [5], следующим способом: если ${{{\mathbf{r}}}_{{\mathbf{n}}}} \notin {\mathbf{\Delta \Omega }}$, то записывается двойная локальная оценка, в противном случае есть два варианта:

1) если в предыдущей точке ${{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}}$ вычислялась двойная локальная оценка, то вклада нет;

2) если в ${{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}}$ вклада нет или в ${{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}}$ вычислялась простая локальная оценка, то записывается простая локальная оценка.

Комбинированную оценку для данного метода обозначим просто $\Psi _{n}^{F}$. Выражение для регистрируемой интенсивности в ${\mathbf{\Delta \Omega }}$ имеет вид

(9)

${{I}_{{{\mathbf{\Delta \Omega }}}}} = PM\left( {\sum\limits_{n = 1}^\infty {\Psi _{n}^{F}} } \right),$

(10)

$P = \sum\limits_q {{{P}_{q}}} ,\,\,\,\,{{P}_{q}} = 4\pi {{V}_{q}}{{\sigma }_{{a,q}}}B({{T}_{q}}),$

где ${{V}_{q}},{{\sigma }_{{a,q}}},{{P}_{q}}$ – объем ячейки, коэффициент поглощения, мощность излучения ячейки $q$ (в Вт/мкм); $M$ – математическое ожидание. При построении промежуточной точки ${\mathbf{r}}{\kern 1pt} ' \in {\mathbf{\Delta \Omega }}$ она распределяется с плотностью (4). Однако, если оптическая длина пути вдоль направления ${\mathbf{\Delta \Omega }}$ мала, то с большой долей вероятности разыгрываемые точки окажутся за пределами расчетной области, что обнулит вклад от двойной локальной оценки. Чтобы избежать этого, воспользуемся модифицированной плотностью, учитывающей вылет частицы из расчетной области:

(11)

$\begin{gathered} \tilde {p}({\mathbf{r}}{\kern 1pt} ') = {{{{\sigma }_{e}}({\mathbf{r}}{\kern 1pt} ')\exp \left( { - \tau ({\mathbf{r}}{\kern 1pt} *,{\mathbf{r}}{\kern 1pt} ')} \right)} \mathord{\left/ {\vphantom {{{{\sigma }_{e}}({\mathbf{r}}{\kern 1pt} ')\exp \left( { - \tau ({\mathbf{r}}{\kern 1pt} *,{\mathbf{r}}{\kern 1pt} ')} \right)} {\left( {1 - \exp \left( { - \tau ({\mathbf{r}}{\kern 1pt} *,{{{\mathbf{r}}}_{{\mathbf{b}}}})} \right)} \right)}}} \right. \kern-0em} {\left( {1 - \exp \left( { - \tau ({\mathbf{r}}{\kern 1pt} *,{{{\mathbf{r}}}_{{\mathbf{b}}}})} \right)} \right)}}, \\ {{\left( {{\mathbf{r}}{\kern 1pt} '\,\, - {\mathbf{r}}{\kern 1pt} *} \right)} \mathord{\left/ {\vphantom {{\left( {{\mathbf{r}}{\kern 1pt} '\,\, - {\mathbf{r}}{\kern 1pt} *} \right)} {\left| {{\mathbf{r}}{\kern 1pt} '\,\, - {\mathbf{r}}{\kern 1pt} *} \right|}}} \right. \kern-0em} {\left| {{\mathbf{r}}{\kern 1pt} '\,\, - {\mathbf{r}}{\kern 1pt} *} \right|}} = {{\left( {{{{\mathbf{r}}}_{{\mathbf{b}}}} - {\mathbf{r}}{\kern 1pt} *} \right)} \mathord{\left/ {\vphantom {{\left( {{{{\mathbf{r}}}_{{\mathbf{b}}}} - {\mathbf{r}}{\kern 1pt} *} \right)} {\left| {{{{\mathbf{r}}}_{{\mathbf{b}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}} \right. \kern-0em} {\left| {{{{\mathbf{r}}}_{{\mathbf{b}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}. \\ \end{gathered} $

При использовании модифицированной плотности распределения промежуточная точка ${\mathbf{r}}{\kern 1pt} '$ всегда принадлежит расчетной области, но в таком случае у нее появится компенсирующий вес $\left( {1 - \exp \left[ { - \tau ({{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}},{{{\mathbf{r}}}_{{\mathbf{b}}}})} \right]} \right)$, который также должен учитываться в выражении для двойной локальной оценки (8). Положение точки ${\mathbf{r}}{\kern 1pt} '$ разыгрывалось на прямой, проходящей из ${\mathbf{r}}{\kern 1pt} *$ по середине ${\mathbf{\Delta \Omega }}$. Плотность (11) можно использовать и при моделировании положений $({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}})$, что увеличит число столкновений, так как частица теперь не покинет расчетную область. Ограничением на расчет в таком случае является достижение весом частицы некоторого порогового значения (обычно ${{10}^{{ - 5}}}$).

При моделировании марковской траектории частицы $({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}})$ длина свободного пробега распределена с плотностью (4) или (11). В очередной точке столкновения при этом должны разыгрываться два типа взаимодействия: рассеяние частицы с вероятностью $\varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$ или поглощение частицы с вероятностью $1 - \varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$. Для такого розыгрыша обычно используется метод русской рулетки, при котором случайное число $\xi $, равномерно распределенное в интервале (0, 1), сравнивается с $\varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$. Если $\xi < \varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$, то разыгрывается рассеяние, иначе – поглощение. Однако наличие $\varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$ в выражениях для простой и двойной локальных оценок показывает, что такой розыгрыш явно не происходит, а во вкладах учитывается сразу все рассеянное излучение. Такой весовой метод называется статистическим весом по поглощению, и с его использованием поглощение не разыгрывается, но вес частицы умножается на альбедо. Данную весовую модификацию можно применять и в методе прямого расчета.

Порядок проведения расчета по прямому уравнению следующий. Расчетная область делится на подобласти с постоянными оптическими характеристиками. Ввиду неравномерности используемых распределений температур и концентраций, а также различных величин объемов подобластей, заключенные в них мощности излучения различны. Выбирается ячейка $q$ с минимальным по области значением испускаемой мощности ${{P}_{{\min }}}$, и задается число ${{N}_{{\min }}}$ частиц с энергией ${{{{P}_{{\min }}}} \mathord{\left/ {\vphantom {{{{P}_{{\min }}}} {{{N}_{{\min }}}}}} \right. \kern-0em} {{{N}_{{\min }}}}}$, которые испускаются из ячейки $q$ изотропно и из равномерно распределенной в ${{V}_{q}}$ точки пространства. Для остальных ячеек число разыгрываемых частиц рассчитывается как ${{N}_{q}} = {{N}_{{\min }}}{{{{E}_{q}}} \mathord{\left/ {\vphantom {{{{E}_{q}}} {{{E}_{{\min }}}}}} \right. \kern-0em} {{{E}_{{\min }}}}}$. Чем большее значение ${{N}_{{\min }}}$ используется, тем статистически точнее результат.

Рассмотрим процедуру вычисления вклада рассеянного излучения при прослеживании траектории одной частицы. Задаются начальный вес $Q = 1$ и положение промежуточной точки ${\mathbf{r}}{\kern 1pt} '$ по плотности (11) для детектируемого направления $\Delta \Omega $. В рамках прослеживания траектории частицы положение ${\mathbf{r}}{\kern 1pt} '$ не меняется, но в то же время оно различно для разных частиц. Впрочем положение ${\mathbf{r}}{\kern 1pt} '$ можно менять при каждом столкновении, но это не приводит к видимым изменениям. Далее последовательно производятся следующие операции:

1) розыгрыш длины свободного пробега по плотности (4) или (11). Если использовалась плотность (11), то вес ${{Q}_{n}}$ умножается на $\left( {1 - \exp \left[ { - \tau ({{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}},{{{\mathbf{r}}}_{{\mathbf{b}}}})} \right]} \right)$;

2) вычисление вклада $\Psi _{n}^{F}$ с использованием комбинации простой и двойной локальных оценок;

3) изменение веса частицы по правилу Q_n =$ = {{Q}_{n}}\varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$, розыгрыш нового направления движения.

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

$\begin{gathered} \Psi _{n}^{{1,F}}({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}}) = {{Q}_{n}}P({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}} \to {\mathbf{\omega }}{\kern 1pt} *){{\chi }_{{\Delta \Omega }}} \times \\ \times \,\,\left( {{{\left( {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right)} \mathord{\left/ {\vphantom {{\left( {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right)} {\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}} \right. \kern-0em} {\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}} \right){{\exp \left( { - \tau ({{{\mathbf{r}}}_{{\mathbf{n}}}},{\mathbf{r}}{\kern 1pt} *)} \right)} \mathord{\left/ {\vphantom {{\exp \left( { - \tau ({{{\mathbf{r}}}_{{\mathbf{n}}}},{\mathbf{r}}{\kern 1pt} *)} \right)} {{{{\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}^{2}}}}} \right. \kern-0em} {{{{\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} *} \right|}}^{2}}}}, \\ \Psi _{n}^{{2,F}}({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}}) = {{Q}_{n}}\left( {1 - \exp \left[ { - \tau ({\mathbf{r}}{\kern 1pt} *,{{{\mathbf{r}}}_{{\mathbf{b}}}})} \right]} \right) \times \\ \times \,\,\varpi ({\mathbf{r}}{\kern 1pt} ')P({{{\mathbf{r}}}_{{\mathbf{n}}}},{{{\mathbf{\omega }}}_{{\mathbf{n}}}} \to {\mathbf{\omega }}{\kern 1pt} ')P({\mathbf{r}}{\kern 1pt} ',{\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }}{\kern 1pt} *)\frac{{\exp \left( { - \tau ({{{\mathbf{r}}}_{{\mathbf{n}}}},{\mathbf{r}}{\kern 1pt} ')} \right)}}{{{{{\left| {{{{\mathbf{r}}}_{{\mathbf{n}}}} - {\mathbf{r}}{\kern 1pt} '} \right|}}^{2}}}}. \\ \end{gathered} $

Таким образом, для каждого направления ${\mathbf{\Delta \Omega }}$ каждой точки на поверхности ${\mathbf{r}}{\kern 1pt} *$ рассчитывается интенсивность (9), по которой затем вычисляется сила излучения (1). В варианте алгоритма прямого расчета Монте-Карло разыгрывается таким же образом распределенное по объему число частиц, но используется только вес по поглощению. В данном методе решение нельзя разделить на прямую и рассеянную компоненты, поэтому сила рассчитывается только для полного излучения. Регистрация вылетающих частиц происходит со всей поверхности, при этом запись осуществляется в локальные системы координат на площадках $d{{S}_{i}}$ при угле $\varphi = 0$ с весом $d\varphi {\text{/}}2\pi $. Получающиеся величины умножаются на энергию одного пакета ${{{{P}_{{\min }}}} \mathord{\left/ {\vphantom {{{{P}_{{\min }}}} {{{N}_{{\min }}}}}} \right. \kern-0em} {{{N}_{{\min }}}}}$ и после деления на величину телесного угла детектируемого направления представляют собой силу излучения.

ЛОКАЛЬНЫЕ ОЦЕНКИ ПРИ РЕШЕНИИ СОПРЯЖЕННОГО УРАВНЕНИЯ ПЕРЕНОСА

При моделировании по данному методу траектория частицы начинается в детекторе. Решение для интенсивности излучения в этом случае представляется в виде

(12)

${{I}_{\varphi }} = \sum\limits_{n = 1}^\infty {\left( {{{{\left( {{{T}^{B}}{{C}^{B}}} \right)}}^{{n - 1}}}{{T}^{B}}\varphi ,C{\text{/}}{{\sigma }_{e}}T\psi } \right)} ,$

где $C_{{}}^{B} = C_{{}}^{\dag }$, а оператор $T_{n}^{B}$ определяется выражением

$\begin{gathered} T_{n}^{B}f\left( {r,\omega } \right) = \mathop \smallint \limits_\mathbb{R} {{\sigma }_{e}}\left( r \right)\frac{{\exp \left( { - \tau \left( {r{\kern 1pt} ',r} \right)} \right)}}{{{{{\left| {r - r{\kern 1pt} '} \right|}}^{2}}}} \times \\ \times \,\,\delta \left( {\omega + \frac{{\left( {r - r{\kern 1pt} '} \right)}}{{\left| {r - r{\kern 1pt} '} \right|}}} \right)f\left( {r{\kern 1pt} ',\omega } \right)dr{\kern 1pt} '. \\ \end{gathered} $

Видно, что оператор $T_{n}^{B}$ аналогичен оператору $T$ с заменой направления движения на противоположное. Согласно (12), производится моделирование траекторий, начинающихся в детекторе, с использованием сопряженных операторов переноса и рассеяния, и на каждом шаге после переноса вычисляется вклад

(13)

$\Psi _{n}^{B} = C{\text{/}}{{\sigma }_{e}}T\psi .$

Так как тепловые источники распределены по всей расчетной области, оценка (13) подразумевает интегрирование по всему пространству: в текущей точке столкновения $({\mathbf{r}},{\mathbf{\omega }})$ производится интегрирование по полному телесному углу приходящего нерассеянного излучения со всех направлений с весами $P({\mathbf{r}},{\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }})$. В это интегрирование входят как распределенные тепловые источники, так и внешняя подсветка. Расчетное выражение имеет вид

(14)

$\Psi _{n}^{B}({\mathbf{r}},{\mathbf{\omega }}) = {{Q}_{n}}\int\limits_{4\pi } {{{I}_{0}}(} {\mathbf{r}},{\mathbf{\omega }}{\kern 1pt} ')P({\mathbf{r}},{\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }})d\omega {\kern 1pt} '.$

Также необходимо соблюдать нормировку $\int_{4\pi } {P({\mathbf{r}},{\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }}} )d\omega {\kern 1pt} ' = 1$, что налагает требования на степень дискретизации полного телесного угла при использовании различных индикатрис рассеяния. Выражение для регистрируемой интенсивности в ${\mathbf{\Delta \Omega }}$ выглядит так

(15)

${{I}_{{{\mathbf{\Delta \Omega }}}}} = M\left( {\sum\limits_{n = 1}^\infty {\Psi _{n}^{B}} } \right).$

В начале моделирования траектории из положения детектора в определенном направлении $({\mathbf{r}}{\kern 1pt} *,{\mathbf{\Delta \Omega }})$ задается вес $Q = 1$, а затем последовательно:

1) разыгрывается длина свободного пробега по плотности (4) или (11). Если использовалась плотность (11), то вес умножается на $\left( {1 - \exp \left[ { - \tau ({{{\mathbf{r}}}_{{{\mathbf{n}} - {\mathbf{1}}}}},{{{\mathbf{r}}}_{{\mathbf{b}}}})} \right]} \right)$;

2) рассчитывается вклад $\Psi _{n}^{B}$;

3) изменяется вес частицы по правилу ${{Q}_{n}} = {{Q}_{n}}\varpi ({{{\mathbf{r}}}_{{\mathbf{n}}}})$, разыгрывается новое направление движения.

РОЗЫГРЫШ ДЛИНЫ ПРОБЕГА И УГЛА РАССЕЯНИЯ

В однородной среде плотность распределения длины свободного пробега $p(l)$ (4) связана с функцией распределения $F(l)$ свободного пробега соотношением [12]

$p(l) = {{dF(l)} \mathord{\left/ {\vphantom {{dF(l)} {dl}}} \right. \kern-0em} {dl}} = \left[ {1 - F(l)} \right]{{\sigma }_{e}}.$

Решение этого уравнения относительно $F(l)$ имеет вид

$F(l) = 1 - \exp \left( { - {{\sigma }_{e}}l} \right).$

Функции распределения $F(l)$ соответствует равномерно распределенная на $(0,1)$ случайная величина $\gamma $, что приводит к решению относительно $l$ в виде

$l = - \sigma _{e}^{{ - 1}}{\kern 1pt} \ln (1 - \gamma ) = - \sigma _{e}^{{ - 1}}{\kern 1pt} \ln (\gamma ).$

В неоднородной среде моделируется не длина пробега, а оптическая толщина, по значению которой затем находится соответствующая длина перемещения. Для модифицированной плотности (11) функция распределения имеет вид

$F(l) = {{\left( {1 - \exp \left[ { - \tau (l)} \right]} \right)} \mathord{\left/ {\vphantom {{\left( {1 - \exp \left[ { - \tau (l)} \right]} \right)} {\left( {1 - \exp \left[ { - \tau ({{r}_{b}})} \right]} \right)}}} \right. \kern-0em} {\left( {1 - \exp \left[ { - \tau ({{r}_{b}})} \right]} \right)}},$

что приводит к выражению для моделируемой оптической толщины

$\tau = - \ln \left[ {1 - \gamma \left( {1 - \exp \left( { - \tau ({{r}_{b}}} \right)} \right)} \right].$

Моделирование угла рассеяния для неизотропных индикатрис производится согласно [13]: $\gamma $ определяет долю от полной площади под графиком индикатрисы, по которой затем находится соответствующий угол рассеяния. Рассеяние описывается в локальной системе координат, связанной с частицей, в которой угол рассеяния $\theta $ отсчитывается от направления движения частицы ${\mathbf{\omega }}$. Перевод направляющего вектора частицы в глобальную систему координат осуществлялся по матрице перехода от локальных координат к глобальным.

АРХИТЕКТУРА ВЫЧИСЛЕНИЙ

Расчеты производились на видеокарте Nvidia GeForce GTX 1660 Super на архитектуре CUDA. Введем управляющие величины, описывающие размерность задачи.

1. Количество ячеек внутри расчетной области по оси X – ${{N}_{x}}$ и по оси Y – ${{N}_{r}}$. Так как распределение термодинамических параметров в струе рассчитывалось в стороннем пакете, эти параметры фиксированы и равны 19. При регистрации сигнала интерес представляют ячейки на поверхности на полосе $\varphi = 0$ – количество таких ячеек ${{N}_{x}} + {\text{\;2\;}}{{N}_{r}}$.

2. Количество направлений, на которое разбивается телесный угол в π радиан над каждой внешней ячейкой на полосе $\varphi = 0$ – ${{N}_{{{\text{zen}}}}}$ и ${{N}_{{{\text{azim}}}}}$. При этом зенитный угол изменяется в диапазоне [0, π/2], а азимутальный – в диапазоне [0, π] (что обусловливается осевой симметрией задачи и позволяет произвести дополнительное осреднение по азимутальному углу в случае метода прямого счета). Параметры ${{N}_{{{\text{zen}}}}}$ и ${{N}_{{{\text{azim}}}}}$ равны 120.

3. Совокупный набор направлений регистрации излучения, построенных над всеми внешними ячейками ${{N}_{1}} = \left( {{{N}_{x}} + {\text{\;}}2{{N}_{r}}{\text{\;}}} \right){{N}_{{{\text{zen}}}}}{{N}_{{{\text{azim}}}}}$ (${{N}_{1}} \approx 820\,800$).

4. Количество фотонов ${{N}_{{{\text{min}}}}}$, разыгрываемых в ячейке с минимальной энергией, задает полное количество фотонов ${{N}_{2}}$ (${{N}_{2}} \approx 2\,{\kern 1pt} 007\,900$ при ${{N}_{{\min }}} = 50$, используемыx в схеме локальных оценок прямого уравнения, и ${{N}_{2}} \approx 80\,316\,000$ при ${{N}_{{\min }}} = 2000$ для схемы прямого расчета), разыгрываемых в расчетной области при применении методов прямого счета и локальных оценок прямого уравнения.

5. Количество траекторий в случае применения локальных оценок сопряженного уравнения ${{N}_{3}}$ ≈ 100.

В случае реализации метода прямого расчета необходимо организовать вычисление в параллельном режиме ${{N}_{2}}$ траекторий фотонов из внутренней области. В случае метода локальных оценок сопряженного уравнения необходимо рассчитывать ${{N}_{1}}{{N}_{3}}$ процессов, так как решения для каждого из ${{N}_{1}}$ направлений являются взаимонезависимыми. В случае метода локальных оценок прямого уравнения возможно несколько подходов к организации вычислений. Можно использовать схему как в случае метода прямого счета, моделируя ${{N}_{2}}$ процессов, и на каждом шаге рассеяния каждого фотона производить ${{N}_{1}}~$ вычислений локальных оценок для всех направлений регистрации. Однако в данном случае таких направлений слишком много, что не позволило организовать эффективный цикл на видеокарте, в том числе при использовании вложенных ядер. Поэтому использовалась схема расчета ${{N}_{1}}{{N}_{2}}$ процессов, при которой для каждого направления из ${{N}_{1}}$ запускается свое вычисление на ${{N}_{2}}$ траекториях из внутренней среды.

Пусть необходимо произвести параллельные вычисления $N$ процессов. Параметры вызова ядра CUDA следующие: используется одномерная сетка и одномерные блоки, количество нитей на блок – 128 (для используемой видеокарты такое количество нитей показало лучшую производительность, чем 256 и 512). Необходимое количество блоков вычисляется по соотношению $(N + 127){\text{/}}128$. Номер нити ${\text{thread}}$ на видеокарте вычисляется как ${\text{threadIdx}}.{\text{x\;}} + {\text{\;blockIdx}}.{\text{x}} \times {\text{blockDim}}.{\text{x}}$, вычисления осуществляются, если ${\text{thread}} < N$. При этом если полученное количество блоков превышает максимально допустимое, определяемое для каждой видеокарты параметром maxGridSize[0], то необходимо разбивать $N$ на части, организуя цикл на CPU для запуска ядра на видеокарте. Организация глубоких циклов (как в случае суммирования по $n$ в (9) и в (15)) может привести к необходимости уменьшения количества блоков. Математическое ожидание в (9) и (15) рассчитывалось с использованием атомарной операции atomiсAdd (в CUDA данная операция определена для переменных типа float), что не уступает по скорости записи данных в массив с последующим сложением на CPU. Генерация случайных чисел осуществлялась на GPU так, что для каждой из $N$ траекторий фотонов инициировался собственный генератор, порождавший уникальный ряд случайных чисел для реализации трассировки фотона согласно уравнениям (5) и (12).

РЕЗУЛЬТАТЫ

Прослежено формирование угловой зависимости силы излучения (1) при изменении положения наблюдателя ${{r}_{{{\text{obs}}}}}$ относительно струи независимо для трех образующих ее поверхностей. Поверхность $A$ соответствует плоскости $x = 0$, и сила излучения формируется на круге области начала струи. Поверхность $B$ соответствует боковой поверхности струи, сила излучения формируется на параболоиде вращения. Поверхность $C$ соответствует плоскости $x = L$ ($L$ – длина струи), и сила излучения формируется на круге конечного сечения струи. Дистанция до объекта (1 км) много больше размера струи, обход осуществляется по круговым траекториям и на всех графиках по оси абсцисс отложено угловое положение наблюдателя. Излучающими компонентами газовой фазы являются вода и углекислый газ. Коэффициенты поглощения получены из базы данных HITRAN [14], и их значения при различных температурах табулированы. От спектральных величин осуществлен переход к интегральным согласно групповой модели [13] в малом спектральном диапазоне. Средний по расчетному объему коэффициент поглощения газовой фазы составляет 0.057 м^–1. Выбраны четыре сечения рассеяния частиц (0.1, 10, 40.5, 80) × 10^–11 м², которые соответствуют среднеобъемным значениям альбедо $\varpi = \left( {0.021,~\,\,0.177,\,\,~0.465,\,\,~0.632} \right)$. На рис. 2 показаны угловые зависимости силы излучения нерассеянной компоненты при различных значениях альбедо.

Рис. 2.

Влияние альбедо среды на изменение угловой зависимости силы прямого (нерассеянного) излучения ${{I}_{0}}$ при использовании изотропной индикатрисы рассеяния для поверхностей A (а), B (б), C (в): 1 – $\varpi = 0.021$, 2 – 0.465, 3 – 0.632.

Наличие криволинейного паттерна в силе излучения с поверхности $B$ связано с количеством ячеек, на которые разбивается расчетная область. При использовании более мелкой сетки по оси $OX$ данный паттерн сглаживается, что требует, однако, большего времени расчета. Для получения общего решения в случае неаналоговых методов нерассеянная компонента складывается с рассеянным излучением. На рис. 3 приведено сравнение угловых зависимостей силы излучения рассеянных компонент, полученных с использованием локальных оценок в решении прямого и сопряженного уравнений при изотропной индикатрисе рассеяния $P({\mathbf{\omega }}{\kern 1pt} ' \to {\mathbf{\omega }}) = {1 \mathord{\left/ {\vphantom {1 {4\pi }}} \right. \kern-0em} {4\pi }}$. Значение ${{N}_{{\min }}}$, используемое в прямом методе, составляет 50.

Рис. 3.

Угловые зависимости силы излучения рассеянной компоненты по алгоритмам решения прямого (а)–(в) и сопряженного (г)–(е) уравнений с изотропной индикатрисой рассеяния для поверхностей A (а), (г), B (б), (д), C (в), (е): 1 – $\varpi = 0.021$, 2 – 0.465, 3 – 0.632.

Полученные по двум различным методам значения силы излучения количественно соответствуют друг другу. Из графиков можно наблюдать увеличение силы излучения рассеянной компоненты при увеличении альбедо. Результаты расчета угловой зависимости силы полного (рассеянного и нерассеянного) излучения по решению сопряженного уравнения в сопоставлении с результатами по прямому методу расчета представлены на рис. 4.

Рис. 4.

Угловые зависимости силы излучения по решению сопряженного уравнения (а)–(в) и прямому методу расчета (г)–(е) с изотропной индикатрисой рассеяния для: 1 – $\varpi = 0.021$, 2 – 0.465, 3 – 0.632.

Значение ${{N}_{{\min }}}$, используемое в методе прямого расчета, составляет 2000. При изотропной индикатрисе рассеяния данный метод демонстрирует несколько большие (примерно на 5%) значения силы излучения с поверхности, но в целом результаты хорошо согласуются друг с другом. Рассмотрим теперь угловые зависимости при использовании неизотропной индикатрисы Хеньи–Гринштейна, которая часто используется для описания индикатрис аэрозолей ${{P}_{{{\text{HG}}}}}\left( \theta \right) = \frac{{\left( {1 - {{g}^{2}}} \right){{{\left( {1 + {{g}^{2}} - 2g{\kern 1pt} \cos (\theta )} \right)}}^{{ - 3/2}}}}}{{4\pi }},$ где $g$ – параметр, определяющий вытянутость индикатрисы (фактор асимметрии). При увеличении рассеяния вперед профиль угловой зависимости силы излучения (рассеянной компоненты) претерпевает существенные изменения. На рис. 5 продемонстрировано изменение профиля силы излучения рассеянной компоненты при расчете по сопряженному уравнению.

Рис. 5.

Угловые зависимости рассеянной компоненты силы излучения по решению сопряженного уравнения и $\varpi = 0.021$ и различных параметрах асимметрии g: 1 – изотропное рассеяние (g = 0), 2 – g = 0.5, 3 – 0.9; для поверхностей A (а), B (б), C (в).

Сравнение различных способов расчета рассеянного излучения при неизотропной индикатрисе показывает плохую адаптируемость метода прямого расчета, занижающего значения по сравнению с методами локальных оценок в 2.5–3 раза при g = 0.9, что продемонстрировано на рис. 6.

Рис. 6.

Сравнение угловых зависимостей силы рассеянной компоненты излучения для трех расчетных методов с индикатрисой Хеньи–Гринштейна при $\varpi = 0.021$ и g = 0.5 (а)–(в) и 0.9 (г)–(е) для поверхностей A (а), (г), B (б), (д), C (в), (е): 1 – метод прямого расчета; 2 – локальная оценка, прямое уравнение; 3 – сопряженное уравнение.

Рассеянная компонента в методе прямого расчета представлена здесь как разность полного излучения и нерассеянного. На графиках отсутствует зависимость рассеянного излучения с боковой поверхности $B$ по данному методу ввиду наличия криволинейного паттерна в расчете нерассеянного излучения (обусловленного степенью детализации расчетной сетки). Методы, основанные на локальных оценках, по-прежнему хорошо согласуются друг с другом качественно, однако стоит отметить возрастающую дисперсию результата для решения прямого уравнения при увеличении неизотропности индикатрисы. Сравнение времен расчета представлено в таблице.

Сравнение времен расчета (в с) при использовании индикатрисы Хеньи–Гринштейна с параметром $g = 0.5$

$\varpi $	Метод прямого расчета, ${{N}_{{{\text{min}}}}} = 2000$	Локальная оценка (прямое уравнение), ${{N}_{{{\text{min}}}}} = 50$	Локальная оценка (сопряженное уравнение)
$0.021$	126.8	125 636	2560.6
$0.465$	126.8	144 719	3186.7
$0.465$	126.9	158 114	3720.5

Вычислительно более “дорогим” для данной задачи является метод локальных оценок при решении прямого уравнения. Следует обратить внимание на вопрос сходимости при вычислении математического ожидания в (9) и в (15). Выбрана площадка в середине поверхности $B$, и проанализирована сходимость значений рассеянного излучения при изменении управляющих параметров (${{N}_{{{\text{min}}}}}$ для локальных оценок прямого уравнения и ${{N}_{3}}$ для локальных оценок сопряженного уравнения). При этом на сетке регистрируемых направлений зенитный угол отсчитывается от нормали к поверхности, таким образом, угол в 0° соответствует направлению внутрь расчетного объема по нормали, а угол в 90° – направлению по касательной с быстрым выходом из объема.

Несмотря на то, что в обоих методах используются одни и те же функции расчета длины свободного пробега и оптической толщины, сходимость метода локальных оценок сопряженного уравнения не зависит от направления регистрации (одинаково низкая дисперсия на графиках 1B–3B, рис. 7), в то время как при решении прямого уравнения сходимость зависит от направления (возрастающая дисперсия на графиках 1D–3D при стремлении зенитного угла к 0°). Плохой сходимостью решения прямого уравнения объясняются пики на представленных ранее зависимостях рассеянного излучения при обходе поверхности. Такое поведение сходимости может обусловливаться и недостаточным качеством примененных алгоритмов трассировки луча.

Рис. 7.

Сходимость метода локальных оценок сопряженного уравнения при зенитном угле 0.375°, азимутальном – 0.75° (а); (б) сравнение зависимостей по зенитному углу методов локальных оценок прямого уравнения (D, ${{N}_{{{\text{min}}}}} = 10{\kern 1pt} 000$) и сопряженного уравнения (B, ${{N}_{3}} = 2000$): 1 – $\varpi = 0.021$, 2 – 0.465, 3 – 0.632.

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

ЗАКЛЮЧЕНИЕ

Таким образом, можно выполнить качественное сравнение рассмотренных методов расчета силы излучения с осесимметричного объема при наличии рассеяния. Методу решения прямого уравнения с использованием локальной оценки и методу прямого расчета сопутствует общее ограничение: выражение (10) для мощности, испускаемой элементарным объемом среды, соответствует предположению объемного высвечивания, т.е. малой оптической толщины такого объема, что может налагать довольно жесткие условия на дискретизацию расчетной области. Метод прямого расчета не подходит при работе с неизотропными индикатрисами, по крайней мере для сеток детектируемых направлений высокой дискретизации. Метод локальных оценок прямого уравнения способен учитывать неизотропность индикатрис рассеяния в среде, однако для малой дисперсии результата может понадобиться достаточно большое число ${{N}_{{\min }}}$. Общим свойством методов локальных оценок является принципиальная независимость от наличия симметрии в среде. Наиболее универсальным оказался метод решения сопряженного уравнения переноса с использованием локальной оценки. Во-первых, энергия, излучаемая объемом среды в данном методе, рассчитывается напрямую, без деления на отдельные элементарные излучающие объемы. Во-вторых, явный вид оценки (14) в данном методе интуитивно понятным образом обеспечивает вклад от всех источников излучения. С ее использованием наиболее просто учитывать вклад внешних источников излучения, видимых из текущей точки под определенным телесным углом.

Продемонстрировано значительное влияние величины альбедо однократного рассеяния на характеристики рассеянного излучения: при увеличении значения альбедо в 30 раз сила излучения от рассеянной компоненты в случае использования однородной индикатрисы возрастает в ~3 раза с поверхности B и в ~1.5 раза с поверхностей A и C. С увеличением доли рассеянного вперед из-лучения профиль силы излучения от поверхности B деформируется, что является важным фактором при оценке максимальной силы излучения рассеивающих струй, так как, несмотря на максимальную площадь проекции струи при положении наблюдателя в 90°, максимум силы излучения достигается при смещенном угле наблюдения (~15° и ~160°). Однако при расчете реальных струй двигателей данные угловые зависимости дополнительно деформируются излучением металлических стенок сопла – расчет подобной задачи также полностью описывается представленными уравнениями.

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

Марчук Г.И., Михайлов Г.А., Назарлиев М.А. и др. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука, 1976.
Yetzer K., Varnai T., Huang S. Intercomparison of 3D Radiation Codes, 2015. http://i3rc.gsfc.nasa.gov/
Михайлов Г.А., Лотова Г.З. Численно-статистическая оценка потока частиц с конечной дисперсией // До-кл. АН. 2012. Т. 447. № 1. С. 18.
Каблукова Е.Г., Каргин Б.А. Эффективные дискретно-стохастические модификации локальных оценок метода Монте-Карло для задач лазерного зондирования рассеивающих сред // Вычислительные технологии. 2012. Т. 17. № 3. С. 70.
Ухинов С.А. Методы Монте-Карло для решения задач теории переноса поляризованного излучения. Дис. … докт. физ.-мат. наук. Новосибирск: Институт вычислительной математики и математической геофизики Сиб. отд. РАН, 2010.
Wang Z., Cui S., Zhang Z. et al. Theoretical Extension of Universal Forward and Backward Monte Carlo Radiative Transfer Modeling for Passive and Active Polarization Observation Simulations // J. Quant. Spectrosc. Radiat. Transf. 2019. V. 235. P. 81.
Wang Z., Cui S., Yang J. et al. A Novel Hybrid Scattering Order-dependent Variance Reduction Method for Monte Carlo Simulations of Radiative Transfer in Cloudy Atmosphere // J. Quant. Spectrosc. Radiat. Transf. 2017. V. 189. P. 283.
Deutschmann T., Beirle S., Frieß U. et al. The Monte Carlo Atmospheric Radiative Transfer Model McArtim: Introduction and Validation of Jacobians and 3D Features // J. Quant. Spectrosc. Radiat. Transf. 2011. V. 112. № 6. P. 1119.
Modest M.F. Backward Monte Carlo Simulations in Radiative Heat Transfer // J. Heat Transfer. 2003. V. 125. № 1. P. 57.
Котов Д.В., Суржиков С.Т. Локальная оценка направленной излучательной способности светорассеивающих объемов методом Монте-Карло // ТВТ. 2007. Т. 45. № 6. С. 885.
Григорьев И.С. Моделирование влияния процессов рассеяния на силу излучения высокотемпературной струи по методу Монте-Карло // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2020. Т. 5. № 92. С. 28.
Соболь И.М. Численные методы Монте-Карло. М.: Наука, 1973. 311 с.
Суржиков С.Т. Тепловое излучение газов и плазмы. М.: МГТУ им. Н.Э. Баумана, 2004. 544 с.
Kochanov R.V., Gordon I.E., Rothman L.S. et al. HITRAN Application Programming Interface (HAPI): A Comprehensive Approach to Working with Spectroscopic Data // J. Quant. Spectrosc. Radiat. Transf. 2016. V. 177. P. 15.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Теплофизика высоких температур

Архивы выпусков Информация о журнале Отправить рукопись в журнал