Реализация бикомпактной схемы для HOLO алгоритмов решения уравнения переноса

Аристова Е.Н.; Караваева Н.И.

doi:https://doi.org/10.20948/prepr-2019-21

№ 21 за 2019 г.

https://doi.org/10.20948/prepr-2019-21

Реализация бикомпактной схемы для HOLO алгоритмов решения уравнения переноса

В работе рассмотрено построение бикомпактных схем высокого порядка аппроксимации для HOLO (High Order – Low Order) алгоритма решения нестационарного уравнения переноса в одномерной плоской геометрии. Достигается четвертый прядок аппроксимации по пространству и третий по времени. Рассмотрены два варианта постановки краевых условий для LO части полной системы уравнений: классический – введением дробно-линейных функционалов для отношения потока и плотности излучения, а также по величине плотности излучения из HO части системы. Исследуется эффективность HOLO алгоритма по сравнению с методом итераций источника. Показано, что второй способ реализации краевых условий приводит к третьему порядку сходимости на гладких решениях, в отличие от первого способа, позволяющего получить только второй порядок; однако для первого способа эффективность HOLO алгоритма может падать.

уравнение переноса, метод квазидиффузии, бикомпактная схема, HOLO алгоритмы решения уравнения переноса, потоковая прогонка, диагонально-неявные методы Рунге–Кутты

Implementation of the bicompact scheme for HOLO algorithms for solving the transport equation

Elena Nikolaevna Aristova
Nataliia Igorevna Karavaeva

The construction of bicompact schemes of high order of approximation for unsteady transport equation numerical solving in 1D geometry by application of HOLO (High Order – Low Order) algorithm is considered. The fourth order of approximation in space and the third in time are achieved. Two variants for the boundary conditions for the LO part are considered: the classical one using fractional-linear functionals for the flux and radiation density ratio, and also by the radiation density value from the HO part of the system. The efficiency of the algorithm is demonstrated in comparison with the source iterations method. It is shown the second variant of the boundary conditions leads to the third order of convergence on smooth solutions, in contrast to the first method, which allows to obtain only the second order; however the efficiency of the HOLO algorithm can decrease.

transport equation, quasi-diffusion method, bicompact scheme, HOLO algorithms for transport equation solving, diagonally implicit Runge-Kutta method

Работа выполнена при поддержке Российского фонда фундаментальных исследований, проект 18-01-00857-а.

1. Введение

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

\begin{array}{l} \frac{1}{v} \frac{\partial I}{\partial t} + μ \frac{\partial I}{\partial x} + (κ_{a} (x, t) + κ_{s} (x, t)) I = S (I) + F (x, t, μ), \\ S (I) = κ_{s} (x, t) \int_{- 1}^{1} \sum_{j = 0}^{\infty} \frac{2 j + 1}{2} w_{j} q_{j} (μ') q_{j} (μ) I (x, t, μ') d μ' . \end{array}

(1)

Здесь введены следующие обозначения: $I (x, t, μ)$ – функция распределения, зависящая от фазовых переменных и времени, $μ$ – косинус угла направления полета частицы, составляемого с направлением оси x, S(I) – член рассеяния, $F (x, t, μ)$ – внешний источник, $κ_{a}, κ_{s}$ – коэффициенты поглощения и рассеяния соответственно, $w_{j}$ – коэффициенты разложения индикатрисы рассеяния по полиномам Лежандра.

Для уравнения (1) ставится начально-краевая задача

I (x, 0, μ) = ψ (x, μ),

{I (x, t, μ) |}_{x = x_{0}} = φ_{0} (t, μ), μ \geq 0, {I (x, t, μ) |}_{x = L} = φ_{L} (t, μ), μ \leq 0.

(2)

В ряде задач ставятся более сложные краевые условия отражения.

Уравнение (1) является интегро-дифференциальным, и обычно по интегральным членам рассеяния (и деления) используется итерационный процесс. Для ряда задач такой итерационный процесс может быть медленно сходящимся.

2. HOLO алгоритмы

Эффективным методом ускорения итераций по членам рассеяния (и деления), а также для организации эффективной связи решения уравнения переноса с другими уравнениями, описывающими состояние среды, через которую проходит излучение или поток частиц, являются HOLO алгоритмы. Их основная идея заключается в совместном решении уравнений высокой размерности (HO – high order) и уравнений низкой размерности (LO – low order) [1–8].

В рассматриваемой задаче уравнением высокой размерности является уравнение переноса (1), а уравнением низкой размерности – система уравнений квазидиффузии. Для ее получения введем следующие величины: $U (x, t)$ – плотность излучения, сопоставляемую с интегралом от функции распределения по угловой переменной $μ = \cos θ$ : $U (x, t) = \int_{- 1}^{1} I (x, t, μ) d μ$ , и $W (x, t)$ – поток, сопоставляемый с аналогичным интегралом с весовой функцией $μ$ : $W (x, t) = \int_{- 1}^{1} μ I (x, t, μ) d μ$ . Проинтегрировав уравнение (1) по переменной $μ$ с весом 1 и $μ$ и использовав ортогональность многочленов Лежандра, получим систему уравнений квазидиффузии:

\begin{array}{l} {\begin{cases} \frac{\partial U}{\partial t} + \frac{\partial W}{\partial x} + κ_{0} (x, t) U = Q_{0} (x, t), \\ \frac{\partial W}{\partial t} + \frac{\partial D U}{\partial x} + κ_{1} (x, t) W = Q_{1} (x, t), \end{cases} \\ Q_{0} (x, t) = \int_{- 1}^{1} F (x, t, μ) d μ, Q_{1} (x, t) = \int_{- 1}^{1} μ F (x, t, μ) d μ, \\ κ_{0} = κ_{a}, κ_{1} = κ_{a} + κ_{s} (1 - w_{1}) . \end{array}

(3)

В уравнениях квазидиффузии введен коэффициент квазидиффузии

D (x, t) = \int_{- 1}^{1} μ^{2} I (x, t, μ) d μ / \int_{- 1}^{1} I (x, t, μ) d μ

,

(4)

зависящий дробно-рациональным образом от искомого решения.

В рамках реализации HOLO алгоритма решения уравнения (1) и системы (3) оказываются тесно связанными. Из системы (3) рассчитываются плотность и поток излучения, их значения передаются в уравнение переноса (1) для вычисления главных трех первых членов интеграла рассеяния. Решение уравнения (1) для ряда угловых направлений, соответствующих выбору узлов квадратурной формулы по углам, позволяет рассчитать интегралы, входящие в коэффициент квазидиффузии (4) и в граничные условия для системы уравнений квазидиффузии. В данном простейшем одномерном случае использовались квадратуры Гаусса отдельно для положительных и отрицательных угловых направлений, т.к. на характеристике $μ = 0$ решение не только может быть недифференцируемым, оно может терпеть разрыв.

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

3. Вывод бикомпактной схемы на примере уравнения переноса

Для решения начально-краевой задачи (1)–(2) в [9–15] были разработаны бикомпактные схемы, обладающие четвертым порядком аппроксимации по пространственным переменным и произвольным (в реализации третьим) по времени на минимальном двухточечном шаблоне по каждой из переменных. В качестве примера построения схемы для расчета, реализующей возможности достижения четвертого порядка аппроксимации по пространству и высокого порядка аппроксимации по времени, рассмотрим вывод схемы для уравнения переноса с изотропным рассеянием, которое в схеме HOLO итераций можно считать уравнением с известной правой частью:

\begin{array}{l} \frac{1}{v} \frac{\partial I}{\partial t} + μ \frac{\partial I}{\partial x} + κ_{t} I = Q, \\ κ_{t} = κ_{a} + κ_{s}, Q = 0.5 κ_{s} U + F . \end{array}

(5)

Здесь величины $v$ – скорость рассматриваемой группы частиц и $μ$ – косинус угла, составляемого направлением движения частицы с положительным направлением оси x, они не зависят от пространственной координаты и времени. Коэффициенты поглощения и рассеяния $κ_{a}, κ_{s}$ и их сумма могут зависеть от времени и пространственной координаты.

Для этого уравнения ставится начально-краевая задача (2).

Для начала введем сетку ${x_{i}, 0 \leq i \leq N_{x}}$ на отрезке $x \in [0, L]$ и пространственный шаг $h_{i} = x_{i + 1} - x_{i}$ . В дальнейшем индекс, указывающий на номер ячейки, будет опускаться, так как разностная схема строится в рамках одной ячейки.

Бикомпактная схема строится методом прямых с выводом полудискретной системы уравнений, при котором пространственные производные аппроксимируются конечными разностями, а производные по времени оставляются в дифференциальном виде. Вывод полудискретной формы бикомпактной схемы для уравнения переноса основан на смешанной FV (Finite Volume) – FD (Finite Difference) технике.

Стандартным образом на двухточечном шаблоне возможно достижение второго порядка аппроксимации, для увеличения порядка до четвертого используется расширение списка неизвестных и включение в него помимо узловых значений искомой функции $I_{i}$ интегральных средних от величины $I$ по ячейке. Введем для них следующее обозначение:

{\bar{I}}_{i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} I d x

.

(6)

В рамках FV подхода уравнение (5) интегрируется по ячейке, в результате получается точное следствие исходного уравнения следующего вида:

$\frac{1}{v} \frac{\partial {\bar{I}}_{i}}{\partial t} + \frac{μ}{h} (I_{i + 1} - I_{i}) + \bar{κ_{t, i} I_{i}} = \bar{Q}$ ,

${\bar{Q}}_{i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} Q d x$ .

(7)

В уравнение (7) входит среднее от произведения заданной и искомой функций $\bar{κ_{t, i} I_{i}}$ . Замена $\bar{κ_{t, i} I_{i}} = {\bar{κ}}_{t, i} {\bar{I}}_{i}$ доставляет лишь второй порядок аппроксимации по пространству. При заданном двойном наборе узловых и средних искомых величин для достижения четвертого порядка аппроксимации по пространству воспользуемся формулой Симпсона. Считая, что значение коэффициента $κ_{t}$ доступно в любой пространственной точке, введем его значение в полуцелом узле. Записывая выражение для интегральных средних $\bar{κ_{t, i} I_{i}}$ и $\bar{I_{i}}$ с помощью формулы Симпсона и исключая неизвестной значение в полуцелом узле $I_{i + 1 / 2}$ , имеем:

\begin{array}{l} \bar{κ_{t, i} I_{i}} = (κ_{t, i} I_{i} + 4 κ_{t, i + 1 / 2} I_{i + 1 / 2} + κ_{t, i + 1} I_{i + 1}) / 6, \\ \bar{I_{i}} = (I_{i} + 4 I_{i + 1 / 2} + I_{i + 1}) / 6. \end{array}

(8)

Таким образом, среднее от произведения примет вид:

\bar{κ_{t, i} I_{i}} = (κ_{t, i} - κ_{t, i + 1 / 2}) I_{i} / 6 + κ_{t, i + 1 / 2} \bar{I_{i}} + (κ_{t, i + 1} - κ_{t, i + 1 / 2}) I_{i + 1} / 6

.

(9)

Уравнения (7) с учетом (8) недостаточно для нахождения всех нужных величин. Вторым шагом (FD) для получения замкнутой системы уравнений к уравнению (7) добавляют проинтегрированное по ячейке первое дифференциальное следствие уравнения (5). Первое дифференциальное следствие представляет собой уравнение, правая и левая часть которого суть производные по x от соответствующей части исходного уравнения:

\frac{1}{v} \frac{\partial}{\partial x} \frac{\partial I}{\partial t} + μ \frac{\partial}{\partial x} \frac{\partial I}{\partial x} + \frac{\partial}{\partial x} (κ_{t} I) = \frac{\partial}{\partial x} Q

.

(10)

Тогда при интегрировании (10) по ячейке получим следующее выражение:

\frac{1}{v} \frac{\partial (I_{i + 1} - I_{i})}{\partial t} + μ ({\frac{\partial I}{\partial x} |}_{i + 1} - {\frac{\partial I}{\partial x} |}_{i}) + {(κ_{t} I)}_{i + 1} - {(κ_{t} I)}_{i} = Q_{i + 1} - Q_{i}

.

(11)

В него, помимо узловых значений $I$ , также входят пространственные производные функции распределения в узлах. Для исключения их вхождения в систему уравнений в [5] предложено использовать формулу Эйлера–Маклорена, обладающую четвертым порядком аппроксимации:

\int_{x_{i}}^{x_{i + 1}} g (x) d x = \frac{h}{2} (g_{i} + g_{i + 1}) - \frac{h^{2}}{12} ({g^{'}}_{i + 1} - {g^{'}}_{i}) .

(12)

С учетом связи разности производных с узловыми значениями и интегральным средним (12) уравнение (11) оказывается возможным переписать следующим образом:

\frac{1}{v} \frac{\partial (I_{i + 1} - I_{i})}{\partial t} + \frac{6 μ}{h} (I_{i + 1} - 2 {\bar{I}}_{i} + I_{i}) + {(κ_{t} I)}_{i + 1} - {(κ_{t} I)}_{i} = Q_{i + 1} - Q_{i} .

(13)

Система полудискретных уравнений (7), (13) замкнута относительно переменных $I, \bar{I}$ , не содержит в себе производных по пространству и может быть проинтегрирована по времени любым удобным способом. Мы будем использовать L–устойчивый диагонально-неявный метод Рунге–Кутты третьего порядка аппроксимации с таблицей Бутчера

(14)

Как было показано в [16], каждая стадия данного метода может быть представлена как неявный метод Эйлера. Соответственно, данный метод определяет базовую схему реализации метода Рунге-Кутты по времени. Система уравнений для неявного метода Эйлера в применении к (7), (13) запишется следующим образом:

$\left\{ \begin{array}{l}\frac{1}{v}\frac{{\bar I_i^{n + 1} - \bar I_i^n}}{{\Delta t}} + \left( {\frac{\mu }{h} + \frac{{({\kappa _{t,i + 1}} - {\kappa _{t,i + 1/2}})}}{6}} \right)I_{i + 1}^{n + 1} + \left( {\frac{{({\kappa _{t,i}} - {\kappa _{t,i + 1/2}})}}{6} - \frac{\mu }{h}} \right)I_i^{n + 1} + \\{\rm{ }} + {\kappa _{t,i + 1/2}}\overline {{I_i}} = \bar Q_i^{n + 1},\\\frac{1}{v}\frac{{(I_{i + 1}^{n + 1} - I_i^{n + 1}) - (I_{i + 1}^n - I_i^n)}}{{\Delta t}} + \frac{{6\mu }}{h}(I_{i + 1}^{n + 1} - 2\bar I_i^{n + 1} + I_i^{n + 1}) + \\{\rm{ }} + ({\kappa _t}I)_{i + 1}^{n + 1} - ({\kappa _t}I)_i^{n + 1} = Q_{i + 1}^{n + 1} - Q_i^{n + 1}.\end{array} \right.$

(15)

Система (15) вместе с соответствующим знаку $μ$ граничным условием решается методом бегущего счета на каждом временном слое: для $μ > 0$ слева направо, для $μ < 0$ справа налево. Более подробно численная реализация схемы для уравнения переноса изложена, например, в [16].

4. Схема для системы уравнений квазидиффузии

Дифференциальная система уравнений низкой размерности (уравнений квазидиффузии) получается интегрированием уравнения переноса с весами 1 и $μ$ по угловой переменной. Поэтому схему для LO системы уравнений можно строить двумя различными путями.

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

2. Пусть у нас задана квадратурная формула для вычисления интегралов от функции распределения вида:

\int_{- 1}^{1} f (I, μ) d μ ≃ \sum_{k} c_{k} f (I_{k}, μ_{k})

,

(16)

где $c_{k}$ – веса квадратурной формулы, $μ_{k}$ – узлы квадратуры. Разностная схема может быть получена из разностной схемы для уравнения переноса суммированием с весами квадратурной формулы каждого из двух уравнений и их же, домноженных на $μ_{k}$ .

Выпишем главные этапы вывода схемы первым из этих методов для системы следующего вида (изотропное рассеяние):

{\begin{cases} \frac{\partial U}{\partial t} + \frac{\partial W}{\partial x} + κ_{a} (x, t) U = Q_{0} (x, t), \\ \frac{\partial W}{\partial t} + \frac{\partial D U}{\partial x} + κ_{t} (x, t) W = Q_{1} (x, t) . \end{cases}

(17)

Система уравнений квазидиффузии (LO) (17) получена из уравнения (5) интегрированием по углам с весом 1 и $μ$ . Начальные условия к ней также получаются интегрированием по углам начального условия (2).

В классическом варианте HOLO алгоритмов граничные условия для системы (17) получаются из решения уравнения переноса введением еще одного класса дробно-линейных функционалов:

\begin{array}{l} {\frac{W - W^{i n}}{U - U^{i n}} |}_{x = 0} = c_{0}, {\frac{W - W^{i n}}{U - U^{i n}} |}_{x = L} = c_{L}, \\ c_{0} = \frac{\int_{- 1}^{0} μ I (t, x, μ) d μ}{\int_{- 1}^{0} I (t, x, μ) d μ}, c_{L} = \frac{\int_{0}^{1} μ I (t, x, μ) d μ}{\int_{0}^{1} I (t, x, μ) d μ} . \end{array}

(18)

Значения падающих плотности и потока излучения на границах области получаются из краевых условий (2):

${U^{i n} |}_{x = 0} = \int_{0}^{1} ϕ_{0} (t, μ) d μ,$ ${W^{i n} |}_{x = 0} = \int_{0}^{1} μ ϕ_{0} (t, μ) d μ,$

${U^{i n} |}_{x = L} = \int_{- 1}^{0} ϕ_{L} (t, μ) d μ,$ ${W^{i n} |}_{x = L} = \int_{- 1}^{0} μ ϕ_{L} (t, μ) d μ .$

(19)

В первую очередь для вывода схемы мы расширяем список неизвестных, добавляя средние по ячейке от обеих неизвестных. Введем обозначения:

\begin{array}{l} {\bar{u}}_{i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} U (x, t) d x, {\bar{w}}_{i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} W (x, t) d x, \\ {\bar{Q}}_{0 i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} Q_{0} (x, t) d x, {\bar{Q}}_{1 i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} Q_{1} (x, t) d x . \end{array}

(20)

Проинтегрировав по ячейке оба уравнения, получим первые два уравнения системы

$\left\{ \begin{array}{l}\frac{{\partial {{\bar u}_i}}}{{\partial t}} + \frac{1}{h}({W_{i + 1}} - {W_i}) + \overline {{\kappa _{a,i}}{u_i}} = \overline {{Q_{0i}}} ,\\\frac{{\partial {{\bar w}_i}}}{{\partial t}} + \frac{1}{h}({(DU)_{i + 1}} - {(DU)_i}) + \overline {{\kappa _{t,i}}{w_i}} = \overline {{Q_{1i}}} .\end{array} \right.$

(21)

Далее находим первые дифференциальные следствия от обоих уравнений и тоже интегрируем их по ячейке

$\left\{ \begin{array}{l}\frac{{\partial ({U_{i + 1}} - {U_i})}}{{\partial t}} + \frac{6}{h}({W_{i + 1}} - 2{{\bar w}_i} + {W_i}) + ({({\kappa _a}U)_{i + 1}} - {({\kappa _a}U)_i}) = {Q_{0i + 1}} - {Q_{0i}},\\\frac{{\partial ({W_{i + 1}} - {W_i})}}{{\partial t}} + \frac{6}{h}({(DU)_{i + 1}} - 2\overline {{{(DU)}_i}} + {(DU)_i}) + ({({\kappa _t}W)_{i + 1}} - {({\kappa _t}W)_i}) = \\{\rm{ }} = {Q_{1i + 1}} - {Q_{1i}}.\end{array} \right.$

(22)

Средние $\bar{κ_{a, i} u_{i}}$ и $\bar{κ_{t, i} w_{i}}$ раскроем с четвертым порядком аппроксимации, аналогично формулам (9). В эти формулы входит также среднее от $\bar{D U}$ .

Коэффициент $\bar{D}$ не является независимым коэффициентом, а вводится как дробно-линейный функционал от интегрального среднего в ячейке $\bar{I}$ при интегрировании по углам (4), (16):

\bar{D} = \sum_{k} c_{k} μ_{k}^{2} {\bar{I}}_{k} / \sum_{k} c_{k} {\bar{I}}_{k} .

(23)

Раскрытие средних ${\bar{(D U)}}_{i}$ по формуле, обладающей формально вторым порядком аппроксимации

{\bar{(D U)}}_{i} = {\bar{D}}_{i} {\bar{u}}_{i}

,

из требования кинетической согласованности схем для уравнений квазидиффузии (HO) со схемой для уравнения переноса (LO) приводит к схеме четвертого порядка аппроксимации, т.к в этом случае при сходимости итераций выполняется точное равенство

\bar{D U} = \sum_{k} c_{k} μ_{k}^{2} {\bar{I}}_{k} = \frac{\sum_{k} c_{k} μ_{k}^{2} {\bar{I}}_{k}}{\sum_{k} c_{k} {\bar{I}}_{k}} \sum_{k} c_{k} {\bar{I}}_{k} = \bar{D} \bar{u} .

Интегрирование по угловым переменным для LO уравнений приводит к тому, что на каждом шаге по времени возникает краевая задача по пространству. Для реализации базовой схемы при интегрировании по времени – неявного метода Эйлера, – в данном случае потребуется использовать метод потоковой прогонки. Подробности реализации схемы можно найти в [19]. Для интегрирования по углам используется метод Гаусса с пятью узлами на каждом из полуинтервалов от [–1,0] и [0,1].

5. Аналитический тест для HOLO алгоритма

Составим аналитическое решение для (5) на $x \in [0, 1]$ , в качестве функции распределения выберем

I (x, t, μ) = (1 + μ + e^{- x}) (1 + μ + e^{- t}) .

Коэффициенты рассеяния и поглощения зададим в виде:

\begin{array}{l} κ_{a} = A_{0}, A_{0} = 0.1, \\ κ_{s} = B_{0} + B_{1} \sin x, B_{0} = 0.2, B_{1} = 0.01. \end{array}

Интегрированием по $μ$ получатся следующие точные решения для плотности и потока излучения

\begin{array}{l} U (x, t) = 2 / 3 + 2 e^{- x - t}, \\ W (x, t) = 2 (e^{- x} + e^{- t}) / 3. \end{array}

В вычислениях потребуются величины F, Q₀, Q₁, , найдем F из (5):

F = \frac{1}{v} \frac{\partial I}{\partial t} + μ \frac{\partial I}{\partial x} + κ_{t} I - \frac{κ_{s}}{2} U .

Величины Q₀, Q₁ вычислим по (3).

Начально-краевая задача ставится в соответствии с (2). В [16] показано, что для достижения третьего порядка аппроксимации по времени для методов Рунге–Кутты краевые условия необходимо ставить, используя производные по времени

{{I^{'}}_{t} (x, t, μ) |}_{x = x_{0}} = {φ^{'}}_{t, 0} (t, μ), μ \geq 0, {{I^{'}}_{t} (x, t, μ) |}_{x = L} = {φ^{'}}_{t, L} (t, μ), μ \leq 0,

(24)

и интегрируя краевые условия (24) в общей схеме применения метода Рунге–Кутты. Это связано со структурой погрешности методов Рунге–Кутты [17–18].

Рассмотрим два варианта нахождения решения задачи (5), (2): решение только на основе уравнения переноса, т.е. классическим методом итераций источника (SIM), и в соответствии с HOLO алгоритмом. Результаты для порядков аппроксимации по времени и пространству представлены в таблицах I,II,III,IV,V,VI. Из таблицы III видно, что порядок аппроксимации по времени падает до второго при подключении к решению системы уравнений квазидиффузии. Это связано с тем, что краевые условия (18) мы не можем представить в виде типа условий (24) для производных по времени от искомых функций, которые необходимы для обеспечения третьего порядка сходимости по времени. В представленных ниже таблицах величина $ε$ определяет точность сходимости итераций по членам рассеяния на каждой стадии используемого метода Рунге-Кутты. Исследовалось либо сгущение сеток при фиксированном гиперболическом числе Куранта, и тогда сходимость определяется членами наименьшего порядка аппроксимации, т.е. по времени, либо фиксировался малый временной шаг и измельчались сетки по пространству.

Таблица I

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

$N_{x}$	$Δ t / h = 2, ε = 10^{- 12}$
$N_{x}$	$\lg (\| \| r \| \|_{C} / \| \| u_{e x} \| \|_{C})$	Порядок схемы
25	–6.20	2.95
50	–7.09	2.97
100	–7.99	2.97
200	–8.91	2.89
400	–9.79	2.36
800	–10.54

Graph1 — *Рис. 1*. Зависимость погрешности численного решения от числа узлов по пространству при фиксированном гиперболическом числе Куранта $С = 2$ для аналитического теста (SIM). Для сравнения проведена прямая с наклоном –3 (порядка аппроксимации по времени).

Таблица II

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

$N_{x}$	$Δ t = 10^{- 3}, N_{t} = 25, ε = 10^{- 12}$
$N_{x}$	$\lg (\| \| r \| \|_{C} / \| \| u_{e x} \| \|_{C})$	Порядок схемы
25	–8.51	3.97
50	–9.71	3.58
100	–10.78

Таблица III

Погрешности численного решения и порядок сходимости по времени при использовании HOLO алгоритма для теста с краевыми условиями (18)–(19)

$N_{x}$	$Δ t / h = 2, ε = 10^{- 12}$
$N_{x}$	$\lg (\| \| r \| \|_{C} / \| \| u_{e x} \| \|_{C})$	Порядок схемы
25	–4.74	1.96
50	–5.33	1.93
100	–5.93	1.96
200	–6.53	1.98
400	–7.13	1.99
800	–7.73

Graph2 — *Рис. 2*. Зависимость погрешности численного решения от числа узлов по пространству при фиксированном гиперболическом числе Куранта $С = 2$ для теста c HOLO алгоритмом. Для сравнения проведена прямая с наклоном –2 (порядка аппроксимации по времени).

Таблица IV

Погрешности численного решения и порядок аппроксимации по пространству для HOLO алгоритма при использовании краевых условий (18)–(19) для теста

$N_{x}$	$Δ t = 10^{- 4.35}, N_{t} = 400, ε = 10^{- 12}$
$N_{x}$	$\lg (\| \| r \| \|_{C} / \| \| u_{e x} \| \|_{C})$	Порядок схемы
25	–8.74	3.95
50	–9.94	2.35
100	–10.56

Как видно из представленных данных, подключение системы LO уравнений (уравнений квазидиффузии) понижает порядок сходимости до второго по времени при использовании классической постановки краевых условий в виде (18), сохраняя четвертый порядок сходимости по пространству. Быстрое падение порядка сходимости в Табл. IV связано как с не слишком малым шагом по времени (и вторым порядком аппроксимации по этой переменной), так и с тем, что константа погрешности по пространству мала (1/1512).

Изучим вопрос о том, как добиться третьего порядка аппроксимации по времени при использовании HOLO алгоритма. Для этого рассмотрим еще один вариант постановки граничных условий. В этом случае мы отключим взаимодействие решения уравнения переноса и системы уравнений квазидиффузии через граничные условия (18), ожидая увеличения количества итераций, необходимого для достижения заданной точности по сравнению с HOLO алгоритмом без этой модификации.

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

{U |}_{x = 0} = \int_{- 1}^{1} I (t, 0, μ) d μ

,

{U |}_{x = L} = \int_{- 1}^{1} I (t, L, μ) d μ

.

(25)

Таблица V

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

$N_{x}$	$Δ t / h = 2, ε = 10^{- 12}$
$N_{x}$	$\lg (\| \| r \| \|_{C} / \| \| u_{e x} \| \|_{C})$	Порядок схемы
25	–5.85	3.00
50	–6.75	3.00
100	–7.66	3.00
200	–8.56	2.99
400	–9.46	2.93
800	–10.34

Graph3 — *Рис. 3*. Зависимость погрешности численного решения от числа узлов по пространству при фиксированном гиперболическом числе Куранта $С = 2$ для теста. Для сравнения проведена прямая с наклоном –3 (порядка аппроксимации по времени).

Таблица VI

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

$N_{x}$	$Δ t = 10^{- 4.35}, N_{t} = 400, ε = 10^{- 12}$
$N_{x}$	$\lg (\| \| r \| \|_{C} / \| \| u_{e x} \| \|_{C})$	Порядок схемы
25	–8.74	3.99
50	–9.95	4.03
100	–11.16	2.50
200	–11.92

Сравним эффективность алгоритмов метода итераций источника (SIM – 1) и два варианта HOLO алгоритма: с использованием классических краевых условий метода квазидиффузии с введением дробно-линейных функционалов (18) (HOLOcl – 2) и с постановкой краевых условий непосредственно интегрированием функции распределения (HOLOmod – 2), сравнив число итераций по всем шагам и по всем стадиям, необходимых для достижения заданной точности. Результаты представлены на Рис. 4.

GraphIterationsAnTestnew — *Рис. 4*. Зависимость суммарного по всем стадиям и всем шагам по времени числа итераций, необходимых для достижения точности $ε$ на каждой стадии каждого шага решения. $N_{x} = 25, N_{t} = 25, Δ t = 0.08$ . Закрашенные круги – метод итераций источника (SIM – 1), полузакрашенные треугольники – HOLOmod –2, полые круги – HOLOcl –2.

6. Задача Рида и ее результаты

Применим описанный метод для решения задачи Рида. Задача Рида является стационарной. Мы стационарное решение получим в пределе $t \to \infty$ . В качестве начального условия возьмем ${I (x, t, μ) |}_{t = 0} = 0$ . Справа поставим граничное условие ${I (x, t, μ) |}_{x = L} = 0, μ < 0$ , а слева – условие зеркального отражения ${I (x, t, μ) |}_{x = 0} = {I (x, t, - μ) |}_{x = 0}$ .

Будем решать уравнение (5) совместно с системой (3). Начальные и граничные условия для нее ставятся следующим образом:

${U (x, t) |}_{t = 0} = 0$ , ${W (x, t) |}_{t = 0} = 0$ ,

${W (x, t) |}_{x = 0} = 0$ , ${W (x, t) |}_{x = L} - c_{L} {U (x, t) |}_{x = L} = 0$ .

(26)

Зададим также правые части и коэффициенты $κ_{t}, κ_{s}, κ_{a}$ в соответствии с Таблицей VII.

Таблица VII

Коэффициенты $κ_{t}, κ_{s}, κ_{a}$ и F в задаче Рида

Слой	x	$κ_{t}$	$κ_{s}$	$κ_{a}$	F
1	(0,2)	50	0	50	25
2	(2,3)	5	0	5	0
3	(3,5)	0	0	0	0
4	(5,6)	1	0.9	0.1	0.5
5	(6,8)	1	0.9	0.1	0

Потребуем сходимость итерационного процесса на каждой из трех стадий используемого метода Рунге–Кутты.

Нами было получено стационарное решение (см. рис. 5 и таблица VIII), соответствующее решению этой задачи в [20].

Была использована схема четвертого порядка по пространству, поэтому можно показать, что полученное решение может быть немонотонным. Проведем расчет на нескольких сетках, изобразим на одном графике решения в узловых точках и в полуцелых, посчитав их значения с четвертым порядком точности через узловые значения и средние по ячейкам при помощи формулы Симпсона. Исследуем, при каких параметрах сетки решение является немонотонным и при каких переходит в монотонное (рис. 5,6,7,8,9,10). Рис. 5 иллюстрирует, что момент перехода немонотонного решения в монотонное соответствует $N_{x} = 200$ , и можно сформулировать следующий критерий: если максимальная оптическая толщина слоя, соответствующего одному шагу по пространству, превышает величину (27), то решение немонотонно, в противном случае – монотонно.

\max l = \frac{κ_{a} h}{\min_{k} μ_{k}} = \frac{50 \cdot 4 \cdot 10^{- 2}}{4.69 \cdot 10^{- 2}} \approx 40.

(27)

Таблица VIII

Экстремальные и промежуточные значения плотности излучения в стационарном решении задачи Рида. DD0 – взвешенная «алмазная» схема DD0 с коррекцией отрицательных потоков “fix-up” [22], LN – LN метод [23], МКЭ – метод конечных элементов [20], HOLOcl – реализуемый в данной работе варианта HOLO алгоритма с использованием классических краевых условий метода квазидиффузии с введением дробно-линейных функционалов при $N_{x} = 440$ , $N_{t} = 400$ , $Δ t = 2$ .

	DD0	LN	МКЭ	HOLOcl
Минимальное значение в слое 2	0.02135	0.03229	0.03035	0.02836
Промежуточное значение в слое 3	1.10715	1.11248	1.10711	1.10649
Максимальное значение в слое 4	1.95216	1.94328	1.95050	1.96628

Порядок аппроксимации в области монотонного решения демонстрирует таблица IX. Порядок аппроксимации по времени падает до первого, так как полученное решение не является гладким.

Таблица IX

Погрешности численного решения и порядок сходимости по времени при фиксированном отношении $h / Δ t = 2$ . За точное взято решение при $N_{x} = 1600$

$N_{x}$	$\| \| r \| \|_{C}$	Порядок схемы
400	4.36·10^–2	1.13
800	1.99·10^–2

Эффективность алгоритма демонстрирует рис. 11. Из графика видно, что суммарное число итераций на всех стадиях и шагах по времени в случае, если использовать разработанный HOLO алгоритм, то есть подключать к решению нестационарного уравнения переноса излучения (5) систему уравнений квазидиффузии (3), примерно вдвое меньше, чем без его использования.

GraphIterationsRid — *Рис. 11*. Зависимость суммарного по всем стадиям и всем шагам по времени числа итераций, необходимых для достижения точности $ε$ на каждой стадии каждого шага решения. $N_{x} = 200, N_{t} = 50, C = 2$ . Закрашенные круги – метод итераций источника (SIM – 1), полузакрашенные треугольники – HOLOmod –2, полые круги – HOLOcl –2.

Отличия значений плотности и потока излучения, получаемых из системы уравнений (3) и как сумматорный интеграл из решения уравнения (5) в стационарном решении соответствуют порядку ошибок округления при вычислении с двойной точностью.

Заключение

Были построены и исследованы бикомпактные схемы для HOLO алгоритма решения нестационарного уравнения переноса в одномерной плоской геометрии (5). Схемы обладают третьим порядком аппроксимации по времени и четвертым по пространству. При реализации граничных условий для части LO (18) порядок аппроксимации по времени падает до второго. Порядок аппроксимации может быть вновь повышен для третьего, если в граничных условиях для системы уравнений квазидиффузии (3) использовать значения плотности излучения, рассчитанные из уравнения переноса (25), однако при этом увеличивается число итераций, которые необходимо совершить для достижения заданной точности, так что эффективность алгоритма падает.