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

Караваева Н.И.

doi:https://doi.org/10.20948/prepr-2023-25

№ 25 за 2023 г.

https://doi.org/10.20948/prepr-2023-25

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

Караваева Н.И.

Предложены и реализованы бикомпактные разностные схемы высокого порядка аппроксимации для одного из этапов HOLO алгоритма для решения уравнения переноса излучения, представляющего собой решение системы уравнений квазидиффузии совместно с уравнением энергии. Схемы обладают четвертым порядком аппроксимации по пространству на двухточечном шаблоне, интегрирование по времени проводится при помощи диагонально-неявного метода Рунге–Кутты третьего порядка аппроксимации, каждая стадия которого может быть представлена как неявный метод Эйлера. Для нахождения решения предложено организовать итерационный процесс, так как коэффициенты системы являются нелинейными функциями температуры. Схемы применены к ряду аналитических тестов, в которых показана сходимость по времени и пространству с заявленными порядками.

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

Bicompact schemes for the joint solution of a one-group system of quasi-diffusion equations with the energy equation

Nataliia Igorevna Karavaeva

High-order bicompact schemes are proposed for one of the stages of the HOLO algorithm for solving the radiative transport equation and implemented. The system of quasi-diffusion equations is solved together with the energy equation. The schemes have a fourth-order approximation in space on a two-point stencil, integration over time is carried out using the diagonal-implicit Runge–Kutta method of the third order of approximation, each stage of the method can be represented as an implicit Euler method. It is proposed to organize an iterative process to find a solution, since the coefficients of the system are non-linear functions of temperature. The schemes are used for a number of analytical tests, the convergence in time and space with the declared orders is shown.

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

1. Введение

Линейное уравнение переноса описывает процесс переноса излучения или незаряженных частиц. К необходимости решения уравнения переноса приводят различные задачи науки и техники от моделирования конструкции и поведения активных зон быстрых реакторов до расчета сложных газодинамических течений, сопровождающихся активным переносом излучения, как, например, в задачах УТС или при входе возвращаемых космических аппаратов в атмосферу Земли.

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

\frac{1}{с} \frac{\partial I_{ν}}{\partial t} + μ \frac{\partial I_{ν}}{\partial x} + {κ^{'}}_{ν} I_{ν} = {κ^{'}}_{ν} I_{ν}^{P l}

,

(1)

\frac{\partial ε}{\partial t} = \int_{0}^{\infty} \int_{- 1}^{1} {κ^{'}}_{ν} (I_{ν} - I_{ν}^{P l}) d ν d μ

.

(2)

Уравнение переноса (1) и уравнение энергии (2) дополняются уравнением состояния

ε = ε (T)

,

(3)

а также начальными и граничными условиями

$I_{ν} (x, 0, μ, ν) = ψ (x, μ, ν)$ ,

${I_{ν} (x, 0, μ, ν) |}_{x = 0} = ϕ_{0} (t, μ, ν)$ при $μ > 0$ ,

${I_{ν} (x, 0, μ, ν) |}_{x = L} = ϕ_{L} (t, μ, ν)$ при $μ < 0$ ,

$T (x, 0) = T_{i n} (x)$ ,

(4)

где $I_{ν} (x, t, μ, ν)$ – спектральная интенсивность излучения, $ε (T)$ – внутренняя энергия вещества, $T (x, t)$ – температура вещества, ${κ^{'}}_{ν} (T, ν)$ – спектральный коэффициент поглощения излучения с поправкой на вынужденное испускание, $I_{ν}^{P l} (T, ν)$ – спектральная интенсивность равновесного излучения, с – скорость света, $μ$ – косинус угла направления полета частицы, составляемого с направлением оси x.

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

Задача (1-4) может быть решена в многогрупповом приближении, тогда уравнения (1-2), начальные и граничные условия (4) примут вид

\frac{1}{с} \frac{\partial I^{p}}{\partial t} + μ \frac{\partial I^{p}}{\partial x} + κ_{U}^{p} I^{p} = κ_{P l}^{p} I_{P l}^{p}

,

(5)

\frac{\partial ε}{\partial t} = \int_{- 1}^{1} \sum_{p = 1}^{Ρ} (κ_{U}^{p} I^{p} - κ_{P l}^{p} I_{P l}^{p}) d μ

.

(6)

$I^{p} (x, 0, μ) = ψ^{p} (x, μ)$ , $p = 1, ..., Ρ$ ,

${I^{p} (x, 0, μ) |}_{x = 0} = ϕ_{0}^{p} (t, μ)$ при $μ > 0$ ,

${I^{p} (x, 0, μ) |}_{x = L} = ϕ_{L}^{p} (t, μ)$ при $μ < 0$ ,

(7)

T (x, 0) = T_{i n} (x)

,

(8)

где p – номер группы,

I^{p} (x, t, μ) = \int_{ν_{p}}^{ν_{p + 1}} I_{ν} (x, t, μ, ν) d ν

,

I_{P l}^{p} (T) = \int_{ν_{p}}^{ν_{p + 1}} I_{ν}^{P l} (T, ν) d ν

,

ψ^{p} (x, μ) = \int_{ν_{p}}^{ν_{p + 1}} ψ (x, μ, ν) d ν

,

ϕ_{0}^{p} (t, μ) = \int_{ν_{p}}^{ν_{p + 1}} ϕ_{0} (t, μ, ν) d ν

,

ϕ_{L}^{p} (t, μ) = \int_{ν_{p}}^{ν_{p + 1}} ϕ_{L} (t, μ, ν) d ν

,

κ_{P l}^{p} (T) = \int_{ν_{p}}^{ν_{p + 1}} {κ^{'}}_{ν} (T, ν) I_{ν}^{P l} (T, ν) d ν / \int_{ν_{p}}^{ν_{p + 1}} I_{ν}^{P l} (T, ν) d ν

,

κ_{U}^{p} (T, θ^{p}) = \int_{ν_{p}}^{ν_{p + 1}} {κ^{'}}_{ν} (T, ν) I_{ν}^{P l} (θ^{p}, ν) d ν / \int_{ν_{p}}^{ν_{p + 1}} I_{ν}^{P l} (θ^{p}, ν) d ν

,

$I^{p} (x, t, μ)$ – групповая интенсивность излучения, $θ^{p}$ – эффективная температура излучения [3,4], $κ_{U}^{p}$ , $κ_{P l}^{p}$ – коэффициенты поглощения в p-й группе.

Одним из подходов для решения задачи является использование HOLO алгоритмов, которые позволяют организовать эффективную связь решения уравнения переноса с параметрами среды [5-13]. Одним из вариантов HOLO алгоритма является метод квазидиффузии В.Я. Гольдина [5,7].

Применение HOLO алгоритма состоит в совместном решении уравнения высокой размерности (HO – high order) и низкой размерности (LO – low order). При решении уравнения переноса излучения (уравнение HO) совместно с уравнением энергии размерность задачи понижается сначала за счет осреднения по угловым координатам. Вводятся групповые плотность $U^{p} (x, t) = \int_{- 1}^{1} I^{p} (x, t, μ) d μ$ и поток излучения $W^{p} (x, t) = \int_{- 1}^{1} μ I^{p} (x, t, μ) d μ$ , уравнение (5) интегрируют по $μ$ с весом 1 и $μ$ . В результате получается многогрупповая система уравнений квазидиффузии (уравнения LO 1), в которую входит групповой коэффициент квазидиффузии

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

,

замыкающий систему уравнений [5]. Далее проводят осреднение по энергии и получают одногрупповую систему уравнений квазидиффузии (уравнения LO 2), из которой могут быть найдены одногрупповая плотность и поток излучения. Постановка задачи будет приведена в следующем разделе.

В данной работе предложены бикомпактные схемы третьего порядка по времени и четвертого порядка аппроксимации по пространству для одного из этапов HOLO алгоритма, в котором решается одногрупповая система уравнений квазидиффузии совместно с уравнением энергии. Бикомпактные схемы для уравнений HO были предложены Б.В. Роговым [14-20]. Разностные схемы для уравнений LO 1 разработаны авторами [21,22] на основе схем Б.В. Рогова. Объединение решения всех этапов HOLO алгоритма [23] позволит применить предложенные разностные схемы для блока решения уравнения переноса совместно с уравнением баланса энергии в задачах радиационной газовой динамики.

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

Эффективная одногрупповая система уравнений квазидиффузии в одномерной плоской геометрии имеет вид [23]

{\begin{cases} \frac{1}{c} \frac{\partial U}{\partial t} + \frac{\partial W}{\partial x} + κ_{U} U = Q_{P l}, \\ \frac{1}{c} \frac{\partial W}{\partial t} + \frac{\partial (D U)}{\partial x} + κ_{R} W + ζ U = Q_{W}, \end{cases}

(9)

Q_{W} = 0

,

U = \sum_{p = 1}^{Ρ} U^{p}

,

W = \sum_{p = 1}^{Ρ} W^{p}

,

κ_{U} = \sum_{p = 1}^{Ρ} κ_{U}^{p} U^{p} / \sum_{p = 1}^{Ρ} U^{p}

,

D = \sum_{p = 1}^{Ρ} D^{p} U^{p} / \sum_{p = 1}^{Ρ} U^{p}

,

κ_{R} = \sum_{p = 1}^{Ρ} κ_{R}^{p} | W^{p} | / \sum_{p = 1}^{Ρ} | W^{p} |

,

ζ = \sum_{p = 1}^{Ρ} (κ_{R}^{p} - κ_{R}) W^{p} / \sum_{p = 1}^{Ρ} U^{p}

,

где U – плотность излучения, W – поток излучения, D – коэффициент квазидиффузии, $κ_{R}^{p}$ – коэффициент поглощения, осредненный по Росселанду, в p-й группе, $ζ$ – функционал, введенный для корректировки осреднения $κ_{R}^{p}$ со знакопеременной весовой функцией $W^{p}$ .

Начальные условия для системы (9), как и величины $W^{i n}$ , $U^{i n}$ , $c_{0}$ , $c_{L}$ , входящие в граничные условия (10), получаются тем же способом, что и сама система – осреднением начальных условий (7) по угловой переменной с весом 1 и $μ$ и последующим суммированием по группам. Классический способ постановки краевых условий в HOLO алгоритмах состоит во введении дробно-линейных функционалов

{\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}

.

(10)

Система (9) решается совместно с уравнением энергии (6), которое записывается в виде

\frac{\partial ε}{\partial t} = κ_{U} U - Q_{P l}

(11)

и дополняется начальным условием (8). Для замыкания системы уравнений (9,11) она дополняется уравнением состояния (3) в виде [23]

ε = a T

.

(12)

Коэффициенты поглощения $κ_{U}$ , $κ_{R}$ и правая часть $Q_{P l}$ являются нелинейными функциями температуры. Уравнение состояния (3) также может представлять собой нелинейную функцию температуры.

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

Для численного решения задачи предлагается использовать бикомпактные схемы, обладающие высоким временным и пространственным порядком аппроксимации. Схемы строятся методом прямых аналогично схемам Б.В. Рогова [14-20] для уравнения переноса.

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

{\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

.

(13)

Проинтегрируем каждое уравнение системы (9) по ячейке и получим

$\left\{ \begin{array}{l}\frac{1}{c}\frac{{\partial {{\bar u}_i}}}{{\partial t}} + \frac{1}{h}({W_{i + 1}} - {W_i}) + \overline {{\kappa _U}{U_i}} = {{\bar Q}_{Pl}}_i,\\\frac{1}{c}\frac{{\partial {{\bar w}_i}}}{{\partial t}} + \frac{1}{h}({(DU)_{i + 1}} - {(DU)_i}) + {\overline {{\kappa _R}W} _i} + \overline {\zeta {U_i}} = {{\bar Q}_W}_i,\end{array} \right.$

(14)

где ${\bar{Q}}_{P l}_{i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} Q_{P l} d x,$ ${\bar{Q}}_{W}_{i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} Q_{W} d x$ .

Для замыкания системы (14) необходимо также проинтегрировать по ячейке дифференциальные следствия уравнений (9). Если исходная система была записана в виде

L_{1} (U, W) \equiv \frac{1}{c} \frac{\partial U}{\partial t} + \frac{\partial W}{\partial x} + κ_{U} U = Q_{P l}

,

L_{2} (U, W) \equiv \frac{1}{c} \frac{\partial W}{\partial t} + \frac{\partial (D U)}{\partial x} + κ_{R} W + ζ U = Q_{W}

,

то первыми дифференциальными следствиями будут уравнения вида

\frac{\partial L_{1} (U, W)}{\partial x} = \frac{\partial Q_{P l}}{\partial x}

,

\frac{\partial L_{2} (U, W)}{\partial x} = \frac{\partial Q_{W}}{\partial x}

.

Интегрируя по ячейке первые дифференциальные следствия и заменяя разность пространственных производных по формуле Эйлера–Маклорена

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

,

получим еще два уравнения, замыкающие систему (14):

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

(15)

В систему (14) входит среднее от произведений $\bar{κ_{U} U_{i}}$ , ${\bar{κ_{R} W}}_{i}$ , $\bar{ζ U_{i}}$ . Для завершения построения разностной схемы предложим способ аппроксимировать выражение $\bar{ζ U_{i}}$ с четвертым порядком по пространству. Воспользуемся формулой Симпсона для $\bar{ζ U_{i}}$ и ${\bar{u}}_{i}$ , получим

\bar{ζ U_{i}} = (ζ_{i + 1} U_{i + 1} + 4 ζ_{i + 1 / 2} U_{i + 1 / 2} + ζ_{i} U_{i}) / 6

,

{\bar{u}}_{i} = (U_{i + 1} + 4 U_{i + 1 / 2} + U_{i}) / 6

и далее исключим $U_{i + 1 / 2}$ из формулы для $\bar{ζ U_{i}}$ :

\bar{ζ U_{i}} = (ζ_{i} - ζ_{i + 1 / 2}) U_{i} / 6 + (ζ_{i + 1} - ζ_{i + 1 / 2}) U_{i + 1} / 6 + ζ_{i + 1 / 2} {\bar{u}}_{i}

.

Аналогично аппроксимируем выражения $\bar{κ_{U} U_{i}}$ , ${\bar{κ_{R} W}}_{i}$ .

Среднее от произведения ${\bar{(D U)}}_{i}$ , которое входит в систему (15), при решении модельной задачи аппроксимируем аналогично. Однако при решении задачи (5-8) необходимо будет использовать более простую аппроксимацию $\bar{(D U)} = \bar{D} \cdot \bar{u}$ , которая в итоге приведет к схеме четвертого порядка аппроксимации [22] в силу кинетической согласованности разностных схем для HO, LO 1 и LO 2 уравнений.

Полученная система (14, 15) не содержит производных по пространству и может быть проинтегрирована по времени различными методами. В данной работе используется L–устойчивый диагонально-неявный метод Рунге–Кутты третьего порядка аппроксимации с таблицей Бутчера:

.

(16)

Преимуществом в реализации метода Рунге-Кутты с такой таблицей Бутчера является то, что каждая стадия этого метода может быть представлена как неявный метод Эйлера [25]. Для решения получаемой системы уравнений на каждом шаге по времени используется эффективный метод потоковой прогонки. Подробности реализации схемы можно найти в [21].

Как было показано в [25], использование граничных условий вида (10) позволяет достичь только второго порядка сходимости по времени вместо третьего для метода (16). В [25] также предлагается способ постановки граничных условий без понижения порядка сходимости метода по времени. Граничные условия задаются в виде

$μ_{0} (t) U (0, t) + τ_{0} (t) W (0, t) = η_{0} (t)$ ,

$μ_{L} (t) U (L, t) + τ_{L} (t) W (L, t) = η_{L} (t)$ ,

(17)

а коэффициенты выбираются следующим образом:

$μ_{0} = 1$ , $τ_{0} = 0$ , $η_{0} = U (0, t) + Δ t {U^{'}}_{e x} (0, t + Δ t)$ ,

$μ_{L} = 1$ , $τ_{L} = 0$ , $η_{L} = U (L, t) + Δ t {U^{'}}_{e x} (L, t + Δ t)$ ,

(18)

где ${U^{'}}_{e x}$ – производная точного решения для плотности излучения.

4. Разностная схема для уравнения энергии. Итерационный процесс

Уравнение баланса энергии (11) интегрируется по времени при помощи того же метода Рунге–Кутты (16), что и одногрупповая система уравнений квазидиффузии (9). Для нахождения узловых значений температуры $T_{i}$ и средних значений температуры по ячейке

{\bar{T}}_{i} = \frac{1}{h} \int_{x_{i}}^{x_{i + 1}} T (x, t) d x

на новом слое по времени необходимо организовать итерационный процесс на каждой из стадий метода Рунге–Кутты, которые могут быть сведены к использованию неявного метода Эйлера в виде

a \frac{{\overset{⌢}{T}}_{i} - T_{i}}{Δ t} = {\overset{⌢}{κ}}_{U}_{i} {\overset{⌢}{U}}_{i} - {\overset{⌢}{Q}}_{P l}_{i}

,

(19)

a \frac{{\overset{⌢}{\bar{T}}}_{i} - {\bar{T}}_{i}}{Δ t} = {\overset{⌢}{\bar{κ}}}_{U}_{i} {\overset{⌢}{\bar{U}}}_{i} - {\overset{⌢}{\bar{Q}}}_{P l}_{i}

,

(20)

где ${\overset{⌢}{T}}_{i} = T_{i}^{n + 1}$ , ${\overset{⌢}{κ}}_{U}_{i} = κ_{U} (T_{i}^{n + 1})$ , ${\overset{⌢}{U}}_{i} = U_{i}^{n + 1}$ , ${\overset{⌢}{Q}}_{P l}_{i} = Q_{P l} (T_{i}^{n + 1})$ и ${\overset{⌢}{\bar{T}}}_{i} = {\bar{T}}_{i}^{n + 1}$ , ${\overset{⌢}{\bar{κ}}}_{U}_{i} = κ_{U} ({\bar{T}}_{i}^{n + 1})$ , ${\overset{⌢}{\bar{U}}}_{i} = {\bar{U}}_{i}^{n + 1}$ , ${\overset{⌢}{\bar{Q}}}_{P l}_{i} = Q_{P l} ({\bar{T}}_{i}^{n + 1})$ , а также учтено (12).

Рассмотрим способ линеаризации величин, входящих в правую часть уравнений (19, 20), на примере уравнения (19). Будем опускать пространственный индекс, так как он одинаков для всех величин, входящих в (19, 20). Пусть s – номер итерации, тогда

{\overset{⌢}{κ}}_{U}^{(s)} = {\overset{⌢}{κ}}_{U}^{(s - 1)} + {\frac{δ κ_{U}}{δ T} |}_{T = {\overset{⌢}{T}}^{(s - 1)}} \cdot ({\overset{⌢}{T}}^{(s)} - {\overset{⌢}{T}}^{(s - 1)})

,

{\overset{⌢}{Q}}_{P l}^{(s)} = {\overset{⌢}{Q}}_{P l}^{(s - 1)} + {\frac{d Q_{P l}}{d T} |}_{T = {\overset{⌢}{T}}^{(s - 1)}} \cdot ({\overset{⌢}{T}}^{(s)} - {\overset{⌢}{T}}^{(s - 1)})

.

Уравнение (19) примет вид

a \frac{{\overset{⌢}{T}}_{i}^{(s)} - T_{i}}{Δ t} = {({\overset{⌢}{κ}}_{U} \overset{⌢}{U} - {\overset{⌢}{Q}}_{P l})}_{i}^{(s - 1)} + {(\frac{δ κ_{U}}{δ T} \overset{⌢}{U} - \frac{d Q_{P l}}{d T})}_{i}^{(s - 1)} ({\overset{⌢}{T}}_{i}^{(s)} - {\overset{⌢}{T}}_{i}^{(s - 1)})

,

(21)

из которого могут быть получены значения температуры на новом слое по времени на новой итерации в каждой точке сетки по пространству при известных значениях ${\overset{⌢}{U}}_{i}^{(s - 1)}$ . Полученные значения ${\overset{⌢}{T}}_{i}^{(s)}$ будут использоваться для нахождения ${\overset{⌢}{U}}_{i}^{(s)}$ , ${\overset{⌢}{W}}_{i}^{(s)}$ из разностной схемы для одногрупповой системы уравнений квазидиффузии, которая примет вид

{\begin{cases} \frac{1}{c} \frac{\partial U^{(s)}}{\partial t} + \frac{\partial W^{(s)}}{\partial x} + κ_{U}^{(s)} U^{(s)} = Q_{P l}^{(s)}, \\ \frac{1}{c} \frac{\partial W^{(s)}}{\partial t} + \frac{\partial (D^{(s)} U^{(s)})}{\partial x} + κ_{R}^{(s)} W^{(s)} + ζ^{(s)} U^{(s)} = Q_{W}^{(s)} . \end{cases}

(22)

Коэффициенты системы и правые части рассчитываются по ${\overset{⌢}{T}}_{i}^{(s)}$ , ${\overset{⌢}{\bar{T}}}_{i}^{(s)}$ , найденным из (19, 20). Итерационный процесс, состоящий из поочередного нахождения поля температур ${\overset{⌢}{T}}_{i}^{(s)}$ , ${\overset{⌢}{\bar{T}}}_{i}^{(s)}$ из (19, 20) и нахождения плотности ${\overset{⌢}{U}}_{i}^{(s)}$ и потока излучения ${\overset{⌢}{W}}_{i}^{(s)}$ из (22), повторяется до сходимости с заданной точностью $ε$ . На каждом слое по времени на нулевой итерации значения температуры ${\overset{⌢}{T}}_{i}^{(0)}$ и плотности излучения ${\overset{⌢}{U}}_{i}^{(0)}$ приравняем к значениям с предыдущего слоя по времени.

5. Аналитические тесты

Пусть $x \in [0, 1]$ , $c = 2$ , $a = 1$ , $ε = 10^{- 13}$ , предложим несколько тестов для демонстрации порядков сходимости по времени и пространству.

Тест 1. Плотность и поток излучения зададим формулами $U (x, t) = A_{1} t^{3} + A_{2} x^{3}$ и $W (x, t) = B_{1} t^{3} + B_{2} x^{3}$ соответственно, а коэффициент поглощения $κ_{U} = T^{2}$ . Коэффициенты $κ_{R}$ и $ζ$ в тесте 1 выберем не зависящими от температуры, так как их линеаризация в решении не проводится. Пусть $κ_{R} = 0.3 + 0.1 \sin x$ , $ζ = 1 + \sin x$ . Зависимость коэффициента квазидиффузии $D$ от пространственной переменной выберем в виде $D = 1 / 3 + 0.01 {(x - 0.5)}^{2}$ . Правые части системы уравнений (9) получим по формулам

Q_{P l} = \frac{1}{c} \frac{\partial U}{\partial t} + \frac{\partial W}{\partial x} + κ_{U} U

,

Q_{W} = \frac{1}{c} \frac{\partial W}{\partial t} + \frac{\partial (D U)}{\partial x} + κ_{R} W + ζ U

.

(23)

Граничные условия поставим в соответствии с (18). Табл. I демонстрирует достижение третьего порядка сходимости по времени и четвертого по пространству для решения, полученного для плотности излучения $U$ и потока излучения $W$ . Сходимость третьего порядка по времени полученного решения для температуры $T$ продемонстрируем при помощи процесса Эйткена (см. Табл. II).

Таблица I

Погрешности численного решения и порядок сходимости по времени и пространству для теста 1 для плотности $U$ и потока излучения $W$

Сходимость по времени при фиксированном $Δ t / h = 1.04,$ $t \in [0, 0.12],$ $A_{1} = B_{1} = 1,$ $A_{2} = B_{2} = 0$ .				Сходимость по пространству при $Δ t = 7 \cdot 10^{- 4},$ $t \in [0, 0.035],$ $A_{1} = B_{1} = 0,$ $A_{2} = B_{2} = 1$ .
$N_{x}$	$N_{t}$	$\| \| {[u]}_{τ h} - y \| \|_{C}$	p	$N_{x}$	$\| \| {[u]}_{τ h} - y \| \|_{C}$	p
520	60	6.25·10^-10	2.97	640	7.86·10^-08	4.29
1040	120	7.97·10^-11	2.97	1280	4.03·10^-09	3.98
2080	240	1.02·10^-11		2560	2.56·10^-10

Таблица II

Порядок сходимости по времени для теста 1 для температуры $T$ при фиксированном $Δ t / h = 1.04,$ $t \in [0, 0.12],$ $A_{1} = B_{1} = 1,$ $A_{2} = B_{2} = 0$ .

$N_{x}$	$N_{t}$	$\| \| y_{h} - y_{h / 2} \| \|_{C}$	p
520	60	9.31·10^-11	3.00
1040	120	1.16·10^-11

Тест 2. Пусть $U = \exp (x + t),$ $W = \exp (x + t) / 3,$ $κ_{U} = T^{2},$ $κ_{R} = T^{2},$ $D = 1 / 3$ , а $ζ = 1 + \sin x,$ $T_{i n} = 1$ . Правые части уравнений получим по формулам (23). Зависимость температуры от x и t может быть получена аналитически. Уравнение энергии (11) примет вид

a \frac{\partial T}{\partial t} = - \frac{1}{c} \exp (x + t) - \frac{1}{3} \exp (x + t)

(24)

и может быть проинтегрировано. Пусть $α = - \frac{c + 3}{3 a c}$ , уравнение (24) можно преобразовать к виду $\frac{\partial T}{\partial t} = α \exp (x + t)$ . Тогда $T = T_{i n} + α \exp (x) \cdot (\exp (t) - 1)$ .

Сходимость по времени третьего порядка и сходимость по пространству четвертого порядка при использовании граничных условий (18) показывают Табл. III и Рис. 1.

Таблица III

Порядок сходимости по времени и по пространству для теста 2

Сходимость по времени при фиксированном $Δ t / h = 0.25,$ $t \in [0, 0.25]$ .				Сходимость по пространству при $Δ t = 8 \cdot 10^{- 5},$ $N_{t} = 10$ .
$N_{x}$	$N_{t}$	$\| \| {[u]}_{τ h} - y \| \|_{C}$	p	$N_{x}$	$\| \| {[u]}_{τ h} - y \| \|_{C}$	p
50	50	6.09·10^-08	2.91	30	8.84·10^-09	3.95
100	100	8.08·10^-09	2.95	60	5.71·10^-10	4.15
200	200	1.05·10^-09	2.96	120	3.22·10^-11	3.96
400	400	1.34·10^-10	2.96	240	2.07·10^-12
800	800	1.72·10^-11

C:\МОИ (Н. К. И.)\МФТИ\Аристова\СТАТЬИ\Статья 9 КД+энергия препринт\Graph4_.jpg — А)

C:\МОИ (Н. К. И.)\МФТИ\Аристова\СТАТЬИ\Статья 9 КД+энергия препринт\Graph5_.jpg — А)

6. Заключение

Для численного решения уравнения переноса излучения совместно с уравнениями, описывающими состояние среды, могут применяться HOLO алгоритмы. Одним из этапов HOLO алгоритма будет совместное решение одногрупповой системы уравнений квазидиффузии (9) совместно с уравнением энергии (11). Для одногрупповой системы уравнений квазидиффузии (9) были предложены бикомпактные разностные схемы на минимальном двухточечном шаблоне по пространству. За счет расширения списка неизвестных и добавления в него интегральных средних по ячейке достигается четвертый порядок аппроксимации по пространству. Так как схемы строятся методом прямых, то интегрирование по времени можно производить при помощи различных методов. В данной работе используется L–устойчивый диагонально-неявный метод Рунге–Кутты третьего порядка аппроксимации, каждая стадия которого может быть сведена к использованию неявного метода Эйлера. Уравнение энергии интегрируется по времени при помощи того же метода Рунге–Кутты. Так как коэффициенты системы (9) и правые части нелинейно зависят от температуры, то для решения предлагается организовать итерационный процесс. Производится линеаризация коэффициента поглощения $κ_{U}$ и правой части $Q_{P l}$ . Бикомпактные схемы для HOLO алгоритма применены к набору аналитических тестов, в которых продемонстрированы третий порядок сходимости по времени и четвертый по пространству. Предложенные разностные схемы можно будет применить для полного HOLO алгоритма для решения задачи (1-4), являющейся одним из блоков метода расщепления при решении задач радиационной газовой динамики.