Журнал вычислительной математики и математической физики, 2022, T. 62, № 5, стр. 861-871

Аналитические решения модельных кинетических уравнений переноса излучения и уравнения энергии

Н. Я. Моисеев^1, *, В. М. Шмаков^1, **

¹ ФГУП “РФЯЦ-ВНИИТФ им. акад. Е.И. Забабахина”
456770 Челябинская обл., Снежинск, а/я 245, ул. Васильева, 13, Россия

^* E-mail: nik.moiseev.43@mail.ru
^** E-mail: v.m.shmakov@vniitf.ru

Поступила в редакцию 11.03.2021
После доработки 12.10.2021
Принята к публикации 14.01.2022

EDN: DFDCGS
DOI: 10.31857/S004446692205009X

Аннотация

Рассмотрена модельная система нестационарных кинетических уравнений переноса теплового излучения и уравнения энергии в многогрупповом изотропном приближении в средах с постоянными коэффициентами поглощения и кусочно-постоянной функцией Планка в группах. Для модельной системы уравнений получены аналитические решения для плоской геометрии и для шара. Аналитические решения получены путем перехода от решения сложной системы уравнений к решению более простой, которая имеет известное аналитическое решение. Обратный переход дает решение сложных уравнений. Приведены решения тестовых задач для плоского слоя и для шара. Библ. 6. Фиг. 10. Табл. 2.

Ключевые слова: кинетическое уравнение переноса теплового излучения и уравнения энергии, аналитические решения.

1. ВВЕДЕНИЕ

Кинетическое нестационарное уравнение переноса теплового излучения является интегро-дифференциальным уравнением. В работах [1]–[3] рассматриваются подходы к решению этого уравнения в различных предположениях относительно коэффициентов, что позволяет получать упрощенные уравнения переноса. Так, если коэффициенты в уравнении переноса положить равными нулю, то уравнение записывается в простейшей форме и моделирует перенос излучения в вакууме. Для уравнения переноса излучения в вакууме А.В. Вронским были получены аналитические решения, которые используются для отладки численных методов. Однако, если решаются уравнения переноса теплового излучения и уравнения энергии в средах с поглощением и переизлучением, то этих решений для отладки численных методов недостаточно.

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

2. ТОЧНЫЕ РЕШЕНИЯ МОДЕЛЬНЫХ КИНЕТИЧЕСКИХ НЕСТАЦИОНАРНЫХ УРАВНЕНИЙ ПЕРЕНОСА ИЗЛУЧЕНИЯ И УРАВНЕНИЯ ЭНЕРГИИ

2.1. Постановка задачи

Предположим, что коэффициенты рассеяния и функция источника равны нулю, плотность вещества $\rho = 1$. Систему кинетических нестационарных уравнений переноса теплового излучения и уравнения энергии в многогрупповом изотропном приближении запишем в виде [2]

$\frac{1}{c}\frac{{\partial {{I}_{g}}}}{{\partial t}} + {\mathbf{\Omega }}\nabla {{I}_{g}} + {{\alpha }_{g}}{{I}_{g}} = 0.5{{\alpha }_{g}}{{B}_{g}},$

(2.1)

$\begin{gathered} \frac{{dE}}{{dt}} = \mathop \sum \limits_{g = 1}^G {{\alpha }_{g}}\left( {{{U}_{g}} - {{B}_{g}}} \right), \\ {{U}_{g}} = \mathop \smallint \limits_{ - 1}^1 {{I}_{g}}d\mu , \\ \end{gathered} $

${\mathbf{\Omega }}\nabla = \left\{ {\begin{array}{*{20}{l}} {\mu \frac{\partial }{{\partial r}} + \frac{{1 - {{\mu }^{2}}}}{r}\frac{\partial }{{\partial \mu }}--{\text{для\;сферически - симметричной\;геометрии}},} \\ {\mu \frac{\partial }{{\partial x}} - {\text{для\;плоской\;геометрии}}.~~~~~~~~~~~~~~~~~~~~} \end{array}} \right.$

Здесь $x$, $r$ и $t$ – независимые переменные по пространству и времени соответственно, $\mu $ – косинус угла между направлением движения частиц и осью $x$ или $r$, $ - 1 \leqslant \mu \leqslant 1$, $g$ – индекс энергетической группы, ${{I}_{g}}$, $E$ – неизвестные функции, $E\left( {t,x} \right)$ – удельная внутренняя энергия вещества, $T\left( {t,x} \right)$ – температура вещества, $с$ – скорость света, ${{I}_{g}}\left( {t,x,\mu } \right)$ – спектральная интенсивность энергии излучения в группе $g$, ${{U}_{g}}\left( {t,x} \right)~$ – спектральная плотность энергии излучения, умноженная на скорость света, ${{\alpha }_{g}}\left( {x,T} \right)$ – коэффициент поглощения,

${{B}_{g}}\left( T \right) = \frac{{8\pi }}{{{{c}^{2}}{{{\bar {h}}}^{3}}}}\mathop \smallint \limits_{{{\varepsilon }_{g}}}^{{{\varepsilon }_{{g + 1}}}} \frac{{{{\varepsilon }^{3}}}}{{{\text{exp}}\left( {\varepsilon {\text{/}}T} \right) - 1}}d\varepsilon $

есть интенсивность равновесного излучения (функция Планка), умноженная на скорость света, $\bar {h}$ – постоянная Планка. Разностная сетка по энергии фотонов включает $G$ независимых групп с энергиями ${{\varepsilon }_{1}}, \ldots ,{{\varepsilon }_{g}}, \ldots ,{{\varepsilon }_{G}}$. Разностная сетка по переменной $\mu $ с центрами ${{\mu }_{m}}$ = $0.5\left( {{{\mu }_{{m + 1/2}}} + {{\mu }_{{m - 1/2}}}} \right)$ и шагами $\Delta \mu = {{\mu }_{{m + 1/2}}}$ – ${{\mu }_{{m - 1/2}}}$, $m = 1, \ldots ,M$, включает $M$ направлений движения частиц (фотонов).

Система уравнений (2.1) замыкается уравнением состояния вещества в форме $E = E\left( T \right)$. Начальные условия: $T\left( {0,x} \right) = {{T}_{0}}\left( x \right)$, ${{I}_{g}}\left( {0,x,\mu } \right)$ = ${{I}_{{g,0}}}\left( {x,\mu } \right)$. Требуется найти решение задачи Коши для $t > 0$ в областях $D = \{ {{x}_{L}} \leqslant x \leqslant {{x}_{R}}$, $ - 1 \leqslant \mu \leqslant 1\} $, $D = \{ 0 \leqslant r \leqslant {{r}_{0}}$, $ - 1 \leqslant \mu \leqslant 1\} $ в случае плоской или сферически-симметричной геометрии соответственно.

2.2. Точные решения модельного уравнения переноса частиц в средах с поглощением

Для упрощения исследования системы уравнений (2.1) коэффициенты поглощения положим равными некоторым постоянным ${{\alpha }_{g}}\left( {x,T} \right) = {\text{const}}$. Функцию Планка аппроксимируем кусочно-постоянной функцией ${{B}_{g}}\left( T \right) = {\text{const}}$ в каждой группе. Уравнение состояния вещества возьмем в форме $E = 0.81T$. Систему кинетических уравнений переноса излучения и уравнения энергии с постоянными коэффициентами поглощения и кусочно-постоянной функцией Планка будем называть модельной системой уравнений. Для модельной системы уравнений рассмотрим решение задач Коши в плоской и сферически-симметричной геометриях. Для плоской геометрии начальные данные задаются в слое $\left| x \right| \leqslant {{x}_{0}}$, для сферической – в шаре с радиусом ${{r}_{0}}$. В первом случае задачу будем называть задачей о плоском слое, во втором – задачей о шаре. Пусть в начальный момент времени на всей оси x задана температура вещества $T\left( {0,x} \right)$ = ${{T}_{0}}\left( x \right)$ и задано однородно распределение частиц (фотонов)

${{I}_{g}}\left( {0,x,\mu } \right) = \left\{ {\begin{array}{*{20}{l}} {{{I}_{{g,0}}}\left( {x,\mu } \right),\quad \left| x \right| \leqslant {{x}_{0}},\quad {{I}_{{g,0}}} \geqslant 0.5{{B}_{g}}\left( {{{T}_{0}}} \right),} \\ {0.5{{B}_{g}}\left( {{{T}_{0}}} \right),\quad \left| x \right| > {{x}_{0}}.} \end{array}} \right.$

Требуется определить интенсивность излучения, плотность энергии излучения и температуру вещества для $t > 0$.

Аналитические решения модельной системы уравнений (2.1) для решения поставленной задачи получим следующим образом. Выполнив в системе уравнений (2.1) преобразование

(2.2)

$\begin{gathered} {{I}_{g}} = {{J}_{g}}\exp \left( { - ct{{\alpha }_{g}}} \right) + 0.5{{B}_{g}}, \\ \widehat {{{U}_{g}}} = \mathop \smallint \limits_{ - 1}^1 {{J}_{g}}d\mu \\ \end{gathered} $

получим эквивалентную систему уравнений

(2.3)

$\begin{gathered} \frac{1}{c}\frac{{\partial {{J}_{g}}}}{{\partial t}} + {\mathbf{\Omega }}\nabla {{J}_{g}} = 0, \\ \frac{{dE}}{{dt}} = \mathop \sum \limits_{g = 1}^G {{\alpha }_{g}}\widehat {{{U}_{g}}}\exp \left( { - ct{{\alpha }_{g}}} \right) \\ \end{gathered} $

с начальными данными для уравнения переноса в (2.3)

${{J}_{g}}\left( {0,x,\mu } \right) = \left\{ {\begin{array}{*{20}{l}} {{{I}_{{g,0}}} - 0.5{{B}_{g}},\quad \left| x \right| \leqslant {{x}_{0}},} \\ {~0,\quad \left| x \right| > {{x}_{0}}.} \end{array}} \right.$

Уравнение переноса в (2.3) можно интерпретировать как уравнение переноса излучения в вакууме, точные решения которого для интенсивности и плотности энергии излучения записываются в виде

(2.4)

${{J}_{g}}\left( {t,x,\mu } \right) = {{J}_{{g,0}}}\xi \left( {t,x} \right),\quad \widehat {{{U}_{g}}}\left( {t,x} \right) = {{U}_{{g,0}}}\xi \left( {t,x} \right)$

соответственно. Подставив (2.4) в выражения (2.2), получим для интенсивности и плотности энергии излучения решения модельного уравнения переноса в (2.1), которые запишем в виде

${{I}_{g}}\left( {t,x,\mu } \right) = \gamma \xi \left( {t,x} \right){{I}_{{g,0}}} + 0.5\left( {1 - \gamma \xi \left( {t,x} \right)} \right){{B}_{g}},$

(2.5)

${{U}_{g}}\left( {t,x} \right) = \gamma \xi \left( {t,x} \right){{U}_{{g,0}}} + \left( {1 - \gamma \xi \left( {t,x} \right)} \right){{B}_{g}},$

$\gamma = {\text{exp}}\left( { - ct{{\alpha }_{g}}} \right)$

соответственно. Отметим, что решения уравнения переноса излучения в вакууме находятся в плоскости $\left( {t,x} \right)$ в областях, которые отделены друг от друга характеристиками

$dx{\text{/}}dt = \pm c,$

выходящими из точек ${{x}_{0}}~$ и $ - {{x}_{0}}$. Функция $\xi \left( {t,x} \right)$ вычисляется в задаче о плоском слое в шести областях следующим образом. Если $ct \leqslant {{x}_{0}}$, то

$\xi \left( {t,x} \right) = \left\{ \begin{gathered} 0,\quad - {\kern 1pt} \infty < x < - {{x}_{0}} - ct,\quad {{x}_{0}} + ct < x < \infty ,\quad {\text{области}}\;I\;{\text{и}}\;II, \hfill \\ \frac{{{{x}_{0}} + ct + x}}{{2ct}},\quad - {\kern 1pt} {{x}_{0}} - ct \leqslant x < - {{x}_{0}} + ct,\quad {\text{область}}\;IV, \hfill \\ 1,\quad - {\kern 1pt} {{x}_{0}} + ct \leqslant x \leqslant ~~{{x}_{0}} - ct,\quad {\text{область}}\;III,~ \hfill \\ \frac{{{{x}_{0}} + ct - x}}{{2ct}},\quad {{x}_{0}} - ct < x \leqslant {{x}_{0}} + ct,\quad {\text{область}}\;V. \hfill \\ \end{gathered} \right.$

Если $ct > {{x}_{0}}$, то

$\xi \left( {t,x} \right) = \left\{ \begin{gathered} 0,\quad - {\kern 1pt} \infty < x < - {{x}_{0}} - ct,\quad {{x}_{0}} + ct < x < \infty ,\quad {\text{области}}\;I\;{\text{и}}\;II, \hfill \\ \frac{{{{x}_{0}} + ct + x}}{{2ct}},\quad - {\kern 1pt} {{x}_{0}} - ct \leqslant x < ~~{{x}_{0}} - ct,\quad {\text{область}}\;IV, \hfill \\ \frac{{{{x}_{0}}}}{{ct}},\quad {{x}_{0}} - ct \leqslant x \leqslant ~~ - {{x}_{0}} + ct,\quad {\text{область}}\;VI, \hfill \\ \frac{{{{x}_{0}} + ct - x}}{{2ct}},\quad - {\kern 1pt} {{x}_{0}} + ct < x \leqslant {{x}_{0}} + ct,\quad {\text{область}}\;V. \hfill \\ \end{gathered} \right.$

Функция $\xi \left( {t,r} \right)$ вычисляется в задаче о шаре радиуса ${{r}_{0}}$ в четырех областях следующим образом. Если $ct \leqslant {{r}_{0}}$, то

$\xi \left( {t,x} \right) = \left\{ {\begin{array}{*{20}{l}} {~1,\quad 0 \leqslant r < {{r}_{0}} - ct,\quad {\text{область}}\;I,~} \\ {~\frac{{r_{0}^{2} - {{{\left( {r - ct} \right)}}^{2}}}}{{4rct}},\quad {{r}_{0}} - ct \leqslant r < {{r}_{0}} + ct,\quad {\text{область}}\;III,} \\ {~0,\quad {{r}_{0}} + ct \leqslant r < \infty ,\quad {\text{область}}\;II.} \end{array}} \right.$

Если $ct > {{r}_{0}}$, то

$\xi \left( {t,x} \right) = \left\{ {\begin{array}{*{20}{l}} {~0,\quad 0 \leqslant r < - {{r}_{0}} + ct,\quad {\text{область}}\;IV,} \\ {~\frac{{r_{0}^{2} - {{{\left( {r - ct} \right)}}^{2}}}}{{4rct}},~\quad - {\kern 1pt} {{r}_{0}} + ct \leqslant r < {{r}_{0}} + ct,\quad {\text{область}}\;III,} \\ {~0,\quad {{r}_{0}} + ct \leqslant r < \infty ,\quad {\text{область}}\;II.} \end{array}} \right.$

2.3. Аналитические решения модельного уравнения энергии в средах с поглощением

Аналитические решения модельного уравнения энергии получим следующим образом. Подставив выражение (2.5) для плотности энергии излучения в уравнение энергии в (2.1), получим обыкновенное дифференциальное уравнение в виде

(2.6)

$\frac{{dE}}{{dt}} = \mathop \sum \limits_{g = 1}^G {{\alpha }_{g}}\left( {{{U}_{{g,0}}} - {{B}_{g}}} \right)\xi \left( {t,x} \right){\text{exp}}\left( { - ct{{\alpha }_{g}}} \right).$

Проинтегрировав уравнение (2.6) на интервале $\left[ {0,t} \right]$, получим для вычисления удельной внутренней энергии вещества выражение

(2.7)

$E\left( {t,x} \right) = E\left( {0,x} \right) + \mathop \sum \limits_{g = 1}^G \,\mathop \smallint \limits_0^t {{\alpha }_{g}}\left( {{{U}_{{g,0}}} - {{B}_{g}}} \right)\xi \left( {\tau ,x} \right)\exp \left( { - c\tau {{\alpha }_{g}}} \right)d\tau .$

Если $\xi \left( {t,x} \right) = 0$, то удельная внутренняя энергия вещества остается постоянной $E\left( {t,x} \right)$ = $E\left( {0,x} \right)$ в задаче о плоском слое в областях I и II, в задаче о шаре – в областях II и IV. Если $\xi \left( {t,x} \right) = 1$, то интеграл в (2.7) берется в квадратурах точно. Взяв интеграл, получим для вычисления удельной внутренней энергии вещества в задачах о плоском слое в области III и о шаре в области I выражение

(2.8)

$E\left( {t,x} \right) = E\left( {0,x} \right) + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \left( {1 - \gamma } \right)\left( {{{U}_{{g,0}}} - {{B}_{g}}} \right),\quad \gamma = \exp \left( { - ct{{\alpha }_{g}}} \right).$

Если $\xi \ne 0$, то удельная внутренняя энергия вещества вычисляется в задаче о плоском слое в областях $I$, $V$ и VI из выражений

(2.9)

$E\left( {t,x} \right) = E\left( {0,x} \right) + 0.5\,\mathop \sum \limits_{g = 1}^G {{\alpha }_{g}}\left( {{{U}_{{g,0}}} - {{B}_{g}}} \right)\mathop \smallint \limits_0^t \left( {{{x}_{0}} + c\tau \pm x} \right)\frac{{\exp \left( { - c\tau {{\alpha }_{g}}} \right)}}{{c\tau }}d\tau ,~$

(2.10)

$E\left( {t,x} \right) = E\left( {0,x} \right) + \mathop \sum \limits_{g = 1}^G \,{{\alpha }_{g}}\left( {{{U}_{{g,0}}} - {{B}_{g}}} \right)\mathop \smallint \limits_0^t {{x}_{0}}\frac{{\exp \left( { - c\tau {{\alpha }_{g}}} \right)}}{{c\tau }}d\tau ~$

соответственно. Здесь плюс берется в области IV, минус – в области V. Удельная внутренняя энергия вещества вычисляется в задаче о шаре в области $III$ из выражения

(2.11)

$E\left( {t,x} \right) = E\left( {0,x} \right) + \mathop \sum \limits_{g = 1}^G \,{{\alpha }_{g}}\left( {{{U}_{{g,0}}} - {{B}_{g}}} \right)\mathop \smallint \limits_0^t \frac{{r_{0}^{2} - {{{\left( {r - c\tau } \right)}}^{2}}}}{{4r}}\frac{{\exp \left( { - c\tau {{\alpha }_{g}}} \right)}}{{c\tau }}d\tau .~$

Интегралы в выражениях (2.9)–(2.11) являются интегралами экспоненциального типа, не интегрируются в элементарных функциях, и имеют устранимую особенность в нуле. Пусть удельная внутренняя энергия вещества в задаче о плоском слое вычисляется в точке $A\left( {t,x} \right)$, которая находится в области $VI$, где $t > {{t}_{0}}$. Если $\left| x \right| \leqslant {{x}_{0}}$ или $\left| x \right| > {{x}_{0}}$, то интеграл в уравнении (2.9) можно представить в виде суммы трех интегралов в областях III, V, VI или II, V, VI на интервалах по времени $\left[ {0,{{t}_{1}}} \right]$, $\left[ {{{t}_{1}},{{t}_{2}}} \right]$, $\left[ {{{t}_{2}},t} \right]$ соответственно. Здесь ${{t}_{1}} = \left( {{{x}_{0}} - x} \right){\text{/}}c$, ${{t}_{2}} = \left( {{{x}_{0}} + x} \right){\text{/}}c$ – это времена прихода возмущений по характеристикам из точек ${{x}_{0}}$ и $ - {{x}_{0}}$ в точку $x$ соответственно. При таком подходе к вычислению интегралов особенность в нуле устраняется за счет того, что функция $\xi \left( {t,x} \right) = 1$ в области $III$ и $\xi \left( {t,x} \right) = 0$ в областях $I,II$. Особенности в нуле при вычислении интегралов в уравнениях (2.10) и (2.11) устраняются аналогичным образом. Экспоненциальные интегралы вычисляются по какой-либо квадратурной формуле с заданной точностью. Если интегралы вычислять по квадратурным формулам правых прямоугольников, которые обозначим как (С1), методом средних (С2) или по двухточечной формуле Гаусса (С3), то погрешности вычисления будут порядка $O\left( \tau \right)$, $O({{\tau }^{2}})$, $O({{\tau }^{3}})$ соответственно.

3. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ РЕШЕНИЙ МОДЕЛЬНЫХ ЗАДАЧ

Аналитические решения модельных уравнений (2.1) проверялись путем решения простейших задач по переносу излучения в задачах о плоском слое и шаре. По энергии фотонов используется 15-групповое приближение из работы [6] с умноженными на 10 границами интервалов: [0, 3, 6, 8, 12, 15, 18, 24, 27, 30, 40, 50, 70, 90, 110, 150]. Средние значения энергии в группах относятся к центрам интервалов и равны полусумме значений энергий на границах интервалов.

В расчетах использовалась модельная функция Планка: [0.029, 0.202, 0.391, 0.609, 0.813, 0.927, 1.0, 0.977, 0.926, 0.762, 0.489, 0.2057, 0.051, 0.010, 0.001]. Модельная функция Планка построена путем усреднения обычной функции Планка в диапазоне температур $\left[ {0,{{T}_{{{\text{max}}}}}} \right]$ для ${{T}_{{{\text{max}}}}} = 10$ кэВ и нормированием по максимальному значению. Для счета задач нормированная функция Планка умножается на положительный коэффициент. В дальнейшем такие функции будем различать по этим коэффициентам. Например, запись ${{B}_{g}} = 10$ определяет функцию Планка, которая получена из нормированной функции, умноженной на 10.

Решения находились на равномерной разностной сетке с шагом интегрирования по пространству равным $h = 0.002$ в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2. Удельная внутренняя энергия вещества вычислялась с контролем точности по квадратурным формулам С1, С2, С3. Контроль точности осуществлялся по результатам двух расчетов с шагами интегрирования по времени $\tau $ и $0.5\tau $ из анализа условия

$\left| {I\left( \tau \right) - I\left( {0.5\tau } \right)} \right| < \varepsilon .$

Здесь $\varepsilon $ – константа сходимости. Если условие выполняется, то интеграл полагается равным $I\left( {0.5\tau } \right)$. Если условие не выполняется, то шаг интегрирования по времени уменьшается в два раза и вычисление интеграла повторяется [4]. Начальный шаг интегрирования ${{t}_{0}} = {{x}_{0}}{\text{/}}c$ = = $6.6\left( 6 \right){\text{е}} - 5$.

3.1. Перенос излучения в оптически прозрачных средах с поглощением. Задача о плоском слое

Задача 1. Коэффициенты поглощения равны ${{\alpha }_{g}} = 0$ и ${{\alpha }_{g}} = 10$ при переносе излучения в вакууме и в среде с поглощением соответственно. Модельные функции Планка в группах равны ${{B}_{g}} = 0$. Начальные данные: если $\left| x \right| \leqslant {{x}_{0}} = 0.2$ см, то ${{I}_{0}} = 500$ и ${{U}_{0}} = 1000$, иначе ${{I}_{0}} = 0$, ${{U}_{0}} = 0$, температура вещества ${{T}_{0}} = 0.001$ кэВ. Найти распределения температуры вещества и плотности энергии излучения на интервале [–0.8, 0.8] в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2.

На фиг. 1–10 введены обозначения: V – обозначение вакуума, t в легендах на фигурах – отношение $t{\text{/}}{{t}_{0}}$.

Фиг. 1.

Зависимости от $x{\text{/}}{{x}_{0}}$ температуры излучения: (а) – в среде с поглощением, (б) – в вакууме.

Фиг. 2.

Зависимости от $x{\text{/}}{{x}_{0}}$ плотности энергии излучения: (а) – в среде с поглощением, (б) – в вакууме.

Фиг. 3.

Зависимости от $x{\text{/}}{{x}_{0}}$ температуры вещества в среде с коэффициентом поглощением ${{\alpha }_{g}} = 10$: (а) – ${{B}_{g}} = 0$, (б) – ${{B}_{g}} = 100$.

Фиг. 4.

Зависимости от $\tau {\text{/}}{{\tau }_{0}}$ удельной внутренней энергии вещества: кружки – С1, квадраты – С2, треугольники – С3.

Фиг. 5.

Зависимости от $x{\text{/}}{{x}_{0}}$: (а) – удельной внутренней энергии вещества, (б) – температуры вещества.

Фиг. 6.

Зависимости от $\tau {\text{/}}{{\tau }_{0}}$ удельной внутренней энергии вещества.

Фиг. 7.

Профили температур излучения: (а) – в среде с поглощением, (б) – в вакууме.

Фиг. 8.

Профили плотности энергии излучения: (а) – в среде с поглощения, (б) – в вакууме.

Фиг. 9.

Профили температур вещества в среде с коэффициентом поглощения ${{\alpha }_{g}} = 10$ и функцией Планка ${{B}_{g}} = 0$.

Фиг. 10.

Зависимости от $x{\text{/}}{{x}_{0}}$ температуры излучения: 1–3 – в вакууме, 4–6 – в среде с поглощением. Сплошные линии – точные решения, штриховые – численные решения.

На фиг. 1 представлены зависимости от $x{\text{/}}{{x}_{0}}~$ температуры излучения в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2 в среде с поглощением и в вакууме.

На фиг. 2 представлены зависимости от $x{\text{/}}{{x}_{0}}$ плотности энергии излучения в среде с поглощением и в вакууме в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2.

Из анализа графиков на фиг. 1 и 2 следует, что положения фронтов тепловых волн излучения в вакууме и в среде с поглощением совпадают между собой. Температуры и плотности энергии излучения существенно различаются.

На фиг. 3 представлены зависимости от $x{\text{/}}{{x}_{0}}$ температуры вещества в среде с коэффициентом поглощения в группах ${{\alpha }_{g}} = 10$ и функциями Планка ${{B}_{g}} = 0$, ${{B}_{g}} = 100$ в моменты времени $t{\text{/}}{{t}_{0}}$ = = 0.5, 1, 1.5, 2.

Точные максимальные значения температуры, рассчитанные в задачах с функциями Планка ${{B}_{g}} = 0$, ${{B}_{g}} = 100$ в момент времени $t{\text{/}}{{t}_{0}} = 0.5$, равны $T\left( {{{B}_{g}} = 0} \right)$ = 3.90297875820097 и $T\left( {{{B}_{g}} = 100} \right)$ = = 2.89551187902263 соответственно.

На фиг. 4 представлены зависимости от $\tau {\text{/}}{{\tau }_{0}}$ удельной внутренней энергии вещества в расчетах на сходимость в момент времени $t{\text{/}}{{t}_{0}} = 0.5$ при вычислении интегралов по квадратурным формулам С1, С2, С3.

В табл. 1 приведены удельные внутренние энергии вещества, шаги интегрирования и число шагов по времени при вычислении интегралов в точке $x = 0.011$ с константами сходимости $\varepsilon = 0.1$, $\varepsilon = 0.0001$ и функцией Планка ${{B}_{g}} = 0$ в момент времени $t{\text{/}}{{t}_{0}} = 0.5$.

Таблица 1

Метод	E $\varepsilon = 0.1$	${{\tau }_{1}}$	Число шагов	E $\varepsilon = 0.0001$	${{\tau }_{2}}$	Число шагов
С1	3.11147	1.4167 е-6	32	3.16050	2.034е-9	16 384
С2	3.12792	1.66(2.6) е-5	2	3.16057	5.208е-7	64
С3	3.15989	3.33(2.3) е-5	1	3.16059	8.333е-6	4

3.2. Перенос излучения в оптически плотных средах с поглощением. Задача о плоском слое

Задача 2. Коэффициенты поглощения в группах равны ${{\alpha }_{g}} = 10\,000$, функция Планка – ${{B}_{g}} = 0$, удельная внутренняя энергия вещества – ${{E}_{0}} = 0.00081$, температура – ${{T}_{0}} = 0.001$ кэВ. Плотность энергии излучения: если $\left| x \right| \leqslant {{x}_{0}} = 0.2$ см, то ${{I}_{0}} = 500$ и ${{U}_{0}} = 1000$, иначе ${{I}_{0}} = 0$, ${{U}_{0}} = 0$, удельная внутренняя энергия излучения – ${{E}_{{0,U}}} = 5$. Требуется определить температуру вещества и излучения в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2.

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

На фиг. 5 представлены зависимости от $x{\text{/}}{{x}_{0}}$ удельной внутренней энергии и температуры вещества в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2.

Результаты расчетов показали, что удельная внутренняя энергия излучения полностью перешла (преобразовалась) в удельную внутреннюю энергию вещества. Удельная внутренняя энергия вещества, рассчитанная по точным формулам (2.6) без учета ${{E}_{0}}$, равна ${{E}_{0}} = 5$, температура вещества – T = 6.17383950617284.

На фиг. 6 представлены зависимости от $\tau {\text{/}}{{\tau }_{0}}$ удельной внутренней энергии вещества в расчетах на сходимость в момент времени $t{\text{/}}{{t}_{0}} = 0.5$ при вычислении интегралов по квадратурным формулам С1, С2, С3.

В табл. 2 приведены удельные внутренние энергии вещества, шаги интегрирования и число шагов по времени при вычислении интегралов в точке $x = 0.011$ с константами сходимости $\varepsilon = 0.1$, $\varepsilon = 0.0001$ и функцией Планка ${{B}_{g}} = 0$ в момент времени $t{\text{/}}{{t}_{0}} = 0.5$.

Таблица 2

Метод	E $\varepsilon = 0.1$	${{\tau }_{1}}$	Число шагов	E $\varepsilon = 0.0001$	${{\tau }_{2}}$	Число шагов
С1	4.92409	1.0 е-9	32 768	4.92409	5.0е-10	65 536
С2	4.95067	1.63е-8	2048	4.99980	5.0е-10	32 768
С3	4.9850	6.5е-8	512	4.99999	8.0е-9	4096

Вычисление интегралов по квадратурным формулам С1, С2, С3 можно интерпретировать как численное решение уравнения энергии по разностным схемам Эйлера первого порядка, предиктор-корректор второго порядка и по схеме повышенного третьего порядка точности соответственно. Из анализа графиков на фиг. 4, 6 и из табл. 1, 2 следует, что для достижения заданной точности с константами сходимости $\varepsilon = 0.1$ и $\varepsilon = 0.0001$ по схеме С3 число шагов в 2, 16 и 4, 8 раз меньше, чем по схеме С2, и в 32, 4096 и в 64, 16 раз меньше, чем по схеме С1 в задачах 1 и 2 по переносу излучения в прозрачных и плотных средах соответственно. По схеме С2 число шагов в 16, 256 и в 16, 2 раза меньше, чем по схеме С1. Поэтому разностные схемы Эйлера первого порядка точности для решения уравнения энергии не эффективны и малопригодны для серийного счета задач. Схемы С2, С3 существенно эффективнее, чем схемы первого порядка точности: шаг интегрирования по времени на один, два порядка больше, чем в схеме первого порядка. Обе схемы позволяют получить результаты с заданной точностью за приемлемое время счета.

3.3. Перенос излучения в оптически прозрачных средах с поглощением. Задача о шаре

Задача 3. Начальные данные такие же, как в задаче 1 и заданы в шаре $r \leqslant {{r}_{0}} = 0.2$ см. Найти распределения температуры вещества и плотности энергии излучения в шаре радиуса $r = 0.6$ см в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2. Шаг интегрирования по пространству $h = 0.001$. Начальный шаг интегрирования ${{t}_{0}} = 6.6\left( 6 \right){\text{е}} - 5$.

На фиг. 7 представлены профили температур излучения в среде с коэффициентом поглощения ${{\alpha }_{g}} = 10$ и функцией Планка ${{B}_{g}} = 0$ и в вакууме в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2.

На фиг. 8 представлены профили плотности энергии излучения в среде с коэффициентом поглощения ${{\alpha }_{g}} = 10$ и функцией Планка ${{B}_{g}} = 0$ и в вакууме в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2.

На фиг. 9 представлены профили температур вещества при переносе излучения в среде с коэффициентом поглощения ${{\alpha }_{g}} = 10$ и функцией Планка ${{B}_{g}} = 0$ в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5, 2.

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ МОДЕЛЬНОЙ ЗАДАЧИ МОДИФИЦИРОВАННЫМ МЕТОДОМ РАСЩЕПЛЕНИЯ

В разделе для сравнения с точными решениями приведены численные решения модельной задачи 1, которые получены модифицированным методом расщепления, [5]. Численные решения уравнения переноса излучения получены на фиксированной разностной сетке с шагом по пространству $h = 0.002$ и с коэффициентом запаса устойчивости $k = 0.8$.

На фиг. 10 представлены профили температуры излучения в вакууме и в среде с поглощением в моменты времени $t{\text{/}}{{t}_{0}}$ = 0.5, 1, 1.5.

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

5. ЗАКЛЮЧЕНИЕ

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

Аналитические решения модельных задач подтверждены численными решениями, которые получены модифицированным методом расщепления.

Показано, что разностные схемы Эйлера первого порядка точности для решения уравнения энергии не эффективны и малопригодны.

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

Бай Ши-и. Динамика излучающего газа. М.: Мир, 1968.
Четверушкин Б.Н. Математическое моделирование задач динамики излучающего газа. М.: Наука, 1985.
Сушкевич Т.А. Математические модели переноса излучения. М.: Бином, 2006.
Каханер Д., Моулер К., Неш С. Численные методы и программное обеспечение. М.: Мир, 2001.
Моисеев Н.Я. Модифицированный метод расщепления по физическим процессам для решения уравнений радиационной газовой динамики // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 2. С. 303–315.
Гольдин В.Я. Квазидиффузионный метод решения кинетического уравнения // Ж. вычисл. матем. и матем. физ. 1964. Т. 4. № 6. С. 1070–1087.

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

Инструменты

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

Журнал вычислительной математики и математической физики

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