Журнал вычислительной математики и математической физики, 2021, T. 61, № 6, стр. 990-1018

ВВЕДЕНИЕ

Алгоритм ускорения итераций по интегралу рассеяния является существенным элементом численной методики решения уравнения переноса, основанной на использовании метода дискретных ординат. Расчеты полей излучения в активной зоне и радиационной защите ядерно-технических установок в неодномерных геометриях требуют значительных затрат процессорного времени. В настоящее время такие расчеты проводятся, в основном, в 3D геометриях. Однако расчеты в 2D геометриях также являются востребованными, позволяя провести оценку радиационных полей за относительно небольшие времена с использованием мультигрупповых константных систем. Поэтому, с практической точки зрения, разработка эффективного алгоритма ускорения для 2D геометрий с использованием нодальных схем высокого порядка точности является актуальной задачей, позволяя улучшить аппроксимацию уравнения переноса и сократить время расчета типичного варианта в 3–10 и более раз. В данной работе нами будет построена $K{{P}_{1}}$ схема ускорения внутренних итераций для уравнения переноса в 2D $r,z$, $x,z$ и $r,\vartheta $ геометриях, согласованная с нодальными Linear Discontinues (LD) и Linear Best (LB) схемами 3 -го и 4-го порядков точности.

Предлагаемый алгоритм представляет собой реализацию $K{{P}_{1}}$ схемы (см. [1]) (известной также как DSA (см. [2], [3]) или ${{P}_{1}}SA$ схема (см. [4]–[6])) ускорения внутренних итераций для практически важного случая LD и LB схем (см. [7]–[13]) в 2D $r,z$, $x,z$ и $r,\vartheta $ геометрий, позволяющих обеспечить быструю сходимость разностного решения при сгущении пространственной сетки задачи. Согласованность построенной $K{{P}_{1}}$ схемы ускорения с разностной аппроксимацией уравнения переноса, как и в случае взвешенной алмазной (WDD) схемы из [2], обеспечивает ее устойчивость и слабую зависимость эффекта ускорения от выбора пространственной сетки. Данный алгоритм реализован в 2D ${{S}_{n}}$ программе КАСКАД-С-3.5 из пакета программ CNCSN (см. [14]). Он представляет собой развитие $K{{P}_{1}}$ схемы ускорения, согласованной со взвешенной алмазной (WDD) схемой, предложенной ранее в [6], [15]–[17] и реализованной в 1D, 2D и 3D ${{S}_{n}}$ программах из пакета CNCSN, а также $K{{P}_{1}}$ схемы ускорения, согласованной со взвешенной WLD-WLB схемой в 1D геометриях (см. [18], [19]).

Обзор предшествующих работ по $K{{P}_{1}}$ схеме ускорения и другим вариантам схем ускорения внутренних итераций имеется в [1], [3], [6], [16]. Построение $K{{P}_{1}}$ схемы ускорения мы проведем на примере уравнения переноса в $r,z$ геометрии.

Последовательность изложения в данной работе следующая. В разд. 1 и 2 рассмотрено построение WLD-WLB схемы и согласованной $K{{P}_{1}}$ схемы ускорения для уравнения переноса в $r,z$ геометрии. В разд. 3 сформулирован метод расщепления (МР) для решения ${{P}_{1}}$ системы для поправок, включая метод прогонки для решения вспомогательных двухточечных уравнений в $r,z$ геометрии. В разд. 4 рассмотрен алгоритм оценки границ спектра радиальной ($\hat {R}$) и аксиальной ($\hat {Z}$) компонент ${{P}_{1}}$ оператора, которые используются для выбора шагов МР для уравнения переноса в $r,z$ геометрии. В разд. 5 рассмотрен вопрос о выборе итерационных параметров циклического МР. В разд. 6 указаны изменения, которые нужно внести в алгоритм ускорения для случая $x,z$ геометрии. Расчетные формулы для $r,\vartheta $ геометрии могут быть получены аналогично случаю $r,z$ геометрии. В разд. 7 представлены численные примеры использования $K{{P}_{1}}$ схемы ускорения при решении задач переноса излучения в двумерных геометриях, а также приведены численные результаты, позволяющие оценить скорость сходимости DD, LD и LB схем в $x,z$ и $r,z$ геометриях при сгущении пространственной сетки задачи.

1. WLB-WLD СХЕМА 2–4 -ГО ПОРЯДКОВ ТОЧНОСТИ ДЛЯ УРАВНЕНИЯ ПЕРЕНОСА В ДВУМЕРНОЙ $r,z$ ГЕОМЕТРИИ

В $r,z$ геометрии уравнение переноса имеет вид

где $\xi $, $\eta $ и $\mu $ – направляющие косинусы единичного вектора ${\mathbf{\Omega }}$ направления скорости частицы: который с учетом имеющейся в данной геометрии симметрии: $\psi (r,z,\mu ,\varphi ) = \psi (r,z,\mu ,2\pi - \varphi )$, изменяется в пределах четырех октантов:$ - 1 \leqslant \xi $, $\mu \leqslant 1$, $0 \leqslant \eta \leqslant 1$, $0 \leqslant \varphi \leqslant \pi $; пространственные переменные изменяются в пределах: $0 \leqslant {{r}_{{\operatorname{int} }}} \leqslant r \leqslant {{r}_{{{\text{out}}}}}$, ${{z}_{{{\text{bot}}}}} \leqslant z \leqslant {{z}_{{{\text{top}}}}}$. Правую часть уравнения (1.1), учитывая имеющуюся симметрию, можно представить в виде где $Y_{l}^{m}(\mu ,\varphi )$ – сферические гармоники: $\Phi _{l}^{m}(r,z)$ – угловые моменты потока:

Для аппроксимации интеграла рассеяния вводится квадратура, которая разбивает рассматриваемые четыре октанта $ - 1 \leqslant \xi $, $\mu \leqslant 1$, $0 \leqslant \theta \leqslant \pi $, $0 \leqslant \varphi \leqslant \pi $ на $L$ полос по переменной $\mu $ веса ${{w}_{l}}$, а каждая полоса $l$, l = 1, 2, …, L, дополнительно еще разбивается на ${{M}_{l}}$ секторов веса ${{w}_{{l,m}}}$, где индекс $m$ на $l$-м слое пробегает значения m = 1, 2, …, M_l. Нумерация узлов по $\varphi $ на $l$-м слое соответствует возрастанию ${{\xi }_{{l,m}}}$ (убыванию ${{\varphi }_{{l,m}}}$):

В (1.4) $M_{l}^{ - }$и $M_{l}^{ + }$ задают, соответственно, число узлов с ${{\varphi }_{{l,m}}} \in [{\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2},\pi ]$ и ${{\varphi }_{{l,m}}} \in [0,{\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2}]$ на слое ${{\mu }_{l}}$; угловые направления на слое нумеруются в направлении убывания $\varphi $, $M$ – полное число узлов квадратуры. В $r,z$ геометрии дополнительно вводятся узлы ${{\varphi }_{{l,1}}} = \pi $ с нулевыми весами ${{w}_{{l,1}}} = 0$. В $r,z$ геометрии эти направления используются как граничные для начала расчета слоя с $\mu = {{\mu }_{l}}$, m =1, 2, …, M_l.

Для граничного направления $\varphi = \pi $, используемого для организации расчета граничных ячеек в $r,z$ геометрии, дивергентная форма уравнения (1.1) имеет вид

Уравнение баланса нулевого порядка получается путем интегрирования уравнения (1.1) по разностной ячейке $({{r}_{{i - 1/2}}},{{r}_{{i + 1/2}}}) \times ({{z}_{{k - 1/2}}},{{z}_{{k + 1/2}}}) \times ({{\varphi }_{{l,m - 1/2}}},{{\varphi }_{{l,m + 1/2}}}) \times ({{\mu }_{{l - 1/2}}},{{\mu }_{{l + 1/2}}})$, i = 1, 2, …, I, K = 1, 2, …, K:

где Здесь ${{\psi }_{{R(L)}}}$ и ${{\psi }_{{T(B)}}}$ – нулевые пространственные моменты потока на пространственных гранях ячейки.

Для построения нодальной WLB-WLD схемы 2–4 -го порядков точности (см. [7], [10], [12]), кроме уравнения баланса нулевого порядка (1.6), нам потребуются также уравнения баланса первого порядка по ячейке, получаемые интегрированием уравнения (1.1) по разностной ячейке с весами $2(r - {{r}_{i}} - \delta _{i}^{c}){\text{/}}\Delta r$ и $2(z - {{z}_{k}}){\text{/}}\Delta z$, где $\delta _{i}^{c}$ – малое отклонение от центра ячейки, выбираемое из условия $\int_{{{r}_{{i - 1/2}}}}^{{{r}_{{i + 1/2}}}} {(r - {{r}_{i}} - \delta _{i}^{c})rd} r = 0$ (см. [9]–[11]):

где $\psi $ – среднее значение (нулевой пространственный момент) потока в ячейке, ${{\psi }^{r}}$ и ${{\psi }^{z}}$ – первые пространственные моменты потока по переменным $r$ и $z$: (аналогично определяются и пространственные моменты источника $S$, ${{S}^{r}}$ и ${{S}^{z}}$), ${{\psi }_{{R\{ L\} }}}$ и ${{\psi }_{{T(B)}}}$ – нулевые, а $\psi _{{R(L)}}^{z}$ и $\psi _{{T(B)}}^{r}$ – первые пространственные моменты потока на пространственных гранях ячейки: где ${{\psi }_{{m \pm 1/2}}}$, $\psi _{{m \pm 1/2}}^{r}$ и $\psi _{{m \pm 1/2}}^{z}$ – нулевой и первые пространственные моменты решения на входящей ($m - 1{\text{/}}2$) и выходящей ($m + 1{\text{/}}2$) угловой грани ячейки:

Выбор нормировочного коэффициента $V_{{i,k}}^{r}$ в уравнении (1.12) в определении первого радиального пространственного момента ${{\psi }^{r}}$ (1.10) гарантирует выполнение следующего условия: дополнительное уравнение LB схемы из [9] ${{\psi }_{{i + 1/2}}} = 2{{\psi }^{r}} + {{\psi }_{{i - 1/2}}}$ по $r$ удовлетворяется точно на линейных по $r$ решениях вида $\psi (r) = a + br$ из [11], что улучшает аппроксимацию LB схемы.

Для получения разностной схемы к трем точным балансным уравнениям (1.6), (1.8) и (1.9) следует добавить семь дополнительных уравнений. WLB-WLD схеме соответствует следующий выбор этих уравнений (см. [7], [10], [12], [13]):

Веса ${{P}_{t}}$ и ${{Q}_{t}}$, $t = r,z$, в уравнениях (1.13) изменяются в пределах (см. [12]) где

Случай ${{P}_{t}} = {{Q}_{t}} = 1$, ${{T}_{t}} = 0$ соответствует LB схеме 4 -го порядка точности, случай ${{P}_{t}} = 0$, ${{Q}_{t}} = 1$, ${{T}_{t}} = 0$ соответствует LD схеме 3 -го порядка точности.

Как и в случае 1D криволинейных геометрий дополнительное уравнение WDD схемы, связывающее входящие и выходящие радиальные пространственные моменты нулевого порядка ${{\psi }_{{m \pm 1/2}}}$ на гранях ${{\xi }_{{l,m \pm 1/2}}}$ разностной ячейки, включает член ~${{\psi }^{r}}$, позволяющий учесть различие в пространственном весовом множителе в определении $\psi $ и ${{\psi }_{{m \pm 1/2}}}$. Отметим, что дополнительные уравнения (1.14) выполняются точно на постоянных по углу $\varphi $ и линейных по пространственным переменным $r$ и $z$ решениях:

В дальнейшем ограничим выбор весов $P_{\xi }^{r}$ и $P_{\xi }^{z}$ в экстраполяционных соотношениях для $\psi _{{m + 1/2}}^{r}$ и $\psi _{{m + 1/2}}^{z}$ (1.14) условием

Система дополнительных уравнений LD схемы для случая $x,y$ геометрии была предложена в [8]. Система дополнительных уравнений WLB-WLD схемы (1.13)–(1.15) представляет собой обобщение указанной системы на случай семейства взвешенных нодальных схем в $r$, $z$ геометрии. Отметим, что альтернативный вариант построения LD схемы в $r,z$ геометрии рассмотрен в [20].

Подстановка дополнительных уравнений (1.13), (1.14) в балансные соотношения (1.6), (1.8) и (1.9) приводит к следующей системе уравнений относительно пространственных моментов решения в ячейке $\psi $, ${{\psi }^{r}}$ и ${{\psi }^{z}}$:

где элементы матриц $A$ и $B$ имеют вид Решение системы (1.19) находится по формулам Крамера.

Система балансных уравнений для граничных ячеек с $m = 1{\text{/}}2$ (${{\xi }_{{l,1/2}}} = - \sqrt {1 - \mu _{l}^{2}} $), используемых для нахождения граничных значений потока ${{\psi }_{{i,k,l,1/2}}}$, $\psi _{{i,k,l,1/2}}^{r}$ и $\psi _{{i,k,l,1/2}}^{z}$ по угловой переменной $\xi $, состоит из уравнения баланса нулевого порядка (1.5), в котором полагается

учитывающее линейную аппроксимацию решения в ячейке (1.17), и уравнений баланса первого порядка, которые получаются путем интегрирования уравнения (1.5) по граничной ячейке с весами $2(r - {{r}_{i}} - \delta _{i}^{c}){\text{/}}\Delta r$ и $2(z - {{z}_{k}}){\text{/}}\Delta z$: В случае WLB-WLD схемы эти уравнения дополняются четырьмя дополнительными уравнениями (1.13). Подставляя их в балансные соотношения (1.6), (1.23) и (1.24), получаем систему вида (1.19) со следующими коэффициентами:

Следует отметить, что расчет ячейки с использованием WLB-WLD схемы происходит аналогично случаю WDD схемы. Сначала из решения системы вида (1.19) находятся нулевой и первые пространственные моменты решения в ячейке, а затем по явным формулам (1.13), (1.14) находятся граничные потоки и соответствующие пространственные моменты граничных потоков на выходящих гранях ячейки.

Как показывает геометрическая интерпретация WLB-WLD схемы из [12], путем соответствующего выбора весов в дополнительных уравнениях (1.13), (1.14) можно обеспечить положительность схемы.

2. ПОСТРОЕНИЕ $K{{P}_{1}}$ СХЕМЫ УСКОРЕНИЯ ВНУТРЕННИХ ИТЕРАЦИЙ В $r$, $z$ ГЕОМЕТРИИ

Перепишем уравнение баланса и систему дополнительных уравнений WLD-WLB схемы в $r$, $z$ геометрии на $(n + 1)$-й внутренней итерации в более удобном для последующего изложения виде:

Дополним их системой дополнительных уравнений Веса ${{T}_{t}}$, $t = r,z$ отличны от нуля только при использовании очень редких пространственных сеток. Их введение усложняет систему дополнительных уравнений для поправок. Однако, если пространственная сетка задачи не слишком грубая, такая коррекция, как показал численный эксперимент, вообще говоря, не требуется. В данной работе вышеуказанные весовые коэффициенты полагаются равными нулю: ${{T}_{t}} = 0$, $t = r,z$.

Дополнительные уравнения для пространственных моментов на гранях ячейки в случае ${{T}_{r}} = {{T}_{z}} = 0$ можно записать в виде

Дополнительное WDD уравнение по угловой переменной можно записать в виде К системе уравнений (2.1)–(2.8) следует также добавить граничные условия при $r = {{r}_{{1/2}}} = {{r}_{{\operatorname{int} }}} \geqslant 0$, $r = {{r}_{{I + 1/2}}} = {{r}_{{{\text{ext}}}}}$, $z = {{z}_{{1/2}}} = {{z}_{{{\text{bot}}}}}$, $z = {{z}_{{K + 1/2}}} = {{z}_{{{\text{top}}}}}$: Для расчета весового коэффициента ${{P}_{\xi }}$ в дополнительном WDD уравнении по угловой переменной используется адаптивная WDD (AWDD) схема (см. [10], [21]–[23].

Для ускорения сходимости внутренних итераций в $K{{P}_{1}}$ схеме используются линейные поправки к нулевому и первым угловым моментам решения следующего вида:

где $\alpha $ и $\beta $ – индексы пространственных моментов. Уравнение (2.11) содержит девять поправок, уравнения (2.12) – восемь поправок.

Для получения системы уравнений для определения ускоряющих поправок ${{f}^{{\alpha ,0}}}$, ${{f}^{{\alpha ,r}}}$ и ${{f}^{{\alpha ,z}}}$ воспользуемся следующей процедурой, близкой к “4-step” процедуре Ларсена (см. [3]). Прежде всего подействуем на уравнение баланса (1.11.3) операторами проектирования ${{\hat {L}}_{0}}$, ${{\hat {L}}_{r}}$, и ${{\hat {L}}_{z}}$:

где интегрирование проводится по четырем октантам. В предположении, что используемая квадратура достаточно точно интегрирует соответствующие сферические гармоники, получим где

Величины $\pi _{{l,m \pm 1/2}}^{r}$, $\bar {\pi }_{{l,m \pm 1/2}}^{r}$ в (2.18) могут быть вычислены с помощью следующих рекуррентных соотношений:

В систему из девяти уравнений (2.14)–(2.16) входит 17 неизвестных величин: $f_{{i,k}}^{{\alpha ,0}}$, $f_{{i,k}}^{{\alpha ,r}}$, $f_{{i,k}}^{{\alpha ,z}}$, $\alpha = 0,r,z$; $f_{{i \pm 1/2,k}}^{{\beta ,0}}$, $f_{{i \pm 1/2,k}}^{{\beta ,r}}$, $\beta = 0,z$; $f_{{i,k \pm 1/2}}^{{\beta ,0}}$, $f_{{i,k \pm 1/2}}^{{\beta ,z}}$, $\beta = 0,r$. Для замыкания системы (2.14) к ней необходимо добавить восемь дополнительных уравнений, а также граничные условия. Искомые дополнительные уравнения получаются применением операторов ${{\hat {L}}_{0}}$ и ${{\hat {L}}_{r}}$ к уравнениям (2.4), (2.6) и операторов ${{\hat {L}}_{0}}$ и ${{\hat {L}}_{z}}$ к уравнениям (2.5), (2.7): где Для получения граничных условий для ${{P}_{1}}$ системы для ускоряющих поправок воспользуемся следующей процедурой. Домножим уравнения (2.9), определяющие закон отражения при $r = {{r}_{{\operatorname{int} }}}$ и $r = {{r}_{{{\text{ext}}}}}$, на ${{w}_{{l,m}}}{{\xi }_{{l,m}}}$ и просуммируем их по направлениям ${{\xi }_{{l,m}}} > 0$ и ${{\xi }_{{l,m}}} < 0$ соответственно. С учетом вида поправок (2.12) получим где Для определения граничных условий для поправок на границах $z = {{z}_{{{\text{bot}}}}}$ и $z = {{z}_{{{\text{top}}}}}$ домножим уравнения (2.10), определяющие закон отражения на этих границах, на ${{w}_{{l,m}}}{{\mu }_{l}}$ и просуммируем их по направлениям ${{\mu }_{l}} > 0$ и ${{\mu }_{l}} < 0$ соответственно. С учетом вида поправок (2.12) получим где Отметим, что в уравнениях (2.25) и (2.27) величины $\lambda _{{j,k}}^{{\beta ,r}}$ и $\lambda _{{i,j}}^{{\beta ,z}}$ стремятся к нулю по мере сходимости итераций. По этой причине возможно использование также и упрощенных граничных условий на границах $r = {{r}_{{{\text{ext}}}}}$ и $z = {{z}_{{{\text{top}}}}}$, отличающихся от вышеприведенных тем, что правые части в них полагаются равными нулю: Использование упрощенных граничных условий, в свою очередь, позволяет упростить постановку граничных условий в методе расщепления, используемом для нахождения решения ${{P}_{1}}$ системы для ускоряющих поправок.

Учитывая, что в 2D и 3D геометриях в задачах с существенной ролью анизотропии рассеяния, когда выполняется неравенство

сходимость $K{{P}_{1}}$ схемы нарушается (см. [1], [24]), для обеспечения сходимости $K{{P}_{1}}$ схемы в задачах переноса нейтрального излучения нами используется следующий простой прием (см. [6], [17]): на четных внутренних итерациях ускоряющая коррекция в соответствии с (2.12) проводится для четырех угловых моментов решения, включая нулевой, на нечетных – только для скалярного потока. Ясно, что такое решение не является оптимальным, но оно обеспечивает сходимость $K{{P}_{1}}$ схемы для широкого класса практических задач.

3. АЛГОРИТМ РЕШЕНИЯ ${{P}_{1}}$ СИСТЕМЫ ДЛЯ УСКОРЯЮЩИХ ПОПРАВОК В $r,z$ ГЕОМЕТРИИ

${{P}_{1}}$ система для $K{{P}_{1}}$ схемы в $r,z$ геометрии состоит из девяти балансных уравнений (2.14)–(2.16) и восьми дополнительных уравнений: (2.20)–(2.24) и граничных уравнений (2.25), (2.27). Эта система может быть решена итерационно с использованием, например, метода расщепления (МР) (см. [1], [25]).

В соответствии с общим подходом МР оператор ${{P}_{1}}$ системы $\hat {A}$ для ускоряющих поправок

представляется в виде суммы: где $\hat {R}$ и $\hat {Z}$ – компоненты оператора $\hat {A}$, определяющие изменение решения ${{P}_{1}}$ системы по переменным $r$ и $z$ соответственно. Более подробно о структуре операторов $\hat {R}$ и $\hat {Z}$ будет сказано ниже. Предположим также, что указанные операторы являются положительно-определенными. В МР сложился следующий конструктивный способ построения итерационной схемы. Для приближенного решения уравнения (3.1) рассмотрим неявную двухслойную итерационную схему вида (см. [1], [25]) где ${{\tau }_{{s + 1}}}$, $\tau _{\alpha }^{{(s + 1)}}$, $\alpha = r,z$ – итерационные параметры, $B$ – регуляризирующий оператор, $s$ – номер итерации. Если операторы $\hat {R}$ и $\hat {Z}$ коммутируют, имеют положительные собственные значения и итерационные параметры не зависят от номера итерации $s$, то достаточным условием сходимости итерационного процесса (3.3) является выполнение неравенства (см. [1], [25]) Стандартный выбор параметров ${{\tau }_{\alpha }} = \tau {\text{/}}2$ обеспечивает второй порядок аппроксимации по времени (см. [25]), если рассматривать итерационный процесс (3.3) как решение на установление нестационарной задачи Наряду с параметрами ${{\tau }_{\alpha }}$ введем также величины $\omega _{\alpha }^{{(s)}} \equiv {1 \mathord{\left/ {\vphantom {1 {\tau _{\alpha }^{{(s)}}}}} \right. \kern-0em} {\tau _{\alpha }^{{(s)}}}}$. Используя эти величины, перепишем уравнение (3.5) в виде Отметим, что при $\tau _{\alpha }^{{(s)}} = {{\tau }^{{(s)}}}{\text{/}}2$ величина ${{\tilde {\tau }}_{s}} = 2{{\omega }^{{(s)}}}$. Схема реализации алгоритма (3.6) имеет следующий вид (см. [1], [25], [27]):

Подчеркнем, что, исключая из оператора $\hat {A}$ с помощью дополнительных уравнений (2.20)–(2.21), (2.23)–(2.24) и граничных условий (2.25) и (2.27) поправки, относящиеся к граням ячеек, можно представить его в виде, в котором присутствуют только значения поправок, отнесенные к центру ячейки. Соответственно, в (3.7) можно считать, что единичный оператор $\hat {E}$ действует только на поправки, отнесенные к центру ячейки.

Учитывая вид поправок (2.12), используемых в $K{{P}_{1}}$ алгоритме, и применяя двухшаговую процедуру метода расщепления (3.7), получаем следующую систему для поправок на первом шаге $\tau _{r}^{{(s + 1)}}$:

Отметим, что в системе уравнений (3.8)–(3.10) и ниже, как и в реализации МР для WDD схемы из [17], используется симметричная декомпозиция членов ∼${{\sigma }^{{00}}}$ и ${{\sigma }^{{11}}}$ по операторам $\hat {R}$ и $\hat {Z}$ (шагам $\tau _{r}^{{(s + 1)}}$ и $\tau _{z}^{{(s + 1)}}$). Это позволяет сохранить вид операторов $\hat {R}$ и $\hat {Z}$ подобным виду P1 оператора для ускоряющих поправок в одномерной геометрии (см. [15], [19]).

Далее мы будем говорить, что величина относится к слою $s$, $s + 1{\text{/}}2$ и т.д., если она имеет соответствующий индекс. Отметим, что в рассматриваемом нами случае 2D геометрии величинa $\zeta $ для дробного шага $s + 1{\text{/}}2$ не аппроксимируeт исходную нестационарную задачу (3.5) (см. [26]). Для получения граничных условий по радиальной переменной для поправок $\zeta $ отметим, что из граничных условий (2.25) для величин $f$ в предположении (2.29), а также уравнения (3.7) следует, что

Так как на шаге $s = 0$ поправки $f$ равны нулю, а на последующих шагах граничные соотношения для поправок $f$ выполняются точно, можно считать, что

Система (3.8)–(3.10) для девяти уравнений, дополненная двумя дополнительными уравнениями для величин ${{\zeta }^{{r,0,s + 1/2}}}$и ${{\zeta }^{{r,r,s + 1/2}}}$, аналогичным (2.20), и двумя дополнительными уравнениями (2.23) для величин ${{\zeta }^{{z,0,s + 1/2}}}$ и ${{\zeta }^{{z,r,s + 1/2}}}$, а также граничными условиями (3.12), может быть решена прямым методом в предположении, что величины с индексом $s$ известны с предыдущего шага. Данная система состоит из двух подсистем, которые могут быть решены независимо методом прогонки по переменной $r$, а также трех уравнений – последних из подсистем уравнений (3.8)–(3.10), которые могут быть явно разрешены относительно переменных ${{\zeta }^{{0,z,s + 1/2}}}$, ${{\zeta }^{{r,z,s + 1/2}}}$ и ${{\zeta }^{{z,z,s + 1/2}}}$.

Первая подсистема относительно величин ${{\zeta }^{{0,0,s + 1/2}}}$, ${{\zeta }^{{0,r,s + 1/2}}}$, ${{\zeta }^{{r,0,s + 1/2}}}$, ${{\zeta }^{{r,r,s + 1/2}}}$, $\zeta _{{i \pm 1/2}}^{{0,r,s + 1/2}}$ и $\zeta _{{i \pm 1/2}}^{{0,0,s + 1/2}}$ состоит из двух первых уравнений в подсистемах уравнений (3.8) и (3.9) (всего четыре уравнения), дополнительных уравнений (2.20) для величин ${{\zeta }^{{r,0,s + 1/2}}}$и ${{\zeta }^{{r,r,s + 1/2}}}$ и граничных условий (3.12). Для ее решения мы должны исключить из нее величины ${{\zeta }^{{0,0,s + 1/2}}}$, ${{\zeta }^{{0,r,s + 1/2}}}$, ${{\zeta }^{{r,0,s + 1/2}}}$ и ${{\zeta }^{{r,r,s + 1/2}}}$, чтобы получить двухточечную систему для величин $\zeta _{{i \pm 1/2}}^{{0,r,s + 1/2}}$ и $\zeta _{{i \pm 1/2}}^{{0,0,s + 1/2}}$ (см. [19]). Исключая из первых двух уравнений подсистемы (3.9) посредством дополнительных уравнений (2.20) величины ${{\zeta }^{{r,0,s + 1/2}}}$и ${{\zeta }^{{r,r,s + 1/2}}}$, получаем следующую систему линейных уравнений относительно ${{\zeta }^{{0,0,s + 1/2}}}$ и ${{\zeta }^{{0,r,s + 1/2}}}$:

где Из системы (3.13) может быть найдена явная зависимость величин ${{\zeta }^{{0,0}}}$ и ${{\zeta }^{{0,r}}}$ от $\zeta _{{i \pm 1/2}}^{{0,r}}$ и $\zeta _{{i \pm 1/2}}^{{0,0}}$: где коэффициенты $x_{1}^{j}$ и $x_{2}^{j}$, $j = \overline {1,5} $, находятся путем решения пяти систем из двух линейных уравнений для каждой пространственной ячейки расчетной области: Массивы коэффициентов $x_{1}^{j}$ и $x_{2}^{j}$, $j = \overline {1,5} $, хранятся. Подставляя уравнения (3.14) (отнесенных к слою $s + 1{\text{/}}2$) в первые два из уравнений подсистемы (3.8), получаем следующую двухточечную систему из двух уравнений для определения величин $\zeta _{{i \pm 1/2,k}}^{{0,s + 1/2}}$ и $\zeta _{{i \pm 1/2,k}}^{{r,s + 1/2}}$: где Для каждого значения $k$ система (3.16), (3.17) решается методом прогонки (см. [1], [15], [27]): где коэффициенты ${{{\mathbf{\xi }}}_{i}} = \operatorname{col} \{ \xi _{i}^{0},\xi _{i}^{r}\} $ и ${{\eta }_{i}} = \operatorname{col} \{ \eta _{i}^{0},\eta _{i}^{r}\} $ вычисляются посредством рекуррентных соотношений К рекуррентным соотношениям (3.18), (3.19) следует также добавить граничные условия на внешней $r = {{r}_{{{\text{ext}}}}}$ и внутренней $r = {{r}_{{\operatorname{int} }}}$ границах области: Достаточным условием устойчивости этого алгоритма прогонки является выполнение неравенства $\left| {\xi _{i}^{0}} \right| < 1$. После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{i \pm 1/2}}^{{0,r,s + 1/2}}$ и $\zeta _{{i \pm 1/2}}^{{0,0,s + 1/2}}$, поправки ${{\zeta }^{{0,0,s + 1/2}}}$, ${{\zeta }^{{0,r,s + 1/2}}}$, ${{\zeta }^{{r,0,s + 1/2}}}$ и ${{\zeta }^{{r,r,s + 1/2}}}$ для средних значений потока вычисляются по явным формулам (3.14) и (2.20). Величины ${{\zeta }^{{0,z,s + 1/2}}}$ и ${{\zeta }^{{r,z,s + 1/2}}}$ находятся по явным формулам из последних двух уравнений подсистем (3.8) и (3.9):

Вторая подсистема для величин ${{\zeta }^{{z,0,s + 1/2}}}$, ${{\zeta }^{{z,r,s + 1/2}}}$, $\zeta _{{i \pm 1/2}}^{{z,0,s + 1/2}}$ и $\zeta _{{i \pm 1/2}}^{{z,r,s + 1/2}}$ состоит из двух первых уравнений в подсистеме уравнений (3.10), дополнительных уравнений (2.23) и граничных условий (3.12). Для ее решения исключим из нее посредством дополнительных уравнений (2.23) величины ${{\zeta }^{{z,0,s + 1/2}}}$ и ${{\zeta }^{{z,r,s + 1/2}}}$, чтобы получить двухточечную систему для величин $\zeta _{{i \pm 1/2}}^{{z,0,s + 1/2}}$ и $\zeta _{{i \pm 1/2}}^{{z,r,s + 1/2}}$:

где Для каждого значения $k$ система (3.22), (3.23) решается методом прогонки (3.18) с граничными условиями (3.12). После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{i \pm 1/2}}^{{z,0,s + 1/2}}$ и $\zeta _{{i \pm 1/2}}^{{z,r,s + 1/2}}$, поправки ${{\zeta }^{{z,0,s + 1/2}}}$ и ${{\zeta }^{{z,r,s + 1/2}}}$ вычисляются по явным формулам (2.23). Величина ${{\zeta }^{{z,z,s + 1/2}}}$ находится по явной формуле из последнего уравнения подсистемы (3.10):

На втором шаге $\tau _{z}^{{(s + 1)}}$ алгоритма МР решается следующая система:

Граничные условия для системы (3.25)–(3.27) с учетом предположения (2.29) и последнего из уравнений алгоритма (3.7) имеют вид

Система (3.25)–(3.27) для девяти уравнений, дополненная двумя дополнительными уравнениями для величин ${{\zeta }^{{z,0,s + 1}}}$ и ${{\zeta }^{{z,z,s + 1}}}$, аналогичным (2.21), и двумя дополнительными уравнениями (2.24) для величин ${{\zeta }^{{r,0,s + 1}}}$, ${{\zeta }^{{r,z,s + 1}}}$, ${{\zeta }^{{\vartheta ,0,s + 1}}}$ и ${{\zeta }^{{\vartheta ,z,s + 1}}}$, а также граничными условиями (3.28), может быть решена прямым методом в предположении, что величины с индексом $s + 1{\text{/}}2$ известны с предыдущего шага. Данная система состоит из двух подсистем, которые могут быть решены независимо методом прогонки по переменной $z$, а также трех уравнений – вторых из подсистем уравнений (3.25)–(3.27), которые могут быть явно разрешены относительно переменных ${{\zeta }^{{0,r,s + 1}}}$, ${{\zeta }^{{r,r,s + 1}}}$ и ${{\zeta }^{{z,r,s + 1}}}$.

Первая подсистема относительно величин ${{\zeta }^{{0,0,s + 1}}}$, ${{\zeta }^{{0,z,s + 1}}}$, ${{\zeta }^{{z,0,s + 1}}}$, ${{\zeta }^{{z,z,s + 1}}}$, $\zeta _{{k \pm 1/2}}^{{0,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{0,z,s + 1}}$ состоит из первого и третьего уравнений в подсистемах уравнений (3.25) и (3.27) (всего четыре уравнения), дополнительных уравнений (2.21) для величин ${{\zeta }^{{z,0,s + 1}}}$ и ${{\zeta }^{{z,z,s + 1}}}$ и граничных условий (3.28). Для ее решения мы должны исключить из нее величины ${{\zeta }^{{0,0,s + 1}}}$, ${{\zeta }^{{0,z,s + 1}}}$, ${{\zeta }^{{z,0,s + 1}}}$, ${{\zeta }^{{z,z,s + 1}}}$, чтобы получить двухточечную систему для величин $\zeta _{{k \pm 1/2}}^{{0,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{0,z,s + 1}}$ (см. [19]). Исключая из первого и третьего уравнений подсистемы (3.27) посредством дополнительных уравнений (2.21) величины ${{\zeta }^{{z,0,s + 1}}}$ и ${{\zeta }^{{z,z,s + 1}}}$, получаем следующую систему линейных уравнений относительно ${{\zeta }^{{0,0,s + 1}}}$ и ${{\zeta }^{{0,z,s + 1}}}$:

где Из системы (3.29) может быть найдена явная зависимость величин ${{\zeta }^{{0,0}}}$ и ${{\zeta }^{{0,z}}}$ от $\zeta _{{k \pm 1/2}}^{{0,0}}$ и $\zeta _{{k \pm 1/2}}^{{0,z}}$: где коэффициенты $x_{1}^{i}$ и $x_{2}^{i}$, $i = \overline {1,5} $, находятся путем решения пяти систем из двух линейных уравнений для каждой пространственной ячейки расчетной области: Массивы коэффициентов $x_{1}^{i}$ и $x_{2}^{k}$, $k = \overline {1,5} $, хранятся. Подставляя уравнения (3.30) (отнесенные к слою $s + 1$) в первое и третье уравнения подсистемы (3.25), получаем следующую двухточечную систему из двух уравнений для определения величин $\zeta _{{k \pm 1/2}}^{{0,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{0,z,s + 1}}$: где Для каждого значения $i$ система (3.32), (3.33) решается методом прогонки с граничными условиями на границах $z = {{z}_{{{\text{bot}}}}}$ и $z = {{z}_{{{\text{top}}}}}$ области:

После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{k \pm 1/2}}^{{0,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{0,z,s + 1}}$, поправки ${{\zeta }^{{0,0,s + 1}}}$, ${{\zeta }^{{0,z,s + 1}}}$, ${{\zeta }^{{z,0,s + 1}}}$ и ${{\zeta }^{{z,z,s + 1}}}$ для средних значений потока вычисляются по явным формулам (3.30) и (2.21). Величины ${{\zeta }^{{0,r,s + 1}}}$ и ${{\zeta }^{{z,r,s + 1}}}$ находятся по явным формулам из вторых уравнений подсистем (3.25) и (3.27):

Вторая подсистема для величин ${{\zeta }^{{r,0,s + 1}}}$, ${{\zeta }^{{r,z,s + 1}}}$, $\zeta _{{k \pm 1/2}}^{{r,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{r,z,s + 1}}$ состоит из первого и третьего уравнений в подсистеме уравнений (3.26), дополнительных уравнений (2.24) и граничных условий (3.28). Для ее решения исключим из нее посредством дополнительных уравнений (2.24) величины ${{\zeta }^{{r,0,s + 1}}}$ и ${{\zeta }^{{r,z,s + 1}}}$, чтобы получить двухточечную систему для величин $\zeta _{{k \pm 1/2}}^{{r,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{r,z,s + 1}}$:

где Для каждого значения индекса $i$ система (3.37), (3.38) решается методом прогонки (3.34) с граничными условиями (3.35). После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{k \pm 1/2}}^{{r,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{r,z,s + 1}}$, поправки ${{\zeta }^{{r,0,s + 1}}}$ и ${{\zeta }^{{r,z,s + 1}}}$ вычисляются по явным формулам (2.24). Величина ${{\zeta }^{{r,r,s + 1}}}$ находится по явной формуле из второго уравнения подсистемы (3.26): Для расчета величин $f_{{k \pm 1/2}}^{{\beta ,0,s + 1}}$ и $f_{{k \pm 1/2}}^{{\beta ,z,s + 1}}$, $\beta = 0,r$, в соответствии с уравнением (3.7), используются соотношения Далее, согласно уравнению (3.7), вычисляются величины ${{f}^{{\alpha ,s + 1}}}$ по формулам

Для расчета величин $f_{{i \pm 1/2}}^{{\beta ,\alpha ,s + 1}}$, $\beta = 0,z$, $\alpha = 0,r$, используемых, в частности, на следующем шаге алгоритма МР, воспользуемся системами дополнительных уравнений (2.20) и (2.23). Для определения величин $f_{{i \pm 1/2}}^{{0,0,s + 1}}$ и $f_{{i \pm 1/2}}^{{0,r,s + 1}}$ решается система, следующая из дополнительных уравнений (2.20):

где Для каждой плоскости $(j,k)$ система (3.42), (3.43) решается методом прогонки: К рекуррентным соотношениям (3.44) следует также добавить граничные условия на внешней $r = {{r}_{{{\text{ext}}}}}$ и внутренней $r = {{r}_{{\operatorname{int} }}}$ границах области:

Для определения величин $f_{{i \pm 1/2}}^{{z,0,s + 1}}$ и $f_{{i \pm 1/2}}^{{z,r,s + 1}}$ решается система, следующая из дополнительных уравнений (2.23):

где Для каждой плоскости $(j,k)$ система (3.46), (3.47) решается методом прогонки (3.44) с граничными условиями (3.44).

Эффективность МР весьма существенно зависит от выбора параметров $\omega _{r}^{{(s + 1)}}$ и $\omega _{z}^{{(s + 1)}}$. В отличие от метода переменных направлений для двумерной геометрии (см. [1], [26], [27]), теоретически обоснованный оптимальный выбор этих параметров отсутствует даже для случая, когда операторы $\hat {R}$ и $\hat {Z}$ являются самосопряженными, имеют положительный (отрицательный) спектр, а также перестановочны.

Однако в рассматриваемом нами общем случае гетерогенной задачи операторы $\hat {R}$ и $\hat {Z}$, вообще говоря, не перестановочны. Следовательно, результаты общей теории по выбору параметров $\omega _{r}^{{(s + 1)}}$ и $\omega _{z}^{{(s + 1)}}$ в данном случае будут применимы лишь приближенно. Вместе с тем можно попытаться воспользоваться для выбора указанных параметров алгоритмом, аналогичным использованному ранее в симметризованном ADI алгоритме (см. [16]), который формально применим и в рассматриваемом нами случае 2D геометрии.

В следующих разделах мы рассмотрим алгоритм для оценки минимального и максимального собственных значений операторов $\hat {R}$ и $\hat {Z}$, используемых в этом подходе, и приведем соотношения для расчета итерационных параметров.

4. ОЦЕНКА ГРАНИЦ СПЕКТРА РАДИАЛЬНОЙ И АКСИАЛЬНОЙ КОМПОНЕНТ ${{P}_{1}}$ ОПЕРАТОРА В $r,z$ ГЕОМЕТРИИ

Достаточно рассмотреть случай оператора $\hat {R}$, так как для оператора $\hat {Z}$ используется аналогичная процедура. Прежде всего, оператор $\hat {R}$ заменяется на оператор ${{\hat {R}}_{a}}$, полученный путем усреднения оператора $\hat {R}$ по переменной $z$. Для этого вычисляются средние по переменной $z$ значения коэффициентов ${{a}^{{r,\alpha }}}$, ${{b}^{{r,\alpha }}}$, ${{c}^{{r,\alpha }}}$ и ${{d}^{{r,\alpha }}}$, $\alpha = 0,r$, из систем (3.16)–(3.17), (3.22)–(3.23) в точке ${{\omega }_{r}} = 0$, величин ${{A}^{{r,0}}}$, ${{A}^{{r,1}}}$, ${{A}^{{r,2}}}$, ${{B}^{{r,0}}}$, ${{B}^{{r,1}}}$ и ${{B}^{{r,2}}}$ из уравнений (2.20), величин $d_{0}^{ + }$, $d_{0}^{ - }$, $d_{1}^{ + }$, $d_{1}^{ - }$, $d$, $e_{0}^{ + }$, $e_{0}^{ - }$, $e_{1}^{ + }$, $e_{1}^{ - }$ и $e$ из уравнений (3.14). Учитываются также коэффициенты дополнительных уравнений (2.23). Указанные величины определяют элементы матрицы оператора ${{\hat {R}}_{a}}$, который имеет блочную структуру и может быть представлен в виде

Здесь оператор ${{\hat {C}}_{0}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{0} = \operatorname{col} \{ g_{{i \pm 1/2}}^{{0,0}},g_{{i \pm 1/2}}^{{0,r}}\} $, определенный на гранях ${{r}_{{i \pm 1/2}}}$ ячеек, в вектор ${\mathbf{q}}_{i}^{0}$, состоящий из двух компонент вектора правой части ${\mathbf{q}}_{i}^{0} = \operatorname{col} \{ q_{i}^{{0,0}},q_{i}^{{0,r}}\} $, отнесенных к центрам ячеек ${{\hat {C}}_{0}}{\mathbf{g}}_{{}}^{0} = {\mathbf{q}}_{{}}^{0}$. Матрица оператора ${{\hat {C}}_{0}}$ получается усреднением по переменной $z$ (с весом ${{{v}}_{z}}$) матрицы системы (3.16). Оператор ${{\hat {G}}_{0}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{0}$ в вектор ${\mathbf{f}}_{i}^{0} = \operatorname{col} \{ f_{i}^{{0,0}},f_{i}^{{0,r}}\} $, определенный в центрах ячеек: ${\mathbf{f}}_{{}}^{0} = {{\hat {C}}_{0}}{\mathbf{g}}_{{}}^{0}$. Матрица оператора ${{\hat {G}}_{0}}$ состоит из усредненной по переменной $z$ (с весом ${{{v}}_{z}}$) матрицы системы (3.14). Оператор ${{\hat {G}}_{r}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{0}$ в вектор ${\mathbf{f}}_{i}^{r} = \operatorname{col} \{ f_{i}^{{r,0}},f_{i}^{{r,r}}\} $, определенный в центрах ячеек: ${\mathbf{f}}_{{}}^{r} = {{\hat {C}}_{r}}{\mathbf{g}}_{{}}^{0}$. Матрица оператора ${{\hat {G}}_{r}}$ состоит из усредненной по переменной $z$ (с весом ${{{v}}_{z}}$) матрицы системы (2.20) после подстановки в уравнения (2.20) величин ${{f}^{{0,0}}}$ и ${{f}^{{0,r}}}$, определяемых из усредненной по переменной $z$ матрицы системы (3.14).

Оператор ${{\hat {C}}_{z}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{z} = \operatorname{col} \{ g_{{i \pm 1/2}}^{{z,0}},g_{{i \pm 1/2}}^{{z,z}}\} $, определенный на гранях ${{r}_{{i \pm 1/2}}}$ ячеек, в вектор ${\mathbf{q}}_{i}^{z}$, состоящий из двух компонент вектора правой части ${\mathbf{q}}_{i}^{z} = \operatorname{col} \{ q_{i}^{{z,0}},q_{i}^{{z,z}}\} $, отнесенных к центрам ячеек: ${{\hat {C}}_{z}}{\mathbf{g}}_{{}}^{z} = {\mathbf{q}}_{{}}^{z}$. Матрица оператора ${{\hat {C}}_{z}}$ получается усреднением по переменной $z$ (с весом ${v}_{z}^{1}$) матрицы системы (3.22). Оператор ${{\hat {G}}_{z}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{z}$ в вектор ${\mathbf{f}}_{i}^{z} = \operatorname{col} \{ f_{i}^{{z,0}},f_{i}^{{z,z}}\} $, определенный в центрах ячеек: ${\mathbf{f}}_{{}}^{z} = {{\hat {C}}_{z}}{\mathbf{g}}_{{}}^{z}$. Матрица оператора ${{\hat {G}}_{z}}$ состоит из матрицы системы (2.23).

Для определения минимального и максимального собственных чисел ${{\lambda }_{r}}$ и ${{\Lambda }_{r}}$ оператора ${{\hat {R}}_{a}}$ может быть использован степенной метод (см. [26], [28]):

где Начальное приближение ${{{\mathbf{f}}}^{{(0)}}}$ также предполагается нормированным: $\left\| {{{{\mathbf{f}}}^{{(0)}}}} \right\| = 1$. (В качестве ${{{\mathbf{f}}}^{{(0)}}}$ можно выбрать, например, ${{{\mathbf{f}}}^{{(0)}}} = \frac{{{{{\mathbf{y}}}^{{(0)}}}}}{{\left\| {{{{\mathbf{y}}}^{{(0)}}}} \right\|}}$, где ${\mathbf{y}}_{i}^{{(0)}} = \operatorname{col} \{ 1,1,1,1,1,1\} $).

Для обращения операторов ${{\hat {G}}_{0}}$, ${{\hat {C}}_{0}}$, ${{\hat {G}}_{r}}$, ${{\hat {C}}_{z}}$, ${{\hat {G}}_{z}}$, входящих в состав оператора ${{\hat {R}}_{a}}$, используется прямой метод, аналогичный использованному при решении систем (2.20), (3.14), (3.16); (2.23), (3.22); (2.21), (3.30), (3.32); (2.24), (3.37).

При расчете максимального собственного значения по радиальной переменной на шаге $s + 1$ решается система уравнений, состоящая из двух подсистем.

Первая подсистема состоит из усредненных по переменной $z$ уравнений (2.20), (3.14) и (3.16). При этом, поскольку решается задача на собственное значение, в уравнениях (3.14) постоянные компоненты $d$ и $e$ не учитываются. Исключая из уравнений (2.20) величины ${{f}^{{0,0}}}$ и ${{f}^{{0,r}}}$ посредством соотношений (3.14), получаем двухточечную систему из двух уравнений для определения величин $f_{{i \pm 1/2}}^{{0,0,s + 1}}$ и $f_{{i \pm 1/2}}^{{0,r,s + 1}}$:

где Система (4.4), (4.5) решается методом прогонки (3.18) с граничными условиями (3.20). После ее решения компоненты собственного вектора $f_{i}^{{0,0,s + 1}}$ и $f_{i}^{{0,r,s + 1}}$ находятся по явным формулам с использованием соотношения (3.16) с усредненными по переменной $z$ коэффициентами: Далее определяются компоненты собственного вектора $f_{i}^{{r,0,s + 1}}$ и $f_{i}^{{r,r,s + 1}}$ с помощью явных уравнений (2.20). В $r,z$ геометрии для случая LD схемы для обеспечения сходимости степенного метода оценки максимального значения собственного значения оператора по радиальной переменной вес ${{P}_{r}}$ полагается равным $P_{r}^{0} = 0.01$.

Вторая подсистема состоит из уравнений (2.23) и (3.22). Соответствующая уравнениям (2.23) система из двух двухточечных уравнений для определения величин $f_{{i \pm 1/2}}^{{z,0,s + 1}}$ и $f_{{i \pm 1/2}}^{{z,r,s + 1}}$ имеет вид

где Система (4.7), (4.8) решается методом прогонки (3.18) с граничными условиями (3.20). После ее решения компоненты собственного вектора $f_{i}^{{z,0,s + 1}}$ и $f_{i}^{{z,r,s + 1}}$ находятся по явным формулам с использованием соотношения (2.23) с усредненными по переменной $z$ коэффициентами:

При расчете минимального собственного значения по радиальной переменной на шаге $s + 1$ решается система уравнений, состоящая из двух подсистем.

Первая подсистема состоит из усредненных по переменной $z$ уравнений (3.16):

и дополнительных уравнений (2.20) и (3.14). При этом, поскольку решается задача на собственное значение, в уравнениях (3.14) постоянные компоненты $d$ и $e$ не учитываются. Система (4.10) решается методом прогонки (3.18) с граничными условиями (3.20). После ее решения компоненты собственного вектора $f_{i}^{{r,0,s + 1}}$, $f_{i}^{{r,r,s + 1}}$, $f_{i}^{{0,0,s + 1}}$ и $f_{i}^{{0,r,s + 1}}$ находятся с помощью уравнений (2.20) и (3.14) (без компонент $d$ и $e$).

Вторая подсистема состоит из усредненных по переменной $z$ уравнений (3.22):

Система (4.11) решается методом прогонки (3.18) с граничными условиями (3.20). После ее решения компоненты собственного вектора $f_{i}^{{z,0,s + 1}}$ и $f_{i}^{{z,r,s + 1}}$ находятся с помощью уравнений (2.23).

Компоненты ${{f}^{{0,z}}}$, ${{f}^{{r,z}}}$ и ${{f}^{{z,z}}}$ собственных векторов, отвечающих минимальному и максимальному собственным значениям оператора ${{\hat {R}}_{a}}$, стремятся к нулю и их вкладом можно пренебречь.

При использовании адаптивной WLB-WLD схемы выбор весов схемы уточняется в процессе итераций. Поэтому границы спектра собственных значений операторов $\hat {R}$ и $\hat {Z}$ должны, вообще говоря, переоцениваться на каждой внутренней итерации $n$. В этом случае собственные значения и собственные вектора с предыдущей $n - 1$-й внутренней итерации используются как начальное приближение. В итоге, среднее число итераций степенного метода при точности расчета собственных значений ${{\varepsilon }_{\lambda }} = {{10}^{{ - 3}}}$ обычно не велико.

5. ОПРЕДЕЛЕНИЕ ОПТИМАЛЬНЫХ ПАРАМЕТРОВ МЕТОДА РАСЩЕПЛЕНИЯ

После оценки границ собственных значений операторов $\hat {R}$ и $\hat {Z}$ для циклического ADI метода с длиной цикла $J$ рассчитываются следующие величины (см. [27]):

где ${{\kappa }_{j}} = {{\kappa }_{j}}(\eta ,J)$ – известная функция (см. [27]) (в случае $J = 1$ $\kappa = \sqrt \eta $).

Величины ${{\omega }_{{r,j}}}$ и ${{\omega }_{{z,j}}}$ в уравнении (5.2) определяют искомые оптимальные значения ADI алгоритма.

В случае, когда ${{\lambda }_{r}} = {{\lambda }_{z}} = \lambda $, ${{\Lambda }_{r}} = {{\Lambda }_{z}} = \Lambda $, соотношения (5.1), (5.2) принимают вид

В используемом нами для случая 2D геометрии варианте МР полагается ${{\omega }_{{r,j}}} = {{\omega }_{{z,j}}} \equiv {{\omega }_{j}}$ с $\lambda = {{\lambda }_{{\min }}} = \min ({{\lambda }_{r}},{{\lambda }_{z}})$ и $\Lambda = {{\Lambda }_{{\max }}} = \max ({{\Lambda }_{r}},{{\Lambda }_{z}})$, ${{\tilde {\tau }}_{j}} = 2{{\omega }_{j}}$. Такой выбор шагов, как показывает численный эксперимент, является приемлемым для достаточно широкого класса задач. В качестве мотивации использования МР с симметризованными границами спектра отметим, что в гетерогенных задачах границы спектра оцениваются приближенно для усредненных значений операторов $\hat {R}$ и $\hat {Z}$. Использование симметризованной оценки границ спектра раздвигает эти границы, позволяя гасить более широкий спектр ошибок.

При выборе остальных параметров алгоритма: длины цикла $J$, критерия и точности сходимости итераций МР, точности расчета границ спектра операторов $\hat {R}$ и $\hat {Z}$ и т.д., позволяющих минимизировать время расчета варианта, использовался опыт задания аналогичных параметров для ADI метода в 2D геометрии (см. разд. 3.4 из [16]), а также результаты численного эксперимента, в котором в качестве критерия оптимизации $K{{P}_{1}}$ алгоритма выбиралось уменьшение полного времени расчета варианта. В частности, проведенное исследование показало, что использование поточечного критерия сходимости итераций МР для каждой из компонент ${{f}^{{\alpha ,0}}}$, ${{f}^{{\alpha ,r}}}$ и ${{f}^{{\alpha ,z}}}$, ${{P}_{1}}$ системы является чрезмерным. Более приемлемым является критерий сходимости вида

где ($\tilde {i},\tilde {k}$) – “грубая” сетка с коэффициентом огрубления $N$ (в 2D геометрии обычно $N = 4$) по отношению к сетке $\left( {i,k} \right)$. Блок ячеек, по которому проводится суммирование в (5.4), состоит из ${{N}^{2}}$ ячеек. В уравнении (5.4) ${{\varepsilon }_{{P1}}}$ – заданная точность сходимости итераций МР.

Численный эксперимент показал, что желательно использовать циклический МР с достаточно большой длиной цикла $J$ (хорошим выбором для большого класса задач является значение $J = 16$) с возрастающей последовательностью ${{\kappa }_{j}}$: ${{\kappa }_{j}} < {{\kappa }_{{j + 1}}}$. Известно (см. [4], [29]), что кинетическая часть итерации хорошо подавляет быстроменяющиеся по пространственной переменной компоненты в спектре ошибки. Основное назначение ${{P}_{1}}$ части итерации – это подавление медленно изменяющихся по пространственной переменной компонент спектра ошибки, что согласуется с выбором критерия сходимости итераций МР (5.4) и вышеуказанной последовательностью шагов в цикле МР.

Оптимальный выбор ${{\varepsilon }_{{P1}}}$ зависит от степени монотонности используемой разностной схемы и длины цикла $J$. Численный эксперимент показал, что при $J \geqslant 8$ для LD и LB схем в 2D расчете достаточно задать ${{\varepsilon }_{{P1}}} = 0.4$. Использование более высоких точностей ${{\varepsilon }_{{P1}}}$ может лишь незначительно уменьшить число внутренних итераций, но при этом возрастает число итераций МР и, соответственно, стоимость одной внутренней итерации и полное время расчета.

Численное исследование также показало, что приемлемая точность расчета границ спектра ${{\varepsilon }_{\lambda }} = {{10}^{{ - 3}}}$. Дальнейшее повышение точности расчета границ спектра практически не влияет на число итераций МР. Выход из итерационного цикла происходит либо по достижении заданной точности, либо если число итераций степенного метода превосходит некоторое предельное значение $N_{{\max }}^{\lambda }$ (обычно, $N_{{\max }}^{\lambda } = 200$).

Если число итераций МР превосходит некоторое предельное значение $N_{{\max }}^{{{\text{split}}}}$ (обычно $N_{{\max }}^{{{\text{split}}}} = \max \left( {32,J} \right)$), то осуществляется корректировка выбора шагов алгоритма МР – исходная верхняя граница спектра $\Lambda $ увеличивается в $10 \times {{2}^{{n - 1}}}$ раз, где n = 1, 2, … – номер корректировки границ спектра. После корректировки верхней границы спектра итерации МР продолжаются (с использованием ранее полученного приближения для поправок) до сходимости, если не возникает необходимости в дополнительной корректировке. Корректировка верхней границы спектра проводится до сходимости итераций МР, но не более 10 раз, после чего ускорение в группе отключается и используются чистые кинетические итерации.

Как показал численный эксперимент, в отличие от случая AWDD схемы, потребность в корректировке выбора шагов алгоритма МР при использовании LD или LB схем возникает только при решении специально подобранных задач.

Ускорение в группе отключается также по достижении некоторого максимально допустимого числа внутренних итераций (обычно, 50) в группе.

Для предотвращения возможных переполнений целесообразно также осуществлять проверку значения прогоночного коэффициента ${{\xi }^{0}}$ в уравнениях (3.18), (3.34) и (3.44). Если абсолютное значение этого коэффициента в некоторой ячейке существенно больше единицы (в 2D коде обычно используется критерий $\left| {{{\xi }^{0}}} \right| > 2.0$), то ускоряющая коррекция решения на текущей внутренней итерации не производится.

6. KP₁ СХЕМА УСКОРЕНИЯ ВНУТРЕННИХ ИТЕРАЦИЙ В $x,\;z$ ГЕОМЕТРИИ

Односкоростное уравнение переноса в $x,\;z$ геометрии имеет вид

Здесь $\xi $ и $\mu $ – направляющие косинусы единичного вектора $\Omega $ направления скорости частицы (1.2), который с учетом имеющейся в данной геометрии симметрии: $\psi (x,z,\mu ,\varphi ) = \psi (x,z,\mu ,2\pi - \varphi )$, изменяется в пределах полусферы $ - 1 \leqslant \xi $, $\mu \leqslant 1$, $0 \leqslant \varphi \leqslant \pi $. Правую часть уравнения (6.1), учитывая имеющуюся симметрию, можно представить в виде

Все расчетные формулы KR₁ схемы для $x,z$ геометрии получаются из соответствующих формул для $r,z$ геометрии путем замены

7. ЧИСЛЕННЫЕ ПРИМЕРЫ ИСПОЛЬЗОВАНИЯ $K{{P}_{1}}$ СХЕМЫ УСКОРЕНИЯ ВНУТРЕННИХ ИТЕРАЦИЙ В 2D ГЕОМЕТРИИ

В данном разделе мы приведем численные результаты, позволяющие оценить эффективность вышеописанного варианта согласованной $K{{P}_{1}}$ схемы. Использовался поточечный критерий сходимости внутренних итераций по скалярному потоку $\varepsilon = {{10}^{{ - 4}}}$. Для решения ${{P}_{1}}$ системы использовался циклический МР с длиной цикла $J = 16$.

При расчете тестовых задач 1–3 использовались LD и LB схемы с весом ${{P}_{\xi }} = 1$ в дополнительном уравнении (1.14) по угловой переменной, WLD и WLB схемы с весом ${{P}_{\xi }}$, выбираемым посредством AWDD схемы (см. [17], [22]), обеспечивающей положительность экстраполяции по угловой переменной. Использовался следующий вариант алгоритма выбора веса ${{P}_{\xi }}$:

где величина ${{u}_{\xi }}$ определялась в результате предварительного расчета ячейки с весом ${{P}_{\xi }} = 1$; $b_{\xi }^{'} = \max (1,2{{b}_{\xi }}U_{\xi }^{0})$, где ${{b}_{\xi }} \geqslant 1$ – параметр монотонизации по переменной $\xi $ (использовалось значение ${{b}_{\xi }} = 1$); корректирующая функция $P(U,\delta )$ выбиралась в виде

В качестве критерия сходимости итераций МР использовался поблочный критерий (5.4) с размером ребра блока $N = 4$ и ${{\varepsilon }_{{P1}}} = 0.4$. Расчетные времена приведены для ПК Intel Core i7-6950X.

Прежде всего, представим результаты использования $K{{P}_{1}}$ схемы для расчета двух тестовых задач из [5]. Геометрия и состав зон для этих задач для случая $x,\;y$ геометрии определены на фиг. 1 (задача 1, EIR-2) и фиг. 2 (задача 2, железо-водная композиция).

В табл. 1 приведено число внутренних итераций, требуемое для решения задач 1 и 2 без ускорения и при использовании $K{{P}_{1}}$ схемы в сочетании с циклическим МР для решения ${{P}_{1}}$ системы. В квадратных и круглых скобках приведены, соответственно, среднее (по внутренним итерациям) число итераций МР и процессорное время (мин); значение параметра $n \geqslant 1$ показывает, что задача считалась с увеличенной в $10 \times {{2}^{n}}$ начальной верхней границей спектра $\Lambda $ для достижения устойчивости (или повышения эффективности) МР.

В табл. 2 приведено суммарное (по 26 группам) число внутренних итераций при расчете модели быстрого реактора в $r,z$ геометрии (задача 3) на пространственной сетке (82 × 35) в ${{S}_{8}}{{P}_{1}}$ приближении при точности сходимости внутренних итераций $\varepsilon = {{10}^{{ - 3}}}$. В активной зоне был задан пространственно постоянный источник со спектром деления ${}^{{235}}U$. В круглых скобках приведено процессорное время (мин).

Применим далее рассматриваемый алгоритм к решению сильно-гетерогенной задачи с малым поглощением из [30] в $r,\;z$ геометрии (задача 4). Расчетная область этой задачи представляет собой гетерогенный цилиндр радиусом 25 см и высотой 50 см. Аксиальное сечение задачи изображено на фиг. 3. Внутренний цилиндрический канал заполнен материалом 1. Окружающая канал среда и центральный диск внутри канала заполнены материалом 2. Плотность материала 1 существенно меньше плотности материала 2. Нижний торец канала ($z = 0$, $0 \leqslant r \leqslant 5$ см) облучается изотропным, радиально постоянным источником единичной интенсивности. Имеется также внутренний изотропный источник с интенсивностью ${{10}^{{ - 6}}}$ частиц/(см³ с), распределенный равномерно по всей расчетной области. Для случая $x,z$ геометрии может быть определена аналогичная задача, расчетная область которой представляет собой гетерогенный бесконечный вдоль оси $Oy$ параллелепипед с каналом. Аксиальное сечение этого параллелепипеда в плоскости OXZ совпадает с аксиальным сечением гетерогенного цилиндра (фиг. 3).

В табл. 3 и 4 для некоторого набора значений полных сечений материалов 1 (${{\sigma }_{{t,1}}}$) и 2 (${{\sigma }_{{t,2}}}$), а также величин отношения сечения рассеяния к полному сечению $c = {{{{\sigma }_{s}}} \mathord{\left/ {\vphantom {{{{\sigma }_{s}}} {{{\sigma }_{t}}}}} \right. \kern-0em} {{{\sigma }_{t}}}}$, приведено число внутренних итераций при решении задачи 4 в $r,\;z$ и $x,\;z$ геометриях на пространственной сетке 50 × 100 шагов, с использованием шаговой (Step), AWDD, LD и LB схем в ${{S}_{6}}$ приближении с точностью сходимости итераций ${{10}^{{ - 5}}}$ при точности решения ${{P}_{1}}$ системы для ускоряющих поправок ${{\varepsilon }_{{P1}}}$ = 0.5 для Step и AWDD схем и ${{\varepsilon }_{{P1}}}$ = 0.4 для LD и LB схем.

Как следует из результатов, приведенных в табл. 3 и 4, в $r,\;z$ и $x,\;z$ геометриях эффективность согласованной $K{{P}_{1}}$ схемы ускорения падает по мере увеличения отношения ${{{{\sigma }_{{t2}}}} \mathord{\left/ {\vphantom {{{{\sigma }_{{t2}}}} {{{\sigma }_{{t1}}}}}} \right. \kern-0em} {{{\sigma }_{{t1}}}}}$ и увеличения доли рассеяния, хотя и остается на приемлемом уровне.

В качестве следующего примера (задача 5) приведем число внутренних и внешних итераций и расчетные времена (см. табл. 5) при расчете пространственного распределения потока нейтронов и фотонов в РУ ВВЭР-1000/320 в $r,\;z$ и $r,\;\vartheta $ геометриях, соответствующих продольному и поперечному сечениям расчетной области в области максимума потока нейтронов с энергией E > > 0.5 МэВ. Для задания исходной геометрии задачи методами комбинаторной геометрии использовался геометрический модуль программы MCU (см. [31]). Комбинаторное задание геометрии затем конвертировалось на сетку задачи c использованием метода трэйсинга (tracing) программой ConDat (см. [32]) с поддержанием локального баланса масс в пределах каждой пространственной ячейки за счет введения дополнительных смесей материалов в рамках volume fraction (VF) метода. Расчеты выполнены с использованием мультигрупповой библиотеки сечений V7-200N47G системы SCALE-6.1.2 (см. [33]) без первой группы нейтронов, доля спектра деления для которой в рассматриваемой задаче равна нулю. В качестве источника использовался потвэльный источник для первой половины 5-й кампании 3-го блока Балаковской АЭС. Расчет в $r,z$ геометрии выполнен с пространственной сеткой 190 × 245; в $r,\;\vartheta $ геометрии – для сектора поворотной симметрии 60° с пространственной сеткой 190 × 120. Для нейтронов с E > 3.0 МэВ использовалась квадратура Лебедева ${{L}_{{12,2}}}$ (см. [34]), для нейтронов с E < 3.0 МэВ и для фотонов – квадратура ${{L}_{{8,2}}}$. Использовался поточечный критерий сходимости по скалярному потоку внутренних и внешних (по области термализации) итераций $\varepsilon = {{10}^{{ - 4}}}$ и ${{\varepsilon }_{{{\text{upsc}}}}} = 5 \times {{10}^{{ - 3}}}$ соответственно; анизотропия рассеяния учитывалась в ${{P}_{3}}$ приближении.

Для ускорения внешних итераций по области термализации нейтронов использовалась $K{{P}_{1}}$ схема ускорения внешних итераций (см. [6], [13]). Спектральная зависимость ускоряющих поправок в этой схеме выбиралась на основе анализа Фурье основной гармоники ошибки итерационного метода Гаусса–Зейделя (см. [35]) для гомогенизированной задачи. Предложенный алгоритм решения ${{P}_{1}}$ системы для ускоряющих поправок использовался также и для решения аналогичной ${{P}_{1}}$ системы, возникающей на каждой внешней итерации в алгоритме ускорения внешних итераций. Для обеспечения поточечной сходимости внешних итераций во всех группах из области термализации нейтронов использование заданной точности сходимости внутренних итераций $\varepsilon $ оказывается недостаточным, поэтому по достижении точности сходимости внешних итераций ${{\varepsilon }_{{{\text{upsc}}}}} < 10$ точность сходимости внутренних итераций увеличивалась до значения $a\varepsilon $, где коэффициент $a = 2 \times {{10}^{{ - 2}}}$ и $a = 2 \times {{10}^{{ - 4}}}$ для AWDD и нодальных схем соответственно.

Приведем также результаты (см. фиг. 5), показывающие скорость сходимости ${{k}_{{{\text{eff}}}}}$ (${{k}_{{{\text{eff}}}}}$ определяется итерационно из решения задачи на собственное значение (см. [3], [36]) и определяет степень критичности исследуемой системы) в зависимости от выбора пространственной сетки и разностной схемы при фиксированной квадратуре $E{{S}_{8}}$ для одногрупповой двухзонной задачи из [20] в $r,z$ геометрии, изображенной на фиг. 4 (задача 6), и ее аналога в $x,z$ геометрии. В случае $x,z$ геометрии на левой границе расчетной области использовалось условие зеркального отражения. Данная задача решалась на равномерной пространственной сетке из 10, 20, 40, 80 и 160 интервалов по каждой из переменных. В качестве точного использовалось значение ${{k}_{{{\text{eff}}}}}$, рассчитанное по LB схеме на сетке из 320 × 320 интервалов с точностью сходимости итераций ${{10}^{{ - 10}}}$.

8. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Численный эксперимент показал, что МР в сочетании с AWLD схемой обычно сходится за 5–7 итераций и это число итераций практически не зависит от номера внутренней итерации и номера группы. Это существенно лучший результат по сравнению с использованием ADI алгоритма в сочетании с адаптивной WDD (AWDD) схемой. Использование AWLD схемы позволяет существенно уменьшить разностную ошибку аппроксимации задачи по пространственным переменным, а также оценить пространственное распределение ошибки аппроксимации при решении задачи с использованием AWDD схемы.

Небольшие расчетные времена позволяют использовать разработанный алгоритм для решения практических задач радиационной защиты с использованием детального энергетического разбиения, предоставляемого мультигрупповыми библиотеками V7-200N47G системы SCALE-6.1.2 и 299n+127g системы ABBN-RF/CONSYST.