Журнал вычислительной математики и математической физики, 2019, T. 59, № 5, стр. 796-821

2.1. Алгоритм решения ${{P}_{1}}$-системы для ускоряющих поправок в $r,\;\vartheta ,\;z$ геометрии

${{P}_{1}}$-система для ускоряющих поправок для $K{{P}_{1}}$-схемы в $r,\;\vartheta ,\;z$ геометрии состоит из 16 балансных уравнений (1.2.19)–(1.2.22), шести дополнительных уравнений (1.2.26), (1.2.27) и (1.2.28); 16 дополнительных уравнений (1.2.30)–(1.2.35) и граничных условий (1.2.36)–(1.2.43) работы [1]. Эта система может быть решена итерационно с использованием, например, метода расщепления (МР) [2], [11].

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

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

Следует подчеркнуть, что исключая из оператора $\hat {A}$ с помощью дополнительных уравнений (1.2.26)–(1.2.28), (1.2.30)–(1.2.35) и граничных условий (1.2.36), (1.2.38) и (1.2.42) поправки, относящиеся к граням ячеек, можно представить его в виде, в котором присутствуют только значения поправок, отнесенные к центру ячейки. Соответственно в (2.1.7) можно считать, что единичный оператор $\hat {E}$ действует только на поправки, отнесенные к центру ячейки.

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

Отметим, что в системе уравнений (2.1.8)–(2.1.11) и ниже, как и в реализации МР для WDD-схемы [10], используется симметричная декомпозиция членов $ \sim {\kern 1pt} {{\sigma }^{{00}}}$ и ${{\sigma }^{{11}}}$ по операторам $\hat {R}$, $\hat {\Theta }$ и $\hat {Z}$ (шагам $\tau _{r}^{{(s + 1)}}$, $\tau _{\vartheta }^{{(s + 1)}}$ и $\tau _{z}^{{(s + 1)}}$). Это позволяет сохранить вид операторов $\hat {R}$, $\hat {\Theta }$ и $\hat {Z}$ подобным виду P₁ оператора для ускоряющих поправок в одномерной геометрии [13], [14].

Далее мы будем говорить, что величина относится к слою $s$, $s + 1{\text{/}}3$ и т.д., если она имеет соответствующий индекс. Следует отметить, что в рассматриваемом нами случае 3D геометрии величины $\zeta $ для дробных шагов $s + 1{\text{/}}3$ и $s + 2{\text{/}}3$ не аппроксимируют исходную нестационарную задачу (2.1.5) [15]. Для получения граничных условий по радиальной переменной для поправок $\zeta $ отметим, что из граничных условий (1.2.36) для величин $f$ в предположении (1.2.44), а также уравнения (2.1.7) следует, что:

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

Система (2.1.8)–(2.1.11) для 16 уравнений, дополненная двумя дополнительными уравнениями для величин ${{\zeta }^{{r,0,s + 1/3}}}$и ${{\zeta }^{{r,r,s + 1/3}}}$, аналогичным (1.2.26), и четырьмя дополнительными уравнениями (1.2.30) и (1.2.31) для величин ${{\zeta }^{{\vartheta ,0,s + 1/3}}}$, ${{\zeta }^{{\vartheta ,r,s + 1/3}}}$, ${{\zeta }^{{z,0,s + 1/3}}}$ и ${{\zeta }^{{z,r,s + 1/3}}}$, а также граничными условиями (2.1.13), может быть решена прямым методом в предположении, что величины с индексом $s$ известны с предыдущего шага. Данная система состоит из трех подсистем, которые могут быть решены независимо методом прогонки [12], [2], [8] по переменной $r$, а также 8 уравнений – последних двух из подсистем уравнений (2.1.8)–(2.1.11), которые могут быть явно разрешены относительно переменных ${{\zeta }^{{0,\vartheta ,s + 1/3}}}$, ${{\zeta }^{{0,z,s + 1/3}}}$, ${{\zeta }^{{r,\vartheta ,s + 1/3}}}$, ${{\zeta }^{{r,z,s + 1/3}}}$, ${{\zeta }^{{\vartheta ,\vartheta ,s + 1/3}}}$, ${{\zeta }^{{\vartheta ,z,s + 1/3}}}$, ${{\zeta }^{{z,\vartheta ,s + 1/3}}}$ и ${{\zeta }^{{z,z,s + 1/3}}}$.

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

где Из системы (2.1.14) может быть найдена явная зависимость величин ${{\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} $, хранятся. Подставляя уравнения (2.1.15) (отнесенных к слою $s + 1{\text{/}}3$) в первые два из уравнений подсистемы (2.1.8), получаем следующую двухточечную систему из двух уравнений для определения величин $\zeta _{{i \pm 1/2,j,k}}^{{0,s + 1/3}}$ и $\zeta _{{i \pm 1/2,j,k}}^{{r,s + 1/3}}$: где Для каждой плоскости $(j,k)$ система (2.1.17)–(2.1.18) решается методом прогонки [12], [2], [8]: где коэффициенты ${{\xi }_{i}} = {\text{col}}\{ \xi _{i}^{0},\xi _{i}^{r}\} $ и ${{\eta }_{i}} = {\text{col}}\{ \eta _{i}^{0},\eta _{i}^{r}\} $ вычисляются посредством рекуррентных соотношений: К рекуррентным соотношениям (2.1.19), (2.1.20) следует также добавить граничные условия на внешней $r = {{r}_{{{\text{ext}}}}}$ и внутренней $r = {{r}_{{\operatorname{int} }}}$ границах области: Достаточным условием устойчивости этого алгоритма прогонки является выполнение неравенства: $\left| {\xi _{i}^{0}} \right| < 1$. После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{i \pm 1/2}}^{{0,r,s + 1/3}}$ и $\zeta _{{i \pm 1/2}}^{{0,0,s + 1/3}}$, поправки ${{\zeta }^{{0,0,s + 1/3}}}$, ${{\zeta }^{{0,r,s + 1/3}}}$, ${{\zeta }^{{r,0,s + 1/3}}}$ и ${{\zeta }^{{r,r,s + 1/3}}}$ для средних значений потока вычисляются по явным формулам (2.1.15) и (1.2.26). Величины ${{\zeta }^{{0,\vartheta ,s + 1/3}}}$, ${{\zeta }^{{0,z,s + 1/3}}}$, ${{\zeta }^{{r,\vartheta ,s + 1/3}}}$ и ${{\zeta }^{{r,z,s + 1/3}}}$ находятся по явным формулам из последних двух уравнений подсистем (2.1.8) и (2.1.9):

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

где Для каждой плоскости $(j,k)$ система (2.1.23), (2.1.24) решается методом прогонки (2.1.19) с граничными условиями (2.1.21). После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{i \pm 1/2}}^{{\vartheta ,0,s + 1/3}}$ и $\zeta _{{i \pm 1/2}}^{{\vartheta ,r,s + 1/3}}$, поправки ${{\zeta }^{{\vartheta ,0,s + 1/3}}}$ и ${{\zeta }^{{\vartheta ,r,s + 1/3}}}$ вычисляются по явным формулам (1.2.30). Величины ${{\zeta }^{{\vartheta ,\vartheta ,s + 1/3}}}$ и ${{\zeta }^{{\vartheta ,z,s + 1/3}}}$ находятся по явным формулам из последних двух уравнений подсистемы (2.1.10):

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

где Для каждой плоскости $(j,k)$ система (2.1.26), (2.1.27) решается методом прогонки (2.1.19) с граничными условиями (2.1.21). После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{i \pm 1/2}}^{{z,0,s + 1/3}}$ и $\zeta _{{i \pm 1/2}}^{{z,r,s + 1/3}}$, поправки ${{\zeta }^{{z,0,s + 1/3}}}$ и ${{\zeta }^{{z,r,s + 1/3}}}$ вычисляются по явным формулам (1.2.31). Величины ${{\zeta }^{{z,\vartheta ,s + 1/3}}}$ и ${{\zeta }^{{z,z,s + 1/3}}}$ находятся по явным формулам из последних двух уравнений подсистемы (2.1.11): На 2-м шаге $\tau _{\vartheta }^{{(s + 1)}}$ алгоритма МР решается следующая система: Граничные условия для системы (2.1.29)–(2.1.32) имеют вид

Система (2.1.29)–(2.1.32) для 16 уравнений, дополненная двумя дополнительными уравнениями для величин ${{\zeta }^{{\vartheta ,0,s + 2/3}}}$и ${{\zeta }^{{\vartheta ,\vartheta ,s + 2/3}}}$, аналогичным (1.2.27), и четырьмя дополнительными уравнениями (1.2.32) и (1.2.33) для величин ${{\zeta }^{{r,0,s + 2/3}}}$, ${{\zeta }^{{r,\vartheta ,s + 2/3}}}$, ${{\zeta }^{{z,0,s + 2/3}}}$ и ${{\zeta }^{{z,\vartheta ,s + 2/3}}}$, а также граничными условиями (2.1.33), может быть решена прямым методом в предположении, что величины с индексом $s + 1{\text{/}}3$ известны с предыдущего шага. Данная система состоит из трех подсистем, которые могут быть решены независимо методом прогонки по переменной $\vartheta $, а также 8 уравнений – второго и четвертого из подсистем уравнений (2.1.29)–(2.1.32), которые могут быть явно разрешены относительно переменных ${{\zeta }^{{0,r,s + 2/3}}}$, ${{\zeta }^{{0,z,s + 2/3}}}$, ${{\zeta }^{{r,r,s + 2/3}}}$, ${{\zeta }^{{r,z,s + 2/3}}}$, ${{\zeta }^{{\vartheta ,r,s + 2/3}}}$, ${{\zeta }^{{\vartheta ,z,s + 2/3}}}$, ${{\zeta }^{{z,r,s + 2/3}}}$ и ${{\zeta }^{{z,z,s + 2/3}}}$.

Первая подсистема относительно величин ${{\zeta }^{{0,0,s + 2/3}}}$, ${{\zeta }^{{0,\vartheta ,s + 2/3}}}$, ${{\zeta }^{{\vartheta ,0,s + 2/3}}}$, ${{\zeta }^{{\vartheta ,\vartheta ,s + 2/3}}}$, $\zeta _{{j \pm 1/2}}^{{0,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{0,\vartheta ,s + 2/3}}$ состоит из первого и третьего уравнений в подсистемах уравнений (2.1.29) и (2.1.31) (всего 4 уравнения), дополнительных уравнений (1.2.27) для величин ${{\zeta }^{{\vartheta ,0,s + 2/3}}}$и ${{\zeta }^{{\vartheta ,\vartheta ,s + 2/3}}}$ и граничных условий (2.1.33). Для ее решения мы должны исключить из нее величины ${{\zeta }^{{0,0,s + 2/3}}}$, ${{\zeta }^{{0,\vartheta ,s + 2/3}}}$, ${{\zeta }^{{\vartheta ,0,s + 2/3}}}$, ${{\zeta }^{{\vartheta ,\vartheta ,s + 2/3}}}$, чтобы получить двухточечную систему для величин $\zeta _{{j \pm 1/2}}^{{0,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{0,\vartheta ,s + 2/3}}$ [14]. Исключая из первого и третьего уравнений подсистемы (2.1.31) посредством дополнительных уравнений (1.2.27) величины ${{\zeta }^{{\vartheta ,0,s + 2/3}}}$и ${{\zeta }^{{\vartheta ,\vartheta ,s + 2/3}}}$ получим следующую систему линейных уравнений относительно ${{\zeta }^{{0,0,s + 2/3}}}$ и ${{\zeta }^{{0,\vartheta ,s + 2/3}}}$:

где Из системы (2.1.14) может быть найдена явная зависимость величин ${{\zeta }^{{0,0}}}$ и ${{\zeta }^{{0,\vartheta }}}$ от $\zeta _{{j \pm 1/2}}^{{0,0}}$ и $\zeta _{{j \pm 1/2}}^{{0,\vartheta }}$: где коэффициенты $x_{1}^{k}$ и $x_{2}^{k}$, $k = \overline {1,5} $, находятся путем решения пяти систем из двух линейных уравнений для каждой пространственной ячейки расчетной области: Массивы коэффициентов $x_{1}^{k}$ и $x_{2}^{k}$, $k = \overline {1,5} $, хранятся. Подставляя уравнения (2.1.35) (отнесенные к слою $s + 2{\text{/}}3$) в первое и третье уравнения подсистемы (2.1.29), получим следующую двухточечную систему из двух уравнений для определения величин $\zeta _{{j \pm 1/2}}^{{0,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{0,\vartheta ,s + 2/3}}$: где Для каждой плоскости $(i,k)$ система (2.1.37), (2.1.38) решается методом прогонки: с граничными условиями на границах $\vartheta = {{\vartheta }_{0}}$ и $\vartheta = {{\vartheta }_{{end}}}$ области:

При наличии периодического граничного условия на границах $\vartheta = {{\vartheta }_{0}}$ и $\vartheta = {{\vartheta }_{{end}}}$ (1.2.45) вместо граничных условий (2.1.33) следует воспользоваться периодическими граничными условиями для поправок:

Для решения системы (2.1.37), (2.1.41) воспользуемся методом циклической прогонки. В [17] описан метод циклической прогонки для решения трехточечного уравнения при наличии периодических граничных условий. Воспользуемся его обобщением [16], [10] для решения двухточечной системы (2.1.37), (2.1.41). Для упрощения обозначений ниже индексы $\beta $ и $s + 2{\text{/}}3$ будем опускать. Будем искать решение системы (2.1.37), (2.1.41) в виде линейной комбинации решения системы (2.1.37) с правой частью $q_{{i,j,k}}^{{\vartheta ,\alpha }}$ с граничными условиями зеркального отражения: которое обозначим $\zeta _{{i,j,k}}^{{\vartheta ,a}}$, и решения этой же системы с нулевой правой частью, удовлетворяющее граничным условиям которое мы обозначим $\zeta _{{i,j,k}}^{{\vartheta ,b}}$. Искомое решение имеет вид: где $\lambda $ – параметр, выбираемый из условия Подстановка уравнения (2.1.44) в уравнение (2.1.45) дает для параметра $\lambda $ следующее значение:

Таким образом, для получения искомого решения системы (2.1.37) с граничными условиями (2.1.41) методом циклической прогонки систему (2.1.37) приходится решать дважды с граничными условиями (2.1.42) и (2.1.43) с использованием обычной прогонки для двухточечных систем, а искомое решение получается в результате линейной комбинации (2.1.44) полученных решений. Это приводит к некоторому увеличению объема вычислений по сравнению со случаем наличия условий зеркального отражения, когда достаточно воспользоваться обычной двухточечной прогонкой.

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

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

где Здесь величины $\zeta _{{j \pm 1/2}}^{{0,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{0,\vartheta ,s + 2/3}}$ берутся из решения 1-й подсистемы. Для каждой плоскости $(i,k)$ система (2.1.48), (2.1.49) решается методом прогонки (2.1.39) с граничными условиями (2.1.40). После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{j \pm 1/2}}^{{r,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{r,\vartheta ,s + 2/3}}$, поправки ${{\zeta }^{{r,0,s + 2/3}}}$ и ${{\zeta }^{{r,\vartheta ,s + 2/3}}}$ вычисляются по явным формулам (1.2.32). Величины ${{\zeta }^{{r,r,s + 2/3}}}$ и ${{\zeta }^{{r,z,s + 2/3}}}$ находятся по явным формулам из второго и четвертого уравнений подсистемы (2.1.30):

Третья подсистема для величин ${{\zeta }^{{z,0,s + 2/3}}}$, ${{\zeta }^{{z,\vartheta ,s + 2/3}}}$, $\zeta _{{j \pm 1/2}}^{{z,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{z,\vartheta ,s + 2/3}}$ состоит из первого и третьего уравнений в подсистеме уравнений (2.1.32), дополнительных уравнений (1.2.33) и граничных условий (2.1.33). Для ее решения исключим из нее посредством дополнительных уравнений (1.2.33) величины ${{\zeta }^{{z,0,s + 2/3}}}$ и ${{\zeta }^{{z,\vartheta ,s + 2/3}}}$, чтобы получить двухточечную систему для величин $\zeta _{{j \pm 1/2}}^{{z,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{z,\vartheta ,s + 2/3}}$:

где Для каждой плоскости $(i,k)$ система (2.1.51), (2.1.52) решается методом прогонки (2.1.39) с граничными условиями (2.1.40). После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{j \pm 1/2}}^{{z,0,s + 2/3}}$ и $\zeta _{{j \pm 1/2}}^{{z,\vartheta ,s + 2/3}}$, поправки ${{\zeta }^{{z,0,s + 2/3}}}$ и ${{\zeta }^{{z,\vartheta ,s + 2/3}}}$ вычисляются по явным формулам (1.2.33). Величины ${{\zeta }^{{z,r,s + 2/3}}}$ и ${{\zeta }^{{z,z,s + 2/3}}}$ находятся по явным формулам из второго и четвертого уравнений подсистемы (2.1.32): На 3-м шаге $\tau _{z}^{{(s + 1)}}$ алгоритма МР решается следующая система: Граничные условия для системы (2.1.54)–(2.1.57), с учетом предположения (1.2.44) и последнего из уравнений алгоритма (2.1.7), имеют вид

Система (2.1.54)–(2.1.57) для 16 уравнений, дополненная двумя дополнительными уравнениями для величин ${{\zeta }^{{z,0,s + 1}}}$и ${{\zeta }^{{z,z,s + 1}}}$, аналогичным (1.2.28), и четырьмя дополнительными уравнениями  (1.2.34)  и (1.2.35) для величин ${{\zeta }^{{r,0,s + 1}}}$, ${{\zeta }^{{r,z,s + 1}}}$, ${{\zeta }^{{\vartheta ,0,s + 1}}}$ и ${{\zeta }^{{\vartheta ,z,s + 1}}}$, а также граничными условиями (2.1.58), может быть решена прямым методом в предположении, что величины с индексом $s + 2{\text{/}}3$ известны с предыдущего шага. Данная система состоит из трех подсистем, которые могут быть решены независимо методом прогонки по переменной $z$, а также 8 уравнений – второго и третьего из подсистем уравнений (2.1.54)–(2.1.57), которые могут быть явно разрешены относительно переменных ${{\zeta }^{{0,r,s + 1}}}$, ${{\zeta }^{{0,\vartheta ,s + 1}}}$, ${{\zeta }^{{r,r,s + 1}}}$, ${{\zeta }^{{r,\vartheta ,s + 1}}}$, ${{\zeta }^{{\vartheta ,r,s + 1}}}$, ${{\zeta }^{{\vartheta ,\vartheta ,s + 1}}}$, ${{\zeta }^{{z,r,s + 1}}}$ и ${{\zeta }^{{z,\vartheta ,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}}$ состоит из первого и четвертого уравнений в подсистемах уравнений (2.1.54) и (2.1.57) (всего 4 уравнения), дополнительных уравнений (1.2.28) для величин ${{\zeta }^{{z,0,s + 1}}}$ и ${{\zeta }^{{z,z,s + 1}}}$ и граничных условий (2.1.58). Для ее решения мы должны исключить из нее величины ${{\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}}$ [14]. Исключая из первого и четвертого уравнений подсистемы (2.1.57) посредством дополнительных уравнений (1.2.28) величины ${{\zeta }^{{z,0,s + 1}}}$ и ${{\zeta }^{{z,z,s + 1}}}$, получаем следующую систему линейных уравнений относительно ${{\zeta }^{{0,0,s + 1}}}$ и ${{\zeta }^{{0,z,s + 1}}}$:

где Из системы (2.1.59) может быть найдена явная зависимость величин ${{\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} $, хранятся. Подставляя уравнения (2.1.60) (отнесенные к слою $s + 1$) в первое и четвертое уравнения подсистемы (2.1.54), получим следующую двухточечную систему из двух уравнений для определения величин $\zeta _{{k \pm 1/2}}^{{0,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{0,z,s + 1}}$: где Для каждой плоскости $(i,j)$ система (2.1.62), (2.1.63) решается методом прогонки: с граничными условиями на границах $z = {{z}_{{bot}}}$ и $z = {{z}_{{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}}}$ для средних значений потока вычисляются по явным формулам (2.1.60) и (1.2.28). Величины ${{\zeta }^{{0,r,s + 1}}}$, ${{\zeta }^{{0,\vartheta ,s + 1}}}$, ${{\zeta }^{{z,r,s + 1}}}$ и ${{\zeta }^{{z,\vartheta ,s + 1}}}$ находятся по явным формулам из второго и третьего уравнений подсистем (2.1.54) и (2.1.57):

Вторая подсистема для величин ${{\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}}$ состоит из первого и четвертого уравнений в подсистеме уравнений (2.1.55), дополнительных уравнений (1.2.34) и граничных условий (2.1.58).  Для  ее  решения  исключим из нее посредством дополнительных уравнений (1.2.34) величины ${{\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,j)$ система (2.1.67), (2.1.68) решается методом прогонки (2.1.64) с граничными условиями (2.1.65). После нахождения величин поправок для потоков на гранях ячейки $\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}}}$ вычисляются по явным формулам (1.2.34). Величины ${{\zeta }^{{r,r,s + 1}}}$ и ${{\zeta }^{{r,\vartheta ,s + 1}}}$ находятся по явным формулам из второго и третьего уравнений подсистемы (2.1.55):

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

где Для каждой плоскости $(i,j)$ система (2.1.70), (2.1.71) решается методом прогонки (2.1.64) с граничными условиями (2.1.65). После нахождения величин поправок для потоков на гранях ячейки $\zeta _{{k \pm 1/2}}^{{\vartheta ,0,s + 1}}$ и $\zeta _{{k \pm 1/2}}^{{\vartheta ,z,s + 1}}$, поправки ${{\zeta }^{{\vartheta ,0,s + 1}}}$ и ${{\zeta }^{{\vartheta ,z,s + 1}}}$ вычисляются по явным формулам (1.2.35). Величины ${{\zeta }^{{\vartheta ,r,s + 1}}}$ и ${{\zeta }^{{\vartheta ,\vartheta ,s + 1}}}$ находятся по явным формулам из второго и третьего уравнений подсистемы (2.1.56): Для расчета величин $f_{{k \pm 1/2}}^{{\beta ,0,s + 1}}$ и $f_{{k \pm 1/2}}^{{\beta ,z,s + 1}}$, $\beta = 0,r,\vartheta $, в соответствии с уравнением (2.1.7), используются соотношения: Далее, согласно уравнению (2.1.7), вычисляются величины ${{f}^{{\alpha ,s + 1}}}$ по формулам:

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

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

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

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

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

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

Аналогично, для определения величин $f_{{j \pm 1/2}}^{{0,0,s + 1}}$ и $f_{{j \pm 1/2}}^{{0,\vartheta ,s + 1}}$ решается система, задаваемая дополнительными уравнениями (1.2.27):

где Для каждой плоскости $(i,k)$ система (2.1.83), (2.1.84) решается методом прогонки: с граничными условиями на границах $\vartheta = {{\vartheta }_{0}}$ и $\vartheta = {{\vartheta }_{{end}}}$ области:

При использовании периодических граничных условий:

для решения системы (2.1.83), (2.184) используется метод циклической прогонки.

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

где Для каждой плоскости $(i,k)$ система (2.1.88), (2.1.89) решается методом прогонки (2.1.85) с граничными условиями (2.1.86).

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

где Для каждой плоскости $(i,k)$ система (2.1.90), (2.1.91) решается методом прогонки (2.1.85) с граничными условиями (2.1.86).

Следует отметить, что в случае использования периодического граничного решения по переменной $\vartheta $ перенормировка решения с использованием ускоряющих поправок, в соответствии с уравнениями (1.2.17), должна производиться на обеих границах $\vartheta = {{\vartheta }_{0}}$ и $\vartheta = {{\vartheta }_{{{\text{end}}}}}$ расчетной области:

В случае же использования условий зеркального отражения на границах $\vartheta = {{\vartheta }_{0}}$ и $\vartheta = {{\vartheta }_{{{\text{end}}}}}$ расчетной области достаточно произвести перенормировку решения только на границе $\vartheta = {{\vartheta }_{{{\text{end}}}}}$ согласно второму из уравнений (2.1.92).

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

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

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

2.2. Оценка границ спектра $r$, $\vartheta $ и $z$ компонент ${{P}_{1}}$ оператора

Достаточно рассмотреть случай оператора $\hat {R}$, так как для операторов $\hat {\Theta }$ и $\hat {Z}$ используется аналогичная процедура. Прежде всего оператор $\hat {R}$ заменяется на оператор ${{\hat {R}}_{a}}$, полученный путем усреднения оператора $\hat {R}$ по переменным $\vartheta $ и $z$. Для этого вычисляются средние по переменным $\vartheta $ и $z$ значения коэффициентов ${{a}^{{r,\alpha }}}$, ${{b}^{{r,\alpha }}}$, ${{c}^{{r,\alpha }}}$ и ${{d}^{{r,\alpha }}}$, $\alpha = 0,r$, из систем (2.1.17), (2.1.18), (2.1.23), (2.1.24) и (2.1.26), (2.1.27) в точке ${{\omega }_{r}} = 0$, величин ${{A}^{{r,0}}}$, ${{A}^{{r,1}}}$, ${{A}^{{r,2}}}$, ${{B}^{{r,0}}}$, ${{B}^{{r,1}}}$ и ${{B}^{{r,2}}}$ из уравнений (1.2.26), величин $d_{0}^{ + }$, $d_{0}^{ - }$, $d_{1}^{ + }$, $d_{1}^{ - }$, $d$, $e_{0}^{ + }$, $e_{0}^{ - }$, $e_{1}^{ + }$, $e_{1}^{ - }$ и $e$ из уравнений (2.1.15). Учитываются также коэффициенты дополнительных уравнений (1.2.30) и (1.2.31). Указанные величины определяют элементы матрицы оператора ${{\hat {R}}_{a}}$, который имеет блочную структуру и может быть представлен в виде

Здесь оператор ${{\hat {C}}_{0}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{0} = {\text{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} = {\text{col}}\{ q_{i}^{{0,0}},q_{i}^{{0,r}}\} $, отнесенных к центрам ячеек: ${{\hat {C}}_{0}}{{{\mathbf{g}}}^{0}} = {{{\mathbf{q}}}^{0}}$. Матрица оператора ${{\hat {C}}_{0}}$ получается усреднением по переменным $\vartheta $ и $z$ (с весом ${{v}_{\vartheta }}{{v}_{z}}$) матрицы системы (2.1.17). Оператор ${{\hat {G}}_{0}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{0}$ в вектор ${\mathbf{f}}_{i}^{0} = {\text{col}}\{ f_{i}^{{0,0}},f_{i}^{{0,r}}\} $, определенный в центрах ячеек: ${{{\mathbf{f}}}^{0}} = {{\hat {G}}_{0}}{{{\mathbf{g}}}^{0}}$. Матрица оператора ${{\hat {G}}_{0}}$ состоит из усредненной по переменным $\vartheta $ и $z$ (с весом ${{v}_{\vartheta }}{{v}_{z}}$) матрицы системы (2.1.15). Оператор ${{\hat {G}}_{r}}$ преобразует вектор ${\mathbf{g}}_{{i \pm 1/2}}^{0}$ в вектор ${\mathbf{f}}_{i}^{r} = {\text{col}}\{ f_{i}^{{r,0}},f_{i}^{{r,r}}\} $, определенный в центрах ячеек: ${{{\mathbf{f}}}^{r}} = {{\hat {G}}_{r}}{{{\mathbf{g}}}^{0}}$. Матрица оператора ${{\hat {G}}_{r}}$ состоит из усредненной по переменным $\vartheta $ и $z$ (с весом ${{v}_{\vartheta }}{{v}_{z}}$) матрицы системы (1.2.26) после подстановки в уравнения (1.2.26) величин ${{f}^{{0,0}}}$ и ${{f}^{{0,r}}}$, определяемых из усредненной по переменным $\vartheta $ и $z$ матрицы системы (2.1.15).

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

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

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

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

Для обращения операторов ${{\hat {G}}_{0}}$, ${{\hat {C}}_{0}}$, ${{\hat {G}}_{r}}$, ${{\hat {C}}_{\vartheta }}$, ${{\hat {G}}_{\vartheta }}$, ${{\hat {C}}_{z}}$, ${{\hat {G}}_{z}}$, входящих в состав оператора ${{\hat {R}}_{a}}$, используется прямой метод, аналогичный использованному при решении систем (1.2.26), (2.1.15)–(2.1.17); (1.2.30), (2.1.23); (1.2.31), (2.1.26).

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

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

где Система (2.2.4), (2.2.5) решается методом прогонки (2.1.77) с граничными условиями (2.1.78). После ее решения компоненты собственного вектора $f_{i}^{{0,0,s + 1}}$ и $f_{i}^{{0,r,s + 1}}$ находятся по явным формулам с использованием соотношения (2.1.17) с усредненными по переменным $z$ и $\vartheta $ коэффициентами:

Далее определяются компоненты собственного вектора $f_{i}^{{r,0,s + 1}}$ и $f_{i}^{{r,r,s + 1}}$ с помощью явных уравнений (1.2.26). В $r,\;\vartheta ,\;z$ геометрии для случая LD-схемы для обеспечения сходимости степенного метода оценки максимального значения собственного значения оператора по радиальной переменной вес ${{P}_{r}}$ полагается равным $P_{r}^{0} = 0.01$.

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

где Система (2.2.7), (2.2.8) решается методом прогонки (2.1.77) с граничными условиями (2.1.78). После ее решения компоненты собственного вектора $f_{i}^{{\vartheta ,0,s + 1}}$ и $f_{i}^{{\vartheta ,r,s + 1}}$ находятся по явным формулам с использованием соотношения (2.1.23) с усредненными по переменным $z$ и $\vartheta $ коэффициентами:

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

где Система (2.2.10), (2.2.11) решается методом прогонки (2.1.77) с граничными условиями (2.1.78). После ее решения компоненты собственного вектора $f_{i}^{{z,0,s + 1}}$ и $f_{i}^{{z,r,s + 1}}$ находятся по явным формулам с использованием соотношения (2.1.26) с усредненными по переменным $z$ и $\vartheta $ коэффициентами:

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

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

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

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

Система (2.2.14) решается методом прогонки (2.1.77) с граничными условиями (2.1.78). После ее решения компоненты собственного вектора $f_{i}^{{\vartheta ,0,s + 1}}$ и $f_{i}^{{\vartheta ,r,s + 1}}$ находятся с помощью уравнений (1.2.30).

Третья подсистема состоит из усредненных по переменным $z$ и $\vartheta $ уравнений (2.1.26):

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

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

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

2.3. Определение итерационных параметров циклического МР

В двумерном случае [17, разд. 3.4] после оценки границ собственных значений операторов $\hat {R}$ и $\hat {Z}$ для циклического ADI метода с длиной цикла $J$ рассчитываются величины, определяемые соотношениями:

где ${{\kappa }_{j}} = {{\kappa }_{j}}(\eta ,J)$ – известная функция [12] (в случае $J = 1$ $\kappa = \sqrt \eta $). Величины ${{\omega }_{{r,j}}}$ и ${{\omega }_{{z,j}}}$ в уравнении (2.3.2) определяют искомые оптимальные значения параметров ADI алгоритма.

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

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

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

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

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

Оптимальный выбор ${{\varepsilon }_{{P1}}}$ зависит от степени монотонности используемой разностной схемы и длины цикла $J$. Численный эксперимент показал, что при $J \geqslant 8$ для LD и LB-схем в 3D расчете достаточно задать ${{\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 \cdot {{2}^{{n - 1}}}$ раз, где $n = 1,2, \ldots $ – номер корректировки границ спектра. После корректировки верхней границы спектра итерации МР продолжаются (с использованием ранее полученного приближения для поправок) до сходимости, если не возникает необходимости в дополнительной корректировке. Корректировка верхней границы спектра производится до сходимости итераций МР, но не более 10 раз, после чего ускорение в группе отключается и используются чистые кинетические итерации.

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

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

Для предотвращения возможных переполнений целесообразно также осуществлять проверку значения прогоночного коэффициента ${{\xi }^{0}}$ в уравнениях (2.1.19), (2.1.39), (2.1.64), (2.1.77) и (2.1.85). Если абсолютное значение этого коэффициента в некоторой ячейке существенно больше $1$ (в 3D коде обычно используется критерий $\left| {{{\xi }^{0}}} \right| > 2.5$), то ускоряющая коррекция решения на текущей внутренней итерации не производится.

В $r,\;\vartheta ,\;z$ геометрии на границах расчетной области по переменной $\vartheta $ обычно используется либо условие зеркального отражения, либо периодическое граничное условие. Для обеспечения стабильности алгоритма ускорения в случае LB-схемы вес ${{P}_{\vartheta }}$ по переменной $\vartheta $ выбирается равным $P_{\vartheta }^{0} = 0.99$.

2.4. Численный пример использования MP в 3D геометрии

В качестве примера приведем некоторые результаты расчетов пространственного распределения потока нейтронов с энергией E > 0.5 МэВ в РУ ВВЭР-1000/320. Продольное и поперечное сечения расчетной области приведены на фиг. 1 и 2. Для задания исходной геометрии задачи методами комбинаторной геометрии использовался геометрический модуль программы MCU [20]. Комбинаторное задание геометрии затем конвертировалось на сетку задачи c использованием метода трэйсинга (tracing) программой ConDat [21] с поддержанием локального баланса масс в пределах каждой пространственной ячейки за счет введения дополнительных смесей материалов в рамках volume fraction (VF) метода.

Расчеты выполнены для 22-х верхних нейтронных групп проблемно-ориентированной константной системы BGL1000_B7 с использованием адаптивной WDD (AWDD) схемы и LD-схемы с коррекцией отрицательных потоков путем использования, аналогично [1], адаптивной WDD-схемы по угловой переменной $\varphi $, которую обозначим как WLD-схема. В качестве критерия сходимости итераций циклического МР с длиной цикла $J = 16$ использовался поблочный критерий (2.3.4) с размером ребра блока $N = 3$ и ${{\varepsilon }_{{P1}}} = 0.4$. При решении задачи для энергетических диапазонов $E > 3$ МэВ и $0.5 < E < 3.0$ МэВ использовались квадратуры Карлсона $E{{S}_{N}}$ [22] и Лебедева ${{L}_{{N,M}}}$ [23] на сфере. Квадратура Лебедева ${{L}_{{N,M}}}$ задается посредством двух параметров: N – порядок квадратуры (число параллелей на полусфере), N = 2, 4, 6, 8; M – порядок аксиальной симметрии квадратуры: M = 1, 2, … . В качестве источника использовался потвэльный источник для 2-й компании 2-го блока Ростовской АЭС.

В табл. 1 приведены пространственные сетки, угловые приближения и расчетные времена для решения задачи на ПК на базе Intel Core i7-6950X, использующем 128 Gb оперативной памяти и 500 Gb быстрой внешней памяти на SSD устройстве.

Для WLD-схемы число итераций МР для данной задачи практически постоянно. Оно равно 7 и не зависит от номера итерации и группы. Для AWDD-схемы возникает необходимость в корректировке верхней границы спектра, а число итераций МР колеблется в интервале 11–31 в зависимости от номера внутренней итерации и группы. При использовании WLD-схемы возможно также дальнейшее уменьшение пространственной сетки задачи по радиальной и аксиальной переменным без потери в точности рассчитываемого функционала.

На фиг. 3 приведены азимутальные распределения потока нейтронов с энергией E > 0.5 МэВ при r = 227.50 см (внешний радиус корпуса РУ) и z = 280.06 см (центр АЗ).

Анализ пространственных распределений в разных характерных точках задачи показывает, что оптимальный с точки зрения минимизации расчетного времени при сохранении необходимой точности расчета выбор угловой сетки задачи достигается при использовании комбинации квадратур ${{L}_{{12,2}}}$ при E > 3.0 МэВ и ${{L}_{{10,2}}}$ при $0.5 < E < 3.0$ МэВ.

2.5. Обсуждение результатов

Как показал численный эксперимент, МР в сочетании с WLD-схемой обычно сходится за 6–7 итераций и это число итераций практически не зависит от номера внутренней итерации и номера группы. Эта существенно лучший результат по сравнению с использованием МР в сочетании с адаптивной WDD (AWDD) схемой. Аналогичная ситуация наблюдается и при решении других задач. Вместе с тем использование WLD-схемы предполагает наличие больших объемов оперативной и быстрой внешней памяти для хранения рабочих массивов задачи, особенно при решении задач с большим числом групп.

Предложенный алгоритм решения ${{P}_{1}}$-системы для ускоряющих поправок используется также и для решения аналогичной ${{P}_{1}}$-системы, возникающей на каждой внешней итерации в алгоритме ускорения внешних итераций по области термализации нейтронов, а также при решении подкритической задачи [7], [24].