1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 62, № 3, 2022

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

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

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

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

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

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

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

Аннотация

Представлены полунеявные и полудискретные разностные схемы повышенной точности для решения кинетических уравнений переноса теплового излучения и уравнения энергии модифицированным методом расщепления. Особенность схем состоит в том, что перенос теплового излучения рассчитывается по явным или неявным разностным схемам, которые аппроксимируют обычное уравнение переноса в частных производных. Взаимодействие излучения с веществом рассчитывается по неявным разностным схемам в полунеявных схемах и из аналитических решений обыкновенных дифференциальных уравнений в полудискретных схемах. Разностные схемы повышенной точности построены на основе метода Рунге–Кутты второго порядка. Решения находятся без внешних итераций по интегралу столкновений и без обращения матриц. Алгоритмы решения разностных уравнений позволяют проводить вычисления в параллельном режиме. Метод естественным образом обобщается на решение задач в многомерных пространствах. Библ. 35. Фиг. 3.

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

1. ВВЕДЕНИЕ

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

Широко используемыми методами решения являются: методы дискретных ординат $\left( {Sn,~DSn~} \right)$ (см. [1], [3], [4]), методы Фотрие (см. [5]) и Райбики (см. [6]), методы характеристик (см. [7], [8]), квазидиффузионные методы (см. [9], [10]), итерационные методы (см. [11]–[15]), методы расщепления (см. [16]–[20]) и т.д. Каждый метод обладает своими достоинствами и недостатками, которые известны и описаны во многих публикациях. Так, если систему уравнений решать методами Фотрие или Райбики, то приходится иметь дело с решением систем линейных алгебраических уравнений с полными матрицами взаимодействия между группами (см. [21]). Как следствие, в случае решения реальных задач требуются большие объемы оперативной памяти и большие временные затраты. Если решение находить итерационными методами по неявным разностным схемам, то основные трудности возникают со сходимостью внешних итераций по интегралу столкновений и независимому источнику в виде функции Планка. Итерации, как правило, сходятся медленно. Для повышения эффективности методов и ускорения сходимости итераций постоянно разрабатываются различные подходы и “ускорители” (см., например, [22]–[25]). Следует отметить работы [26]–[28], в которых оператор переноса рассчитывается по явной разностной схеме, а оператор взаимодействия излучения с веществом – по неявной разностной схеме.

В настоящей работе подробно рассмотрен вопрос построения полунеявных и полудискретных разностных схем повышенных порядков (второго и третьего) аппроксимации для решения интегродифференциальных кинетических уравнений переноса теплового излучения и уравнения энергии в многогрупповом приближении на основе модифицированного метода расщепления (см. [20]). Оператор взаимодействия излучения с веществом рассчитывается по разностной схеме, которая построена с использованием аналитических решений обыкновенных дифференциальных уравнений (ОДУ). Оператор переноса может рассчитываться по известным явным или неявным разностным схемам повышенной точности, например, по схемам из [29], [30]. Решения разностных уравнений находятся без итераций по интегралу столкновений и без обращения матриц.

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

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

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

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

$\frac{1}{с}\frac{{\partial {{I}_{g}}}}{{\partial t}} + \mu \frac{{\partial {{I}_{g}}}}{{\partial x}} + ({{\alpha }_{g}} + {{\beta }_{g}}){{I}_{g}} = 0.5{{\alpha }_{g}}{{B}_{g}} + 0.5{{\beta }_{g}}{{U}_{g}} + 0.5{{f}_{g}},$
(1)
${{U}_{g}} = \mathop \smallint \limits_{ - 1}^1 \,{{I}_{g}}d\mu ,$
$\frac{{dE}}{{dt}} = \mathop \sum \limits_{g = 1}^G \,{{\alpha }_{g}}({{U}_{g}} - {{B}_{g}}).$
Здесь $t,\;x$ – независимые переменные по времени и по пространству, ${{\mu }}$ – косинус угла между вектором движения частиц (фотонов) и осью $x$, $ - 1 \leqslant \mu \leqslant 1,$ с – скорость света, $g$ – индекс группы, ${{I}_{g}}(t,x,\mu )$ – интенсивность энергии излучения в группе $g$, ${{U}_{g}}(t,x)$ – плотность энергии излучения, умноженная на скорость света,
${{B}_{g}} = \frac{{8\pi }}{{{{c}^{2}}{{{\bar {h}}}^{3}}}}\mathop \smallint \limits_{{{\varepsilon }_{g}}}^{{{\varepsilon }_{{g + 1}}}} \frac{{{{\varepsilon }^{3}}}}{{{\text{exp(}}\varepsilon {\text{/}}T{\text{)}} - 1}}d\varepsilon $
есть интенсивность равновесного излучения (функция Планка), умноженная на скорость света, $\bar {h}$ – постоянная Планка, ${{\alpha }_{g}},\;~{{\beta }_{g}}$ – коэффициенты поглощения и рассеяния, ${{f}_{g}}$ – функция источника. Разностная сетка по энергии фотонов включает $G$ независимых групп с энергиями ${{\varepsilon }_{1}} \ldots {{\varepsilon }_{g}} \ldots {{\varepsilon }_{G}}$, ${{\varepsilon }_{g}}$ – энергия фотонов. Разностная сетка по переменной ${{\mu }}$ с центрами ${{\mu }_{m}} = 0.5({{\mu }_{{m + 1/2}}} + {{\mu }_{{m - 1/2}}})$ и шагами $\Delta \mu = {{\mu }_{{m + 1/2}}} - {{\mu }_{{m - 1/2}}}$, $m = 1,~2, \ldots ,M$, включает $M$ направлений движения частиц.

Система уравнений (1) замыкается уравнением состояния вещества в форме $E = E(T)$. Предполагаем, что плотность вещества $\rho = 1$. Начальные условия: $T(0,x) = {{T}^{0}}(x)$, ${{I}_{g}}(0,x,\mu ) = I_{g}^{0}(x,\mu )$.

Граничные условия: на левой границе – ${{I}_{g}}(t,{{x}_{L}},\mu ) = {{\varphi }_{L}}(t,{{x}_{L}},\mu )$ для $\mu > 0$, на правой – ${{I}_{g}}(t,{{x}_{R}},\mu ) = {{\varphi }_{R}}(t,{{x}_{R}},\mu )$ для $\mu < 0$, где ${{\varphi }_{L}}(t,{{x}_{L}},\mu )$, ${{\varphi }_{R}}(t,{{x}_{R}},\mu )$ – известные функции. Требуется найти решение системы уравнений (1) в области $D = \left\{ {{{x}_{L}} \leqslant x \leqslant {{x}_{R}},~~ - 1 \leqslant \mu \leqslant 1} \right\}$ для $t > 0$.

2.2. Полунеявные разностные схемы

Разностную схему для решения системы уравнений (1) будем конструировать на основе Sn-метода из [1] и методов из [16], [17], [20]. В пространстве $(t,x)$ построим разностную сетку с шагами интегрирования $\tau ,\;{{h}_{j}}$ вдоль координатных линий $t,\;x$ соответственно. Границы ячеек разностной сетки с центрами в точках ${{x}_{j}}$ обозначим как ${{x}_{{j - 1/2}}},\;{{x}_{{j + 1/2}}}$, где $j$ – индекс ячейки. Введем обозначения: $E_{j}^{{n + 1}},\;E_{j}^{n}$ – удельные внутренние энергии вещества, $T_{j}^{{n + 1}},\;T_{j}^{n}$ – температуры вещества, $I_{{g,j}}^{{n + 1}},\;I_{{g,j}}^{n}$ – интенсивности излучения, $U_{{g,j}}^{{n + 1}},\;U_{{g,j}}^{n}$ – плотности энергии излучения, $B_{{g,j}}^{{n + 1}},\;B_{{g,j}}^{n}$ – функции Планка, $\alpha _{{g,j}}^{{n + 1}},\;\alpha _{{g,j}}^{n}$ – коэффициенты поглощения, $\beta _{{g,j}}^{{n + 1}},\;\beta _{{g,j}}^{n}$ – коэффициенты рассеяния, $f_{{g,j}}^{{n + 1}},\;f_{{g,j}}^{n}$ – источники. Величины с индексами $n + 1$, $n$ и с целым индексом $j$ являются основными и относятся к центрам ячеек в моменты времени ${{t}^{{n + 1}}} = {{t}^{n}} + \tau $ и ${{t}^{n}}$ соответственно. Величины с дробными индексами ${{I}_{{g,j + 1/2}}},$ ${{I}_{{g,j - 1/2}}},$ которые будем называть “большими”, как в схемах С.К. Годунова (см. [31]), являются вспомогательными и относятся к границам ячеек. В дальнейшем индекс $j$ у величин в центрах ячеек опущен.

Проинтегрировав уравнения (1) по ячейке фазового пространства $(t,x)$, заменив интеграл в (1) квадратурной формулой и применив теорему Гаусса–Остроградского, получим для вычисления основных величин систему разностных уравнений

$I_{g}^{{n + 1}} = \left\{ {I_{g}^{n} - \frac{{c\tau \mu }}{{{{h}_{j}}}}({{I}_{{g,j + 1/2}}} - {{I}_{{g,j - 1/2}}})} \right\} + c\tau [0.5(\alpha _{g}^{*}B_{g}^{*} + \beta _{g}^{*}U_{g}^{*} + f_{g}^{*}) - (\alpha _{g}^{*} + \beta _{g}^{*})I_{g}^{*}],$
(2)
$U_{g}^{{n + 1}} = ~~\mathop \smallint \limits_{ - 1}^1 \,I_{g}^{{n + 1}}d\mu = \mathop \sum \limits_{m = 1}^M \,I_{{g,m}}^{{n + 1}}\Delta \mu ,$
${{E}^{{n + 1}}} = {{E}^{n}} + \tau \mathop \sum \limits_{g = 1}^G \,\alpha _{g}^{*}(U_{g}^{*} - B_{g}^{*}).$
Здесь $I_{g}^{*},\;\alpha _{g}^{*},\;B_{g}^{*},\;U_{g}^{*}$ – средние значения функций на интервале ${{t}^{n}},\;{{t}^{{n + 1}}}$. Эти величины относятся к временному слою ${{t}^{n}} + \tau {\kern 1pt} *$. Введем интенсивности $I_{g}^{{n + 1/2}}$, где $0 \leqslant \tau {\kern 1pt} {\text{*}} \leqslant \tau $, которые вычисляются из выражения в фигурных скобках в (2) по разностным уравнениям
(3)
$I_{g}^{{n + 1/2}} = I_{g}^{n} - \frac{{c\tau \mu }}{{{{h}_{j}}}}({{I}_{{g,j + 1/2}}} - {{I}_{{g,j - 1/2}}}).$
Разностная схема (3) консервативная и аппроксимирует дифференциальное уравнение переноса в частных производных

(4)
$\frac{1}{с}\frac{{\partial {{I}_{g}}}}{{\partial t}} + \mu \frac{{\partial {{I}_{g}}}}{{\partial x}} = 0.$

Если “большие” величины ${{I}_{{g,j + 1/2}}},\;{{I}_{{g,j + 1/2}}}$ вычисляются по величинам с нижнего временного слоя ${{t}^{n}}$, то разностная схема (3) будет явной. Шаг интегрирования по времени выбирается в этом случае из условия устойчивости схемы:

(5)
$1 - \frac{{c\tau \mu }}{{{{h}_{j}}}} > 0 \Rightarrow \tau \leqslant k\frac{{{{h}_{j}}}}{{c\mu }},\quad 0 \leqslant k \leqslant 1.$
Здесь $k$ – коэффициент запаса устойчивости. Условие (5) – это обычное ограничение на шаг интегрирования по времени при решении гиперболических уравнений, которое не зависит от коэффициентов поглощения.

Если интенсивности $I_{g}^{{n + 1/2}}$ вычисляются по величинам с верхнего временного слоя ${{t}^{{n + 1}}}$, то схема (3) будет неявной без ограничения на выбор шага интегрирования по времени. Однако известно, что для обеспечения точности в неявных разностных схемах желательно, чтобы число Куранта было близким к единице.

В дальнейшем предполагаем, что интенсивности $I_{g}^{{n + 1/2}}$ находятся по какой-либо известной разностной схеме повышенной точности, которая аппроксимирует уравнение переноса (4).

Подставив выражение для $I_{g}^{{n + 1}}$ из первого уравнения (2) в квадратурную формулу для $U_{g}^{{n + 1}}$ второго уравнения (2), получим с учетом введенной интенсивности $I_{g}^{{n + 1/2}}$ эквивалентную систему разностных уравнений, которую запишем в виде

$\frac{{I_{g}^{{n + 1}} - I_{g}^{{n + 1/2}}}}{{c\tau }} + (\alpha _{g}^{*} + \beta _{g}^{*})I_{g}^{*} = 0.5\alpha _{g}^{*}B_{g}^{*} + 0.5\beta _{g}^{*}U_{g}^{*} + 0.5f_{g}^{*},$
(6)
$\frac{{U_{g}^{{n + 1}} - U_{g}^{{n + 1/2}}}}{{c\tau }} = - \alpha _{g}^{*}(U_{g}^{*} - B_{g}^{*}) + f_{g}^{*},\quad U_{g}^{{n + 1/2}} = \mathop \sum \limits_{m = 1}^M \,I_{g}^{{n + 1/2}}\Delta \mu ,$
${{E}^{{n + 1}}} = {{E}^{n}} + \frac{1}{с}\mathop \sum \limits_{g = 1}^G \,c\tau \alpha _{g}^{*}(U_{g}^{*} - B_{g}^{*}).$

Особенность разностной схемы (6) в том, что плотность энергии излучения $U_{g}^{{n + 1}}$ вычисляется из разностного уравнения, а не по квадратурной формуле, как в системе уравнений (2). Решение системы уравнений (6) в этом случае не вызывает затруднений и находится без итераций по интегралу столкновений и без обращения матриц.

Первые дифференциальные приближения (ПДП) (см. [32]) разностной схемы (6) запишем в виде

$\frac{1}{с}\frac{{d{{I}_{g}}}}{{dt}} + ({{\alpha }_{g}} + {{\beta }_{g}}){{I}_{g}} = 0.5({{\alpha }_{g}}{{B}_{g}} + {{\beta }_{g}}{{U}_{g}} + {{f}_{g}}) + (\tau {\kern 1pt} {\text{*}} - 0.5\tau )\frac{{d{{F}_{1}}}}{{dt}},$
(7)
$\begin{gathered} \frac{{d{{U}_{g}}}}{{dt}} = - {{\alpha }_{g}}({{U}_{g}} - {{B}_{g}}) + ~{{f}_{g}} + (\tau {\kern 1pt} {\text{*}} - 0.5\tau )\frac{{d{{F}_{2}}}}{{dt}}, \\ \frac{{dE}}{{dt}} = \mathop \sum \limits_{g = 1}^G \,{{\alpha }_{g}}({{U}_{g}} - {{B}_{g}}) + \mathop \sum \limits_{g = 1}^G \,(\tau {\kern 1pt} {\text{*}} - 0.5\tau )\frac{{d{{F}_{2}}}}{{dt}}, \\ \end{gathered} $
${{F}_{1}} = 0.5({{\alpha }_{g}}{{B}_{g}} + {{\beta }_{g}}{{U}_{g}} + {{f}_{g}}) - ({{\alpha }_{g}} + {{\beta }_{g}}){{I}_{g}},\quad {{F}_{2}} = - {\kern 1pt} {{\alpha }_{g}}({{U}_{g}} - {{B}_{g}}) + ~{{f}_{g}}.$
Из ПДП (7) следует, что разностная схема (6) аппроксимирует систему ОДУ
$\frac{1}{с}\frac{{d{{I}_{g}}}}{{dt}} + ({{\alpha }_{g}} + {{\beta }_{g}}){{I}_{g}} = 0.5({{\alpha }_{g}}{{B}_{g}} + {{\beta }_{g}}{{U}_{g}} + {{f}_{g}}),$
(8)
$\frac{1}{с}\frac{{d{{U}_{g}}}}{{dt}} = - {{\alpha }_{g}}({{U}_{g}} - {{B}_{g}}) + ~{{f}_{g}},$
$\frac{{dE}}{{dt}} = \mathop \sum \limits_{g = 1}^G \,{{\alpha }_{g}}({{U}_{g}} - {{B}_{g}})$
с первым порядком по времени. Начальные данные для первых двух уравнений в (8) – это решения уравнения переноса (4).

Погрешность аппроксимации разностной схемы (6) зависит от выбора величин $I_{g}^{*},\;\alpha _{g}^{*},\;B_{g}^{*},\;U_{g}^{*}$. Если $\tau {\kern 1pt} * = \tau $, то величины $I_{g}^{*},\;\alpha _{g}^{*},\;B_{g}^{*},\;U_{g}^{*}$ вычисляются по величинам с верхнего временного слоя $I_{g}^{{n + 1}},\;\alpha _{g}^{{n + 1}},\;B_{g}^{{n + 1}},\;U_{g}^{{n + 1}}$. Разностная схема (6) будет неявной схемой Эйлера, которая аппроксимирует ОДУ (8) с первым порядком по времени с погрешностью $O(\tau )$.

Просуммировав по всем группам второе уравнение в (8) и сложив результат с уравнением энергии в (8), получим для суммарной удельной внутренней энергии вещества и энергии излучения ОДУ

$\frac{d}{{dt}}\left( {Е + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,{{U}_{g}}} \right) = \mathop \sum \limits_{g = 1}^G \,{{f}_{g}},$
первый интеграл которого записывается в виде
(9)
$\left( {{{E}^{{n + 1}}} + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,U_{g}^{{n + 1}}} \right) - \left( {{{E}^{n}} + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,U_{g}^{{n + 1/2}}} \right) = \mathop \sum \limits_{g = 1}^G \mathop \smallint \limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {{f}_{g}}(\tau )d\tau .$
Следовательно, в системе ОДУ (8) изменение суммарной удельной внутренней энергии вещества и энергии излучения за промежуток времени $\tau $ равно энергии, вложенной от внешних источников за время $\tau $. Разностная схема (9) является консервативной. Аналогичное уравнение и вывод можно получить и в разностной схеме (6) из анализа двух последних уравнений в (6). Поэтому температуру ${{T}^{{n + 1}}}$ вещества будем находить из решения уравнения (9) методом Ньютона, как в [20].

Вычитая первое уравнение в (8) из второго, умноженного на 0.5, получим ОДУ

$\frac{1}{с}\frac{{d(0.5{{U}_{g}} - {{I}_{g}})}}{{dt}} = - ({{\alpha }_{g}} - {{\beta }_{g}})(0.5{{U}_{g}} - {{I}_{g}}),$
решение которого выписывается в квадратурах. Это решение позволяет получить для вычисления интенсивности энергии излучения дополнительное уравнение
$I_{g}^{{n + 1}} = 0.5U_{g}^{{n + 1}} - (0.5U_{g}^{{n + 1/2}} - I_{g}^{{n + 1/2}}){\text{exp}}\left( { - c\mathop \smallint \limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} \,({{\alpha }_{g}} - {{\beta }_{g}})d\tau } \right),$
которое можно использовать для контроля точности или как основное уравнение для вычисления интенсивности энергии излучения $I_{g}^{{n + 1}}.$ Здесь плотность энергии излучения $U_{g}^{{n + 1}}$ находится из решения второго уравнения в (6).

Если интенсивности $I_{g}^{{n + 1/2}}$ находятся по явной разностной схеме (3) , а интенсивности $I_{g}^{{n + 1}}$ – по неявной разностной схеме (6) , то разностную схему (3), (6) будем называть полунеявной разностной схемой.

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

2.3. Полунеявные разностные схемы повышенной точности

Известно, что разностные схемы Эйлера первого порядка точности для решения “жестких” ОДУ не эффективны и малопригодны (см. [33]). Поэтому для решения ОДУ желательно применять разностные схемы повышенной точности, например, схемы Рунге–Кутты второго порядка, в которых решения находятся по схеме предиктор-корректор. Рассмотрим построение такой схемы для решения системы ОДУ (8). В дальнейшем, не нарушая общности изложения, предположим, что источники ${{f}_{g}} = 0$. Из ПДП (7) следует, что, если величины $I_{g}^{*},\;\alpha _{g}^{*},\;B_{g}^{*},\;U_{g}^{*}$ вычислять в момент времени $\tau {\kern 1pt} * = 0.5\tau ~$, то разностная схема (6) аппроксимирует ОДУ (8) со вторым порядком по времени с погрешностью $O({{\tau }^{2}})$. Для реализации такой схемы температура $T{\kern 1pt} *$ вещества в общем случае предварительно вычисляется на этапе предиктор в момент времени $t{\kern 1pt} * = {{t}^{n}} + \tau {\kern 1pt} *,\;\tau {\kern 1pt} * \geqslant 0.5\tau $ из уравнения энергии (9), которое запишем в виде

(10)
$E{\kern 1pt} * = {{E}^{n}} + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,(1 - {{\xi }_{1}})(U{{_{g}^{*}}^{{,n + 1/2}}} - B_{g}^{*}),\quad {{\xi }_{1}} = \frac{1}{{1 + c\tau {\kern 1pt} *{\kern 1pt} \alpha _{g}^{*}}}.$
Уравнение $\left( {10} \right)$ получено из уравнения энергии в $(9)$ после подстановки выражения для плотности энергии излучения $U_{g}^{{n + 1}}$ из второго уравнения в $(6)$. Вспомогательные плотности энергии излучения $U{{_{g}^{*}}^{{,n + 1/2}}}$ вычисляются приближенно линейной интерполяцией по формуле
$U{{_{g}^{*}}^{{,n + 1/2}}} = (1 - {{\xi }_{3}})U_{g}^{n} + {{\xi }_{3}}U_{g}^{{n + 1/2}},\quad 0 \leqslant {{\xi }_{3}} \leqslant 1,$
в момент времени $t{\kern 1pt} * = {{t}^{n}} + \tau {\kern 1pt} *$. Здесь ${{\xi }_{3}}$ – параметр, который подбирается, исходя из какого-либо условия, например, ${{\xi }_{3}} = \tau {\kern 1pt} *{\kern 1pt} {\text{/}}\tau ~$. Вычислив температуру вещества $T{\kern 1pt} *$ из уравнения (10), удельные внутренние энергии ${{E}^{{n + 1}}}$ и интенсивности $I_{g}^{{n + 1}}$ найдем на этапе корректор в момент времени ${{t}^{{n + 1}}} = {{t}^{n}} + \tau ~$ из разностных уравнений, которые запишем в виде
${{E}^{{n + 1}}} = {{E}^{n}} + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,(1 - {{\xi }_{1}})(U_{g}^{{n + 1/2}} - B_{g}^{*}),$
(11)
$I_{g}^{{n + 1}} = \frac{{I_{g}^{{n + 1/2}} + 0.5c\tau (\alpha _{g}^{*} + \beta _{g}^{*})((1 - {{\xi }_{2}})B_{g}^{*} + {{\xi }_{2}}U_{g}^{{n + 1/2}})}}{{1 + c\tau (\alpha _{g}^{*} + \beta _{g}^{*})}},$
$~{{\xi }_{1}} = \frac{1}{{1 + c\tau \alpha _{g}^{*}}},\quad {{\xi }_{2}} = \frac{{\beta _{g}^{*}}}{{\alpha _{g}^{*} + \beta _{g}^{*}}}.$
Здесь величины $B_{g}^{*},\;\alpha _{g}^{*},\;\beta _{g}^{*}$ вычисляются по температуре $T{\kern 1pt} *$ в момент времени $t{\kern 1pt} * = {{t}^{n}} + \tau {\kern 1pt} *$. Разностная схема $(10)$ , (11) – это схема предиктор-корректор повышенной точности.

Из практики известно, что, если $\tau {\kern 1pt} * = \tau $, то решения, как правило, монотонные. Это следует и из того, что коэффициенты в разностных схемах положительные. Однако, если $\tau {\kern 1pt} * = 0.5\tau $, то решения могут оказаться немонотонными. В этом случае параметр $\tau {\kern 1pt} *$ полагается равным $\tau {\kern 1pt} * = \xi \tau $, где $0.5 < \xi \leqslant 1,$ и подбирается так, чтобы решение стало монотонным. Это известный прием подавления осцилляций в гибридных разностных схемах (см. [34]).

2.4. Полудискретные разностные схемы

Поскольку разностные уравнения (6) аппроксимируют ОДУ (8), то можно рассмотреть подход, в котором интенсивности $I_{g}^{{n + 1}}$ и удельные внутренние энергии вещества ${{E}^{{n + 1}}}$ вычисляются из аналитических решений ОДУ (8). Такой подход к повышению точности разностных схем отмечается в [3], где утверждается, что в этом случае решения оказываются более точными, чем решения, полученные по разностной схеме. Проинтегрировав по времени ОДУ (8), с учетом того, что первые два уравнения линейные относительно неизвестных ${{I}_{g}},\;{{U}_{g}}$, получим для вычисления основных величин ${{I}_{g}},\;{{U}_{g}},\;E$ следующие уравнения:

${{I}_{g}}(t) = {{\gamma }_{2}}{{I}_{g}}({{t}_{0}}) + 0.5{{\gamma }_{2}}\mathop \smallint \limits_{{{t}_{0}}}^t \,c({{\alpha }_{g}} + {{\beta }_{g}})[(1 - \lambda ){{B}_{g}} + \lambda {{U}_{g}}]\exp (c({{\alpha }_{g}} + {{\beta }_{g}})\tau )d\tau ,$
(12)
$\begin{gathered} {{U}_{g}}(t) = {{\gamma }_{1}}{{U}_{g}}({{t}_{0}}) + {{\gamma }_{1}}\mathop \smallint \limits_{{{t}_{0}}}^t \,c{{\alpha }_{g}}{{B}_{g}}{\text{exp}}(c{{\alpha }_{g}}\tau )d\tau , \\ E(t) = E({{t}_{0}}) + \mathop \sum \limits_{g = 1}^G \,\mathop \smallint \limits_{{{t}_{0}}}^t \,{{\alpha }_{g}}({{U}_{g}} - {{B}_{g}})d\tau , \\ \end{gathered} $
${{\gamma }_{1}} = \exp ( - ct{{\alpha }_{g}}),\quad {{\gamma }_{2}} = \exp ( - c({{\alpha }_{g}} + {{\beta }_{g}})t),\quad \lambda = \frac{{{{\beta }_{g}}}}{{{{\alpha }_{g}} + {{\beta }_{g}}}}.$
В общем случае интегралы вычисляются приближенно по квадратурным формулам. В частном случае, если коэффициенты поглощения, рассеяния и функция Планка постоянные, то интегралы в (12) берутся в квадратурах, которые запишем в виде
${{I}_{g}}(t) = {{\gamma }_{2}}I_{g}^{{n + 1/2}} + 0.5({{\gamma }_{1}} - {{\gamma }_{2}})U_{g}^{{n + 1/2}} + 0.5(1 - {{\gamma }_{1}}){{B}_{g}},$
(13)
$\begin{gathered} {{U}_{g}}(t) = {{\gamma }_{1}}U_{g}^{{n + 1/2}} + (1 - {{\gamma }_{1}}){{B}_{g}}, \\ E(t) = E({{t}_{0}}) + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,(1 - {{\gamma }_{1}})(U_{g}^{{n + 1/2}} - {{B}_{g}}), \\ \end{gathered} $
${{\gamma }_{1}} = \exp ( - ct{{\alpha }_{g}}),\quad {{\gamma }_{2}} = \exp ( - ct({{\alpha }_{g}} + {{\beta }_{g}})).$
Аналитические решения (13) будем использовать при построении разностных схем. Пусть $t \equiv {{t}^{{n + 1}}},\;{{t}_{0}} \equiv {{t}^{n}}$. Зафиксируем температуру $T = T{\kern 1pt} *$ в момент времени $t{\kern 1pt} * = {{t}^{n}} + \tau {\kern 1pt} *$, $0 \leqslant \tau {\kern 1pt} * \leqslant \tau $. При построении разностных схем предполагается, что коэффициенты поглощения, рассеяния и функция Планка усредняются по фиксированной температуре $T{\kern 1pt} *$ и остаются постоянными на интервале $[{{t}^{n}},{{t}^{{n + 1}}}]$. Следовательно, численные решения можно найти из уравнений (13) в момент времени ${{t}^{{n + 1}}}$, которые записываются в виде

${{E}^{{n + 1}}} = {{E}^{n}} + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,(1 - {{\gamma }_{1}})(U_{g}^{{n + 1/2}} - B_{g}^{*}),$
(14)
$I_{g}^{{n + 1}} = {{\gamma }_{2}}I_{g}^{{n + 1/2}} + 0.5({{\gamma }_{1}} - {{\gamma }_{2}})U_{g}^{{n + 1/2}} + 0.5(1 - {{\gamma }_{1}})B_{g}^{*},$
${{\gamma }_{1}} = \exp ( - c\tau \alpha _{g}^{*}),\quad {{\gamma }_{2}} = \exp ( - c\tau (\alpha _{g}^{*} + \beta _{g}^{*})).$

Если $\tau {\kern 1pt} * = \tau $, то температура $T{\kern 1pt} * = {{T}^{{n + 1}}}$ находится из уравнения энергии (14) в момент времени ${{t}^{{n + 1}}} = {{t}^{n}} + \tau $. Величины $\alpha _{g}^{*},\;\beta _{g}^{*},\;B_{g}^{*}$ вычисляются по температуре $~{{T}^{{n + 1}}}$. Интегралы в (12) вычисляются по квадратурным формулам методом правых прямоугольников с погрешностью $O(\tau )$.

Если $\tau {\kern 1pt} * = 0.5\tau $, то интегралы в (12) вычисляются по квадратурным формулам методом средних с погрешностью $O({{\tau }^{2}})$. Для реализации такой схемы температура $T{\kern 1pt} * = {{T}^{{n + 1/2}}}$ вещества предварительно вычисляется на этапе предиктор в момент времени ${{t}^{{n + 1}}} = {{t}^{n}} + 0.5\tau $ из уравнения энергии в (14), которое записывается в виде

(15)
$\begin{gathered} {{E}^{{n + 1/2}}} = {{E}^{n}} + \frac{1}{c}\mathop \sum \limits_{g = 1}^G \,(1 - {{\gamma }_{1}})(U{{_{g}^{*}}^{{n + 1/2}}} - B_{g}^{{n + 1/2}}), \\ {{\gamma }_{1}} = \exp ( - c\tau {\kern 1pt} *{\kern 1pt} \alpha _{g}^{{n + 1/2}}),\quad \tau {\kern 1pt} * = 0.5\tau . \\ \end{gathered} $
Величины $\alpha _{g}^{*},\;\beta _{g}^{*},\;B_{g}^{*}$ вычисляются по температуре $T{\kern 1pt} * = {{T}^{{n + 1/2}}}$. Удельные внутренние энергии ${{E}^{{n + 1}}}$ и интенсивности $I_{g}^{{n + 1}}$ находятся из явных разностных уравнений (14) на этапе корректор в момент времени ${{t}^{{n + 1}}} = {{t}^{n}} + \tau $. В результате получаем разностную схему (15), (14) типа предиктор-корректор повышенной точности.

Если для вычисления интегралов в (12) применить двухточечную формулу Гаусса (см. [33]), которая имеет погрешность $O({{\tau }^{3}})$, то на этапе предиктор необходимо вычислить два набора вспомогательных температур $T_{1}^{{n + 1/2}},\;T_{2}^{{n + 1/2}}$ из уравнения (15) в моменты времени

$t_{1}^{{n + 1/2}} = {{t}^{n}} + 0.5\tau - \frac{\tau }{{2\sqrt 3 }},\quad t_{2}^{{n + 1/2}} = {{t}^{n}} + 0.5\tau + \frac{\tau }{{2\sqrt 3 }}$
соответственно. Подынтегральная функция в этом случае усредняется по формуле

$F = 0.5(F(T_{2}^{{n + 1/2}}) + F(T_{1}^{{n + 1/2}})).$

Поскольку разностные схемы (14), (15) получены из аналитических решений (14), то разностные схемы (3), (14), (15) будем называть полудискретными разностными схемами повышенной точности.

При вычислении интегралов в (12) приближенно по квадратурным формулам коэффициенты поглощения, рассеяния и функция Планка полагаются равными усредненным значениям на интервале интегрирования $[{{t}^{n}},{{t}^{{n + 1}}}]$. Очевидно, что, если шаг интегрирования по времени $\tau = {{t}^{{n + 1}}} - {{t}^{n}}$ большой и функции резко изменяются, то усреднение может оказаться грубым. Следовательно, численное решение будет найдено с большой погрешностью. Погрешность можно уменьшить, если интегралы в (12) вычислять по обобщенной квадратурной формуле с шагом интегрирования $\tilde {\tau } < \tau $, разбив интервал интегрирования $[{{t}^{n}},{{t}^{{n + 1}}}]$ на $K$ подынтервалов. В этом случае интегралы в (12) представляются в виде суммы интегралов, которые вычисляются на этих подынтервалах. На этапе предиктор дополнительно вычисляются температуры ${{T}^{{k + 1/2}}}$ в точках $\tau _{k}^{*} = k\tilde {\tau },$ $\tau _{k}^{*} = (k - 0.5)\tilde {\tau }$, $k = 1,~2,...,~K$, в схемах первого и второго порядков точности соответственно.

Для обеспечения гарантированной точности в схемах (9), (14) вычисления в каждой точке желательно проводить с контролем точности. Например, иметь критерий точности или проводить два расчета с шагами интегрирования $\tau $ и $0.5\tau $. Если различие результатов расчетов не превосходит заданной величины, то решение получено с заданной точностью. В противном случае шаг уменьшается в 2 раза и расчет повторяется. В обоих случаях такая технология счета обеспечит гарантированную точность и эффективность расчетов.

3. ВЕРИФИКАЦИЯ ПРОГРАММЫ И МЕТОДИКИ

Работоспособность представленной методики и программы проверялась путем сравнения результатов расчетов модельной задачи на сходимость с результатами, которые получены по методике из [15] на разностной сетке с шагом интегрирования по пространству $h = 0.000125$. Расчеты на сходимость проведены на последовательности сгущающихся разностных сетках с шагами интегрирования по пространству ${{h}_{{i + 1}}} = 0.5{{h}_{i}}$, $~{{h}_{0}} = 0.002$, $~i = 0,~1, \ldots ,~5$, которым соответствует число интервалов 100, 200, 400 и т.д. Шаги интегрирования по времени выбирались из условия устойчивости (5) с коэффициентами запаса $k = 1$, $k = 0.25$. В расчетах перенос излучения рассчитывался по явной разностной схеме повышенной точности из [29]. Коэффициент интерполяции ${{\xi }_{3}} = 0.0625$.

3.1. Постановка модельной задачи

На левую границу $(x = 0)$ плоского слоя толщиной 0.2 см падает планковский поток излучения, соответствующий температуре вещества Т = 10 кэВ. На правой границе $(x = 0.2)$ задается условие свободной поверхности ${{I}_{g}}(t) = 0~$ для $\mu > 0$. Слой состоит из двух физических областей. Коэффициент поглощения взят из [35] и вычисляется в первой области по формуле

${{\alpha }_{g}} = \left\{ \begin{gathered} \frac{{27(1 - {{e}^{{ - {{\varepsilon }_{g}}/T}}})}}{{\varepsilon _{g}^{3}}},\quad 0 \leqslant x \leqslant 0.1,\quad ~{{\varepsilon }_{g}} \leqslant 30, \hfill \\ 10~000.0,\quad 0 \leqslant x \leqslant 0.1,\quad {{\varepsilon }_{g}} > 30, \hfill \\ \end{gathered} \right.$
во второй области –
${{\alpha }_{g}} = \frac{{0.001(1 - {{e}^{{ - {{\varepsilon }_{g}}/T}}})}}{{\varepsilon _{g}^{3}}},\quad 0.1 \leqslant x \leqslant 0.2.$
Коэффициент рассеяния ${{\beta }_{g}} = 0$. Начальная температура вещества Т = 0.001 кэВ. Уравнение состояния вещества Е = 0.81Т. По энергии фотонов взято 15-групповое приближение из [9] с границами интервалов, умноженными на 10: [0, 3, 6, 8, 12, 15, 18, 24, 27, 30, 40, 50, 70, 90, 110, 150]. Средние значения энергии в группах относятся к центрам интервалов и равны полусумме значений энергий на границах интервалов. Надо найти состояние системы в момент времени $t = 0.002$.

3.2. Результаты расчетов по полунеявным и полудискретым разностным схемам

Представим результаты расчетов по полунеявным и полудискретным разностным схемам предиктор-корректор повышенной точности в момент времени $t = 0.002$. Температура вещества ${{T}^{{n + 1/2}}}$ вычисляется в момент времени ${{t}^{{n + 1/2}}} = {{t}^{n}} + {{\tau }^{{n + 1/2}}}$, где ${{\tau }^{{n + 1/2}}} = 0.78\tau $. Параметр ${{\tau }^{{n + 1/2}}}$ подобран так, чтобы решения уравнения энергии оставались монотонными. Результаты расчетов по полунеявным и полудискретным разностным схемам близки между собой. Поэтому на графиках приведены результаты расчетов по полудискретным разностным схемам. Разностные схемы повышенной точности обозначим символом $С_{2}^{3}$, а на фиг. 1–3 – С32. Это означает, что разностная схема аппроксимирует уравнение переноса (4) на гладких решениях с третьим порядком по времени и по пространству, а уравнение энергии со вторым порядком по времени. На фиг. 1, 2 методика из [15] обозначена аббревиатурой “era”.

Фиг. 1.

Зависимости от $x$ температур вещества в расчетах на сходимость: (а) – $k = 1$, (б) – $k = 0.25$.

Фиг. 2.

Зависимости от x температур излучения в расчетах на сходимость: (а) – $k = 1$, (б) – $k = 0.25$.

Фиг. 3.

Зависимости от $\varepsilon $ спектральных потоков излучения: квадраты – входящий поток, кружки и треугольники – выходящий.

Зависимости от $x$ температур вещества в расчетах задачи на сходимость по полудискретной разностной схеме $С_{2}^{3}$ приведены на фиг. 1 в момент времени $t = 0.002$.

Зависимости от $x$ температур излучения в расчетах задачи на сходимость по полудискретной разностной схеме $С_{2}^{3}$ приведены на фиг. 2 в момент времени $t = 0.002$.

Результаты расчетов входящих и выходящих спектральных потоков излучения по полудискретной схеме $С_{2}^{3}$ представлены на фиг. 3 в момент времени $t = 0.002$. Графики потоков, рассчитанные с коэффициентами запаса $k = 1$, $k = 0.25,$ практически совпадают между собой. Поэтому приведены результаты расчетов с коэффициентом запаса $k = 1$.

Из фиг. 1–3 следует, что результаты расчетов температуры вещества, температуры излучения и выходящего спектрального потока излучения по полудискретной разностной схеме (3), (14) удовлетворительно согласуются с результатами расчетов по методике из [15]. Времена расчетов в последовательном режиме по обеим методикам сравнимы между собой.

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

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

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

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

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

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

  1. Карлсон Б., Белл Дж. Решение транспортного уравнения Sn-методом. Физика ядерных реакторов. М.: Атомиздат, 1959. С. 408–432.

  2. Бай Ши-и. Динамика излучающего газа. М.: Мир, 1968.

  3. Четверушкин Б.Н. Математическое моделирование задач динамики излучающего газа. М.: Наука, 1985.

  4. Басс Л.П., Волощенко А.М., Гермогенова Т.А. Методы дискретных ординат в задачах о переносе излучения. М.: ИПМ им. М.В. Келдыша РАН, 1986.

  5. Featrier P.C.R. Sur la resolution numerique de I’equation de transfert // Acad. Sci. Paris. 1964. V. 258. P. 3198–3210.

  6. Rybicki G. A modified Feautrier method // J. of Quantitative Spectroscopy and Radiative Transfer. 1971. V. 11. P. 589–596.

  7. Владимиров В.С. Численное решение уравнения для сферы // Вычисл. матем. М.: ВЦ АН СССР, 1958. С. 3–33.

  8. Никифорова А.В., Тарасов В.А., Трощиев В.Е. О решении кинетических уравнений дивергентным методом характеристик // Ж. вычисл. матем. и матем. физ. 1972. Т. 12. № 4. С. 1041–1048.

  9. Гольдин В.Я. Квазидиффузионный метод решения кинетического уравнения // Ж. вычисл. матем. и матем. физ. 1964. Т. 4. № 6. С. 1070–1087.

  10. Карлыханов Н.Г. Применение метода квазидиффузии для решения задач переноса излучения // Вопросы атомной науки и техники. Сер. Матем. моделирование физ. процессов. 2010. Вып. 1. С. 32–38.

  11. Гусев В.Ю., Козманов М.Ю., Рачилов Е.Б. Метод решения неявных разностных уравнений, аппроксимирующих системы уравнений переноса и диффузии излучения // Ж. вычисл. матем. и матем. физ. 1984. Т. 24. № 12. С. 1842–1849.

  12. Долголёва Г.В. Численное решение системы уравнений, описывающей перенос излучения и взаимодействие его с веществом // Вопросы атомной науки и техники. Сер. Матем. моделирование физ. процессов. 1991. Вып. 1. С. 58–60.

  13. Федотова Л.П., Шагалиев Р.М. Конечно-разностный метод КМ-метод для двумерных нестационарных процессов переноса в многогрупповом кинетическом приближении // Матем. моделирование. 1991. Т. 3. № 6. С. 29–41.

  14. Анистратов Д.Ю., Аристова Е.Н., Гольдин В.Я. Нелинейный метод решения задач переноса излучения в среде // Матем. моделирование. 1996. Т. 8. № 12. С. 3–28.

  15. Карлыханов Н.Г. Построение оптимальных многодиагональных методов решения задач переноса излучения // Ж. вычисл. матем. и матем. физ. 1997. Т. 37. № 4. С. 494–498.

  16. Марчук Г.И., Яненко Н.Н. Решение многомерного кинетического уравнения методом расщепления // Докл. АН СССР. 1964. Т. 157. № 6. С. 1291–1292.

  17. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

  18. Марчук Г.И. Методы расщепления. М.: Наука, 1988.

  19. Самарский А.А., Вабищев П.Н. Аддитивные схемы для задач математической физики. М.: Наука, 2001.

  20. Моисеев Н.Я. Модифицированный метод расщепления по физическим процессам для решения уравнений радиационной газовой динамики // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 2. С. 303–315.

  21. Михалс Д. Звездные атмосферы. М.: Мир, 1982.

  22. Alcouff R.E. A stable diffusion synthetic acceleration method for neutron transport iterations // Trans. Am. Nucl. Soc. 1976. V. 23. P. 203.

  23. Alcouff R.E., McCoy D.R., Larsen E.W. Finite difference effects in the synthetic acceleration method // Trans. Am. Nucl. Soc. 1981. V. 39. P. 462.

  24. Гаджиев А.Д., Селезнёв В.Н., Шестаков A.A. DS-n метод с искусственной диссипацией и ВДМ–метод ускорения итераций для численного решения двумерного уравнения переноса теплового излучения в кинетической модели // Вопросы атомной науки и техники. Сер. Матем. моделирование физ. процессов. 2003. Вып. 4. С. 33–46.

  25. Завьялов В.В., Шестаков А.А. Выделение диагонального элемента для ускорения итераций в кинетическом приближении при расчете теплопереноса // Матем. моделирование. 2010. Т. 22. № 2. С. 93–104.

  26. Klar A. An asymtotic-induced scheme for nonstationary transport equations in the diffusive limit // SIAM J. 1998. Numer. Anal. V. 35. № 3. P. 1073–1078.

  27. Klar A., Unterreiter A. Uniform stability of a finite difference scheme for transport equations in diffusive regimes // SIAM J. 2002. Numer. Anal. V. 40. № 3. P. 891–913.

  28. McClarren R.G., Evans T.M., Lowrie R.B., Densmore J.D. Semi-implicit time integration for ${{P}_{N}}$ thermal radiative transfer // J. of Comp. Phys. 2008. V. 227. Iss. 16. P. 7561–7586.

  29. Моисеев Н.Я., Силантьева И.Ю. Разностные схемы произвольного порядка аппроксимации для решения линейных уравнений переноса с постоянными коэффициентами методом Годунова с антидиффузией // Ж. вычисл. матем. и матем. физ. 2008. Т. 48. № 7. С. 1282–1293.

  30. Моисеев Н.Я. Неявные разностные схемы бегущего счета повышенной точности // Ж. вычисл. матем. и матем. физ. 2011. Т. 51. № 5. С. 920–935.

  31. Годунов С.К., Забродин А.В., Прокопов Г.П., Иванов М.Я. Численное решение многомерных задач газовой динамики. Под ред. С.К. Годунова. М.: Наука, 1976.

  32. Шокин Ю.И. Метод дифференциального приближения. Новосибирск: Наука, 1979.

  33. Каханер Д., Моулер К., Неш С. Численные методы и программное обеспечение. М.: Мир, 2001.

  34. Федоренко Р.П. Применение разностных схем высокой точности для решения гиперболических уравнений // Ж. вычисл. матем. и матем. физ. 1962. Т. 2. № 6. С. 1122–1128.

  35. Fleck J.F., Cummings J.D., Jr., An implicit Monte-Carlo scheme for calculating time and frequency dependent radiation transport // J. Comp. Phys. 1971. V. 8. № 3. P. 313–342.

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

Инструменты

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

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

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