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

1. Введение

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

Численное решение уравнения переноса излучения сопряжено с рядом трудностей, в числе которых высокая размерность задачи, нелинейность и нелокальность взаимодействия излучения с веществом, широкий диапазон, в котором изменяются коэффициенты поглощения [3]. Для эффективного понижения размерности используются HOLO-алгоритмы, в которых организуется совместное решение уравнений высокой размерности (HO – high order) и низкой размерности (LO – low order) [4-12]. Нелинейность и жесткость задачи обуславливают использование неявных схем. Значительное изменение коэффициента поглощения от ячейки к ячейке ведет к необходимости использовать разностные схемы на компактном шаблоне.

Перечисленным требованиям отвечают бикомпактные схемы, предложенные Б.В. Роговым и его коллегами для решения уравнения переноса [13-19], а также разработанные авторами на их основе бикомпактные схемы для уравнений более низкой размерности [20, 21]. Схемы строятся методом прямых на минимальном двухточечном шаблоне, обладают четвёртым порядком аппроксимации по пространству. Для интегрирования по времени может быть использован метод любого порядка аппроксимации. Используемый авторами диагонально-неявный метод Рунге–Кутты третьего порядка аппроксимации хорошо зарекомендовал себя при решении модельной задачи переноса нейтронов [22] и поэтому в рамках данной работы применяется для интегрирования по времени. Бикомпактные схемы для HOLO-алгоритма решения задач переноса излучения были ранее исследованы на аналитических тестах [23].

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

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

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

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

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

Искомыми являются следующие функции: $I_{\nu }(x,t,\mu ,\nu )$ – спектральная интенсивность излучения, проинтегрированная по азимутальному углу, $T(x,t)$ – температура вещества. В (1-4) используются следующие обозначения: $\epsilon (T)$ – внутренняя энергия вещества, ${\kappa }_{\nu }^{\prime }(T,\nu )$ – спектральный коэффициент поглощения излучения с поправкой на вынужденное испускание, $I_{\nu }^{Pl}(T,\nu )$ – спектральная интенсивность равновесного излучения, проинтегрированная по азимутальному углу, с – скорость света, $\mu$ – косинус угла направления полета частицы, составляемого с направлением оси x.

Спектральная интенсивность равновесного излучения, проинтегрированная по азимутальному углу, представляет собой функцию Планка, умноженную на $2\pi$:

где h – постоянная Планка, k – постоянная Больцмана.

3. Многогрупповое приближение

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

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

$I^p(x,t,\mu )$ – групповая интенсивность излучения, ${\theta }^p$ – эффективная температура излучения [26], ${\kappa }_U^p$, ${\kappa }_{Pl}^p$ – коэффициенты поглощения в p-й группе.

Правую часть уравнения (5) для удобства обозначим $Q^p$ и преобразуем к виду

где ${\sigma }^p(T)=\frac{2h\pi }{c^2{T}^4}{\displaystyle \underset{{\nu }_p}{\overset{{\nu }_{p+1}}{\int }}{\nu }^3{\left(e^{h\nu /kT}-1\right)}^{-1}d\nu }$.

4. HOLO-алгоритм

Для понижения размерности задачи (3, 5-8) и получения многогрупповой системы уравнений квазидиффузии (уравнений LO₁) проводится осреднение по угловой переменной. Вводятся групповые плотность $U^p(x,t)={\displaystyle \underset{-1}{\overset{1}{\int }}{I}^p(x,t,\mu )d\mu }$ и поток излучения $W^p(x,t)={\displaystyle \underset{-1}{\overset{1}{\int }}\mu I^p(x,t,\mu )d\mu }$, уравнение (5) интегрируется по $\mu$ с весами 1 и $\mu$, и получается система

где $Q_1^p={\displaystyle \underset{-1}{\overset{1}{\int }}{Q}^{p}d\mu }=2Q^p=4{\kappa }_{Pl}^p{\sigma }^p{T}^4={\kappa }_{Pl}^p{U}_{Pl}^p$, $U_{Pl}^p=4{\sigma }^p{T}^4$, ${\kappa }_R^p$ – коэффициент поглощения, осреднённый по Росселанду [26], в p-й группе с эффективной температурой ${\theta }^p$, $D^p(x,t)$ – групповой коэффициент квазидиффузии, замыкающий систему уравнений [4],

Начальные условия для системы (9) также получаются интегрированием условий (7) по $\mu$ с весами 1 и $\mu$:

Краевые условия в классическом варианте постановки задачи представляются в виде дробно-линейных функционалов

где коэффициенты $c_0^p$, $c_L^p$ рассчитываются по выходящей из области функции распределения

Следующим этапом осреднения является осреднение по энергии, в результате чего получается одногрупповая система уравнений квазидиффузии (уравнения LO₂). Введём одногрупповые плотность излучения $U={\displaystyle \underset{p=1}{\overset{{\rm P}}{\sum }}{U}^p}$ и поток излучения $W={\displaystyle \underset{p=1}{\overset{{\rm P}}{\sum }}{W}^p}$. Просуммируем уравнения (9) по группам и получим

где D – коэффициент квазидиффузии, $\zeta$ – функционал, введённый для корректировки осреднения ${\kappa }_R^p$ с возможно знакопеременной весовой функцией $W^p$, $\sigma$ – постоянная Стефана–Больцмана. Начальные условия для (13) также получаются суммированием групповых значений плотности и потока излучения

Аналогично поступают с граничными условиями. Введём

Тогда граничные условия примут вид

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

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

Коэффициенты поглощения ${\kappa }_U$, ${\kappa }_R$ и правая часть $Q_{Pl}$ являются нелинейными функциями температуры. Уравнение состояния (3) также может имеет другой вид и представлять собой нелинейную функцию температуры.

В рамках HOLO-алгоритма осуществляется эффективное взаимодействие между уравнениями высокого порядка HO (5) и уравнениями низкого порядка LO₁ (9) и LO₂ (13). Решая уравнение переноса (HO) (5), получают групповые коэффициенты граничных условий (12) и групповой коэффициент квазидиффузии (10) для каждой группы для дальнейшего решения многогрупповой системы уравнений квазидиффузии LO₁ (9). После решения системы LO₁ (9) производят суммирование по группам и получают коэффициенты граничных условий (15) и коэффициент квазидиффузии (14), необходимые для решения одногрупповой системы уравнений квазидиффузии LO₂ (13) совместно с уравнением энергии (17). В результате совместного решения уравнений (13) и (17) получают значения температуры, через которую пересчитываются значения коэффициентов и правых частей, входящих в (5,9).

5. Разностная схема для HOLO-алгоритма

Для HOLO-алгоритма предлагается использовать разностную схему, ранее описанную в [23]. Для каждого из этапов HOLO-алгоритма, на которых решаются уравнения HO (5), LO₁ (9) и LO₂ (13), используется бикомпактная схема, построенная на двухточечном шаблоне методом прямых. Порядок аппроксимации по пространству равен четырём, а порядок аппроксимации по времени определяется используемым численным методом интегрирования по времени. В данной работе будет использоваться диагонально-неявный метод Рунге–Кутты третьего порядка аппроксимации с таблицей Бутчера

Каждая стадия метода (19) может быть реализована как неявный метод Эйлера [27].

Построение пространственной аппроксимации состоит из нескольких этапов. Список неизвестных расширяют, и помимо значений искомой функции в узлах в него включают интегральные средние по ячейке. Уравнения и их дифференциальные следствия при дифференцировании по пространственной переменной интегрируют по ячейке. Для замыкания получаемой системы используют формулу Эйлера–Маклорена, что и определяет четвёртый порядок аппроксимации по пространству. Полученная система обыкновенных дифференциальных уравнений интегрируется по времени при помощи метода (19).

Для уравнения HO (5) система уравнений, полученная методом прямых, принимает вид

где ${\bar{I}}_i^p=\frac{1}{h}{\displaystyle \underset{x_i}{\overset{x_{i+1}}{\int }}{I}^{p}dx}$, ${\bar{Q}}_i^p=\frac{1}{h}{\displaystyle \underset{x_i}{\overset{x_{i+1}}{\int }}{Q}^{p}dx}$. Для среднего по ячейке от произведения $\overline{{{\kappa }_U^p}_i{I}_i^p}$ используется формула четвёртого порядка аппроксимации [28]

Интегрирование по времени системы (20) при помощи метода (19) приводит к расчёту неизвестных на новом временном слое в режиме бегущего счёта.

Система дифференциальных уравнений для многогрупповой системы уравнений квазидиффузии LO₁ (9) записывается в виде

где ${\bar{U}}_i^p=\frac{1}{h}{\displaystyle \underset{x_i}{\overset{x_{i+1}}{\int }}{U}^{p}dx}$, ${\bar{W}}_i^p=\frac{1}{h}{\displaystyle \underset{x_i}{\overset{x_{i+1}}{\int }}{W}^{p}dx}$, ${\bar{Q}}_{1,i}^p=\frac{1}{h}{\displaystyle \underset{x_i}{\overset{x_{i+1}}{\int }}{Q}_1^{p}dx}$. Средние от произведений, входящие в первые два уравнения (22), аппроксимируются аналогично (21). Среднее $\overline{D_i^p{U}_i^p}$ аппроксимируется формально со вторым порядком: $\overline{D_i^p{U}_i^p}={\bar{D}}_i^p\cdot {\bar{U}}_i^p$, но ввиду кинетической согласованности схем и того, что коэффициент $D_i^p$ не является независимым коэффициентом, это в итоге не приводит к понижению порядка сходимости [28].

Использование метода (19) для интегрирования по времени системы (22) приводит к системе линейных алгебраических уравнений, которая может быть решена методом потоковой прогонки.

Аналогичный вид имеет система дифференциально-разностных уравнений, получаемая методом прямых для системы (13)

Уравнение баланса энергии интегрируется по времени также методом (19).

Бикомпактные схемы для уравнения переноса были разработаны Б.В. Роговым и его коллегами [13-19], бикомпактные схемы для многогрупповой и одногрупповой систем квазидиффузии были разработаны авторами [20, 21]. Так как схемы построены в рамках одной ячейки, можно использовать неравномерную сетку по пространству, а также менять шаг по времени в процессе расчёта в зависимости от того, сколько итераций необходимо для сходимости итерационных процессов.

6. Итерационный процесс

Используемая неявная разностная схема для HOLO-алгоритма, а также нелинейный вид зависимости коэффициентов системы и правых частей от температуры определяют необходимость организовать итерационный процесс для нахождения решения задачи. В [23, 25] подробно описаны все этапы итерационного процесса, в рамках этой работы приведём его короткое описание и схему (рис. 1).

Расчёт начинается с внутренних квазиньютоновских итераций, в процессе которых определяется следующее приближение к температуре на новом слое по времени с точностью ${\epsilon }_T$. Для этого решается одногрупповая система уравнений квазидиффузии LO₂ (13) совместно с уравнением энергии (17). По новому приближению для температуры рассчитываются коэффициенты и правые части в уравнении переноса HO (5), после чего находится новое приближение для спектральной интенсивности излучения для различных угловых направлений в каждой из групп по энергии. По полученному решению многогруппового уравнения переноса находятся коэффициенты граничных условий и многогрупповые коэффициенты квазидиффузии, которые используются для нахождения нового приближения для групповых плотности и потока излучения. Найденные значения используются для осреднения по энергии и нахождения коэффициентов, необходимых для решения одногрупповой системы уравнений квазидиффузии LO₂ (13).

Всё описанное в абзаце выше представляет собой последовательность действий, которые повторяются во внешнем итерационном процессе. Итерации повторяются до сходимости с заданной точностью $\epsilon$.

Во внутреннем итерационном процессе слагаемые, входящие в правую часть (17), линеаризуют в окрестности температуры, рассчитанной на предыдущей внутренней итерации, считая значения одногрупповой плотности излучения U известными с предыдущей внутренней итерации, при этом производную $\frac{d{Q}_{Pl}}{dt}$ находят аналитически, а производную одногруппового коэффициента поглощения $\frac{\delta {\kappa }_U}{\delta T}$ вычисляют разностно по значениям коэффициента поглощения и температуры, известным из двух предыдущих внешних итераций. В качестве начального приближения для значений на следующем шаге по времени выбираются значения, полученные на текущем шаге.

7. Численное решение задач о взаимодействии излучения с веществом

Предложенная схема ранее была исследована на аналитическом тесте, была показана сходимость схемы с третьим порядком по времени и четвёртым по пространству [23]. Теперь применим разностную схему для решения задач о взаимодействии излучения с веществом. Для удобства будем использовать энергетические единицы для температуры и частоты. Пусть дан одномерный плоскопараллельный слой толщиной L=4 см, на левую границу которого падает излучение с планковским спектром, соответствующим температуре $T=1$ кэВ, коэффициент излучения вещества имеет вид ${\kappa }_{\nu }=\frac{{\kappa }_0}{{\nu }^3}\left(1-e^{-\nu /T}\right)$ см^–1 (значения коэффициента ${\kappa }_0$ см. в Табл. I), уравнение состояния выбрано в форме (18), где $a\text{=8}\text{.1}\cdot {\text{10}}^{-\text{3}}$. Начальная температура слоя $T^{*}={10}^{-4}$.

Пусть $\epsilon ={10}^{-4}$, ${\epsilon }_T={10}^{-6}$. Будем использовать 15-групповое приближение, сетку по частоте (кэВ) зададим набором точек [0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3, 4, 5, 7, 9, 11, 15]. Для дискретизации по угловой переменной выберем сетку из 5 узлов квадратуры Гаусса на отрезках [–1,0] и [0,1]. Будем использовать равномерную сетку с шагом h по пространству, $h=L/N_x$. Для начала применим схему для задачи с непрерывным коэффициентом излучения вещества [24]. В этом случае будет использоваться также равномерная сетка с шагом $\Delta t$ по времени, $\Delta t=t/N_t$.

Первая задача Флека

Полученные профили решений первой задачи Флека представлены на рис. 2. При использовании сетки с $N_x=100$ заметны значительные немонотонности вблизи правой границы отрезка $x\in [0,4]$, а также небольшие немонотонности вблизи левой границы (рис. 2А)). При $N_x=300$ немонотонности остаются только вблизи правой границы отрезка (рис. 2Б)).

$N_t$	$\left|\right|y_{\Delta t}-y_{\Delta t/2}|{|}_C$	Порядок сходимости схемы
15	4.61·10–2	2.03
30	1.13·10–2

$N_x$	$\left|\right|y_h-y_{h/2}|{|}_C$	Порядок сходимости схемы
100	4.21·10–2	3.82
200	2.98·10–3

Порядок сходимости по времени в области монотонного решения равен двум, а порядок сходимости по пространству равен четырём, что демонстрируют Табл. II и Табл. III. Снижение порядка сходимости по времени ожидаемо и связано со структурой погрешности используемого метода Рунге–Кутты [23, 27, 29]. Учёт приграничных областей, в которых значительны немонотонности, снижает порядок сходимости по обеим переменным.

Полученные результаты качественно согласуются с результатами, которые приведены в [24, 25]. Количественное сопоставление проведено с результатами работы [30] и представлено в Табл. IV. Стационарные решения первой задачи Флека, полученные по предложенной в данной работе разностной схеме и по явно-неявной разностной схеме из [30], совпадают с хорошей точностью.

	$x=0.5$	$x=1.0$	$x=1.5$	$x=2.0$	$x=2.5$	$x=3.0$	$x=3.5$
HOLO $ct=6$	0.6925	0.3935	0.1968	0.1127	0.0662	0.0383	0.0234
HOLO $ct=15$	0.8490	0.7683	0.6703	0.5368	0.3784	0.2608	0.1892
HOLO $ct=90$	0.9302	0.9031	0.8764	0.8475	0.8138	0.7705	0.7046
E-I стационар	0.9330	0.9060	0.8792	0.8501	0.8161	0.7726	0.7071

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

Вторая задача Флека

Как уже говорилось выше, из-за высокого порядка аппроксимации схемы может быть значительным влияние немонотонностей, обусловленных как пространственной, так и временной аппроксимацией. Так как во второй задаче Флека коэффициент излучения вещества является разрывной функцией, то влияние немонотонностей на решение становится ещё более выраженным, и интегрирование по времени при помощи метода (19) третьего порядка не позволяет получить решение. Поэтому при решении задач с разрывом коэффициентов для интегрирования по времени используется метод первого порядка – неявный метод Эйлера, который задаётся первой стадией метода (19). Использование метода первого порядка по времени является естественным, так как расщепление по физическим процессам в задачах высокотемпературной радиационной газовой динамики обычно делается с первым порядком аппроксимации по времени. При решении второй задачи Флека, а также предложенных ниже её модификаций используется неравномерная сетка по времени, уменьшение шага $\Delta t$ требуется в те моменты времени, когда температурный фронт входит в область с более высоким коэффициентом поглощения или выходит из неё. На рис. 3Е) можно видеть полученный профиль решения второй задачи Флека, результат хорошо согласуется с результатами, полученными другими методами (рис. 4, Табл. V).

HOLO	Наименьшее значение на отрезке $x\in [0,2]$	Наибольшее значение на отрезке $x\in [2,2.4]$	Значение при $x=2.4$	Значение при $x=3$
$ct=3$	0.031	0.630	–	–
$ct=6$	0.252	0.796	–	–
$ct=9$	0.543	0.835	–	–
$ct=15$	0.837	0.905	–	–
$ct=24$	–	0.959	0.141	–
$ct=36$	–	0.982	0.397	0.031
$ct=400$	–	0.990	0.543	0.507

QD	Наименьшее значение на отрезке $x\in [0,2]$	Наибольшее значение на отрезке $x\in [2,2.4]$	Значение при $x=2.4$	Значение при $x=3$
$ct=3$	0.052	0.637	–	–
$ct=6$	0.317	0.806	–	–
$ct=15$	0.840	0.905	–	–
$ct=24$	–	0.959	0.140	–
$ct=36$	–	0.978	0.309	0.097

Табл. V позволяет провести сравнение результатов решения второй задачи Флека по предложенной бикомпактной схеме с решением квазидиффузионным методом, предложенным в [25]. Значения в таблице для квазидиффузионной схемы получены по рис. 4 из [25] с точностью 0.005. Решение из работы [25] было получено на неравномерной сетке с $h_{\max }=0.4$ и $h_{\min }=0.01$, сетка в области $x\in [2,2.4]$ была более подробной. Поэтому в первую очередь имеет смысл сравнивать значения в области с большим коэффициентом поглощения. Значения решений по предложенной нами схеме и по квазидиффузионной схеме более низкого порядка аппроксимации по пространству в этой области практически не отличаются. В остальных областях результаты качественно не отличаются.

Решения по явно-неявной схеме и MDS_n-схеме взяты из [30,31] (рис. 4). Сравнить результаты расчётов по этим схемам с предложенной нами схемой позволяет рис. 4, из которого видно хорошее совпадение решений в области $x\in [0,2.4]$, за исключением области вблизи $x=2.4$, в которой проявились немонотонности схемы для HOLO-алгоритма. В области $x\in [2.4,4]$ результаты качественно совпадают. Сопоставление с графиками решения по REED-схеме и DDAD-схеме, представленными в [31], также даёт хорошее совпадение. При $x=2.4$ и ct=400 значение по MDS_n-схеме примерно равно 0.516, а значение по схемам REED и DDAD примерно равно 0.535, что ближе к значению 0.543, получаемому в предложенной нами схеме.

Рассмотрим также набор задач, представляющих собой модификацию второй задачи Флека и отличающихся значением коэффициента излучения в центральной области $x\in [2,2.4]$.

Модифицированная вторая задача Флека

Решения модифицированной второй задачи Флека для набора значений ${\rm K}={10}^3,2\cdot {10}^3,4\cdot {10}^3,6\cdot {10}^3,8\cdot {10}^3$ представлены на рис. 3А), 3Б), 3В), 3Г), 3Д). Из рис. 3А) и 3Б) видно, что при значениях коэффициента поглощения ${\rm K}\le 2\cdot {10}^3$ немонотонности проявляются незначительно, а при дальнейшем увеличении величины ${\rm K}$ достигают величины 10–20% от максимального значения температуры и слабо зависят от величины ${\rm K}$.

Рис. 5 демонстрирует, как растёт со временем максимальная температура в оптически толстой области, рис. 6 показывает, как со временем изменяется положение теплового фронта, которое определялось как координата, в которой температура принимает значение 0.5 кэВ. Скорость движения теплового фронта $v_T$ в оптически толстой области зависит от ${\rm K}$, что демонстрирует рис. 7А). Полученная зависимость в двойных логарифмических координатах хорошо приближается линейной функцией (рис. 7Б)). Прямую, построенную по методу наименьших квадратов, также можно видеть на рис. 7Б). После возвращения к исходным координатам полученная аппроксимация зависимости скорости фронта $v_T$ от ${\rm K}$ имеет вид $v_T={10}^{3.150}/{{\rm K}}^{0.633}$.

8. Особенности старта и детали расчёта

Сходимость итерационного процесса зависит от выбора целого набора параметров расчёта. Так как внутренние итерации являются квазиньютоновскими, то важно выбрать не слишком крупный шаг по времени, а также удачное начальное приближение. Модуль производной $\frac{\delta {\kappa }_U}{\delta T}$ необходимо ограничить сверху, так как от итерации к итерации может происходить сильная перестройка спектра.

Несмотря на то что при $t=0$ значения величин, вычисляемых на различных этапах HOLO-алгоритма, полагаются малыми, групповые и осреднённые по энергии значения плотности и потока излучения уже в начальном распределении должны быть согласованными и удовлетворять соотношениям $U={\displaystyle \underset{p=1}{\overset{{\rm P}}{\sum }}{U}^p}$ и $W={\displaystyle \underset{p=1}{\overset{{\rm P}}{\sum }}{W}^p}$. Так как используемые в работе схемы являются кинетически согласованными, то и на следующих шагах по времени одногрупповые плотность и поток излучения равны суммам групповых значений плотности и потока соответственно с хорошей точностью.

В расчётах полагалось, что начальная температура $T(x,0)={10}^{-4}$, начальные значения групповых величин $I^p(x,0,\mu )=2.5\cdot {10}^{-5}$, $U^p(x,0,\mu )=2.5\cdot {10}^{-4}$, $W^p(x,0,\mu )={10}^{-7}$, начальные значения величин, осреднённых по энергии, $U(x,0,\mu )={\rm P}\cdot U^p(x,0,\mu )$, $W(x,0,\mu )={\rm P}\cdot W^p(x,0,\mu )$, производная одногруппового коэффициента поглощения $\left|\frac{\delta {\kappa }_U}{\delta T}\right|\le {10}^{-4}$. На первом шаге по времени и на первой внутренней итерации значения коэффициентов ${\kappa }_U={\kappa }_R={10}^4$, $\zeta =0$. Коэффициенты, входящие в многогрупповую систему уравнений квазидиффузии и уравнение переноса, задавать отдельно не требуется, так как после расчёта температуры во внутреннем итерационном процессе они могут быть рассчитаны по формулам из раздела 2, необходимая для вычисления эффективная температура считалась равной температуре падающего излучения ${\theta }^p=1$ кэВ.

Так как используемая схема является немонотонной, то ещё до достижения сходимости первой внешней итерации могут возникнуть отрицательные значения $I^p$, $U^p$, $U$, T, которые требуется скорректировать. Полагалось, что если $I^p(x,t,\mu )<0$, то $I^p(x,t,\mu )={10}^{-9}$, аналогично если $U^p(x,t)<0$, то $U^p(x,t)={10}^{-8}$, если $U(x,t)<0$, то $U(x,t)={10}^{-8}$. Значения температуры также корректировались: если $T(x,t)<{10}^{-4}$, то $T(x,t)={10}^{-4}$.

Описанный выбор начального приближения и способ корректировки отрицательных значений позволяют достичь сходимости итерационных процессов в первой задаче Флека в достаточно широкой области изменения сеточных параметров. Ещё расширить область значения параметров сетки, при которых имеет место сходимость итерационных процессов, можно, если на первом шаге по времени или первых нескольких шагах по времени использовать для интегрирования по времени вместо метода (19) третьего порядка неявный метод Эйлера первого порядка, который задаётся первой стадией метода (19). Если и это не приводит к сходимости итерационного процесса, то можно сделать расчёт температуры на первом слое по времени в диффузионном приближении. В таком случае расчёт не включает в себя решение уравнения переноса, групповые значения коэффициента квазидиффузии полагаются равными 1/3, коэффициенты граничных условий $c_0^p=-0.5$, $c_L^p=0.5$. Полученные значения температуры используются как нулевое приближение во внутреннем итерационном процессе в расчёте по полному алгоритму. Использование описанных подходов позволяет получить решение первой задачи Флека при числе точек по пространству от 30 до 3200. Расчёт решения при $N_x>3200$ требует большого времени, поэтому более подробные сетки не рассматривались. Шаг по времени обычно выбирался достаточно крупным, порядка 10^–3, и не изменялся в ходе расчёта. На пространственной сетке с крупным шагом сходимость возможна при довольно крупных шагах и по времени, на подробных сетках по пространству расчёт возможен и с более мелким шагом по времени.

При решении второй задачи Флека и её модификаций на всех шагах по времени используется метод Эйлера первого порядка, как и было сказано в предыдущем разделе. При значениях коэффициента поглощения в центральной области ${\rm K}\le {10}^3$ никаких дополнений к описанным выше подходам не требуется, но приходится измельчать шаг по времени при входе тепловой волны в оптически толстую область и выходе из неё. При ${\rm K}>{10}^3$ возникает необходимость начинать расчёт с ${{\rm K}}_0={10}^3$, $\Delta t=5\cdot {10}^{-4}$, а на следующих шагах увеличивать значение ${\rm K}$ и уменьшать $\Delta t$. В данной работе коэффициент ${\rm K}$ увеличивался обычно скачкообразно в два-три этапа. Увеличение до заданного в задаче значения ${\rm K}$ необходимо осуществить к моменту времени, когда тепловая волна достигнет оптически толстого слоя. Область значений сеточных параметров, при которых достигается сходимость итерационных процессов во второй задаче Флека и её модификациях, значительно более узкая. Сетка по пространству также выбиралась равномерной, использование неравномерной сетки, более подробной в центральной области, приводит к более сильному развитию немонотонностей, а иногда и к отсутствию сходимости итерационного процесса.

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

Для решения модельных задач переноса излучения в среде были применены бикомпактные схемы высокого порядка аппроксимации для HOLO-алгоритма решения уравнения переноса.

Схемы были применены для решения первой задачи Флека. Был показан четвёртый порядок сходимости схемы по пространству в области монотонного решения, что соответствует четвёртому порядку аппроксимации схемы. Показано, что порядок сходимости по времени в области монотонного решения равен двум, что меньше третьего порядка аппроксимации схемы по времени. Снижение порядка сходимости связано со структурой погрешности используемого метода Рунге–Кутты и использованием классических граничных условий для одногрупповой и многогрупповой систем уравнений квазидиффузии. Была проведена валидация полученного решения первой задачи Флека по решениям, полученным другими методами, было получено хорошее совпадение результатов по предложенной схеме и по явно-неявной схеме Н.Я. Моисеева.

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

Было предложено использовать для тестирования схем модификацию второй задачи Флека, которая отличается от второй задачи Флека величиной коэффициента поглощения ${\rm K}$ в центральной области. Было показано, что при ${\rm K}\le 2\cdot {10}^3$ решение слабо подвержено немонотонностям, при увеличении ${\rm K}$ немонотонности становятся значительными. Была исследована зависимость скорости движения фронта $v_T$ в оптически толстой области от величины ${\rm K}$ и было получено, что эта зависимость хорошо описывается формулой $v_T={10}^{3.150}/{{\rm K}}^{0.633}$.

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