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

ВВЕДЕНИЕ

Конструирование метода интегрирования по времени уравнений Навье–Стокса обычно осуществляется методом прямых: при известной аппроксимации пространственных производных задача численного интегрирования сводится к решению системы обыкновенных дифференциальных уравнений и к такой задаче применяется какой-либо хорошо известный численный метод, разработанный для интегрирования обыкновенных дифференциальных уравнений.

Ограничение на шаг по времени, связанное с применением сильно неравномерных разностных сеток или, в общем случае, жесткость полученной системы требуют применения неявных методов. Неявный линеаризованный метод Эйлера блестяще зарекомендовал себя при решении стационарных задач методом установления, однако, имеет всего лишь 1-й порядок точности по времени. Для применения неявных методов Рунге–Кутты более высоких порядков необходимо на каждом временном шаге решать систему нелинейных уравнений. Размерность такой системы может быть большой. Для прямого численного моделирования турбулентности или достаточно сложной инженерной задачи в рамках осредненных по Рейнольдсу уравнений размерность системы может достигать несколько десятков миллиардов. Получить решение такой системы (при высоком порядке аппроксимации пространственных производных) с большой точностью итерационным (типа Ньютона) методом сложно. Как показывает опыт, даже в сравнительно простой ситуации – обтекание пластины, достаточно высокое число Рейнольдса, метод трапеций – для уменьшения начальной невязки на 2, 3 порядка требуется 3–5 итераций на каждом временном слое. Уменьшить невязку до уровня значений ${{10}^{{ - 7}}}{\kern 1pt} - {\kern 1pt} {{10}^{{ - 8}}}$ не удается ни за какое, сколько-нибудь разумное, число итераций. Трудность усугубляется тем, что a priori неизвестно, с какой точностью нужно решать систему нелинейных уравнений. Таким образом, интегрирование по времени неявным методом Рунге–Кутты может превратиться в очень трудоемкую процедуру со слабо контролируемой точностью вычислений.

Преодолеть указанную трудность позволяют методы типа Розенброка [1]. В этих методах на каждой стадии решается (без итераций) система линейных уравнений. Поэтому такие методы так же называют линейно неявными методами Рунге–Кутты.

Основная цель настоящей работы – разработка линейно неявного метода 3-го порядка точности по времени для уравнений Навье–Стокса. Предлагаемый подход является по сути вариантом W-метода [2]: в отличие от классического метода Розенброка в W-методе не требуется (точного) вычисления матрицы Якоби $A$. В настоящей работе предлагается выбирать матрицу $A$ так, чтобы соответствующий неявный оператор допускал факторизацию по координантным направлениям. Метод имеет 4 стадии – на каждой стадии решение линейной системы находится скалярными трехточечными прогонками.

Ранее, применительно к трехмерным уравнениям Навье–Стокса для сжимаемого газа, методы типа Розенброка повышенного порядка точности с управляемой длиной шага по времени рассматривались в ряде работ [3]–[6]. В [3]–[5] для аппроксимации пространственных производных использовался разрывный метод Галеркина (DG); в [6] – конечно-разностные (центральные разности с искусственной диссипацией) аппроксимации. Для решения системы линейных уравнений применялся метод обобщенных минимальных невязок (GMRES). На ряде тестовых задач оценивалась эффективность методов типа Розенброка и диагонально неявных методов Рунге–Кутты. Предлагаемый в настоящей работе метод по конструкции наиболее близок к (рассмотренному в [3], [6]) четырехстадийному W-методу ROS34PW2 3-го порядка, разработанному в [7] для интегрирования жестких систем уравнений.

Для оценки точности предложенного метода в качестве тестового примера рассматривается задача с начальными данными (задача Коши) на трехмерном торе. В уравнения Навье–Стокса для сжимаемого газа добавляются источниковые члены таким образом, чтобы предельное при $t \to \infty $ решение являлось ABC-потоком. Для аппроксимации по пространству применяется разностная схема типа Годунова 3-го (на гладких решениях) порядка точности. Численное интегрирование по времени осуществлялось тремя различными способами: линейно неявным методом Эйлера и линейно неявными методами 2-го и 3-го порядков точности. Рассчитанные распределения сопоставляются с решением, задаваемым ABC-потоком.

1. ИНТЕГРИРОВАНИЕ ПО ВРЕМЕНИ

Для системы обыкновенных дифференциальных уравнений $y{\kern 1pt} ' = f(y),\;y({{t}_{0}}) = {{y}_{0}}$ модифицированный W-метод записывается в виде

где $h$ – шаг по времени, $s$ – число стадий, ${{\alpha }_{{ij}}},\;{{\gamma }_{{ij}}},\;{{b}_{i}}$ – коэффициенты метода. В отличие от классического W-метода [2], в (1.1) матрица $A$ предполагается зависящей произвольным (гладким) образом от ${{g}_{i}}$ и ${{\gamma }_{{ii}}}h$. Зависимость от ${{g}_{i}}$ выбирается из соображений устойчивости и для эффективного обращения матрицы неявного оператора (см. разд. 2). Зависимость от ${{\gamma }_{{ii}}}h$ (для автономных систем) возникает из-за следующего обстоятельства. Перепишем первое уравнение в (1.1) в следующем виде: где $I$ – единичная матрица. В таком виде удается избежать умножения на матрицу $A$ и задача сводится к обращению неявного оператора в левой части (1.2). Для эффективного нахождения неизвестных ${{\hat {k}}_{i}}$ выберем матрицу $A$ так, чтобы неявный оператор допускал факторизацию по координантным направлениям: где матрицы ${{A}_{\xi }},{{A}_{\eta }},{{A}_{\zeta }}$ суть функции ${{g}_{i}}$. Из (1.3) видно, что матрица $A$ в общем случае зависит от ${{\gamma }_{{ii}}}h$. Применительно к уравнениям Навье–Стокса явный вид матриц ${{A}_{\xi }},{{A}_{\eta }},{{A}_{\zeta }}$ будет приведен в разд. 2.

Для метода (1.1), (1.3) условия 3-го порядка примут следующий вид:

где суммирование происходит по всем индексам от 1 до $s$ и ${{\alpha }_{{ij}}} = 0,\;j \geqslant i,\;{{\gamma }_{{ij}}} = 0,\;j > i$ . Первые три условия $(A1){\kern 1pt} - {\kern 1pt} (A3)$ в (1.4) соответствуют порядку 2. Условия $(A1){\kern 1pt} - {\kern 1pt} (A3),\;(B1){\kern 1pt} - {\kern 1pt} (B5)$ совпадают с условиями для W-метода, приведенными в [2] (см. также [8, c. 135]). Условие $(C1)$ возникает из-за зависимости $A$ от ${{g}_{i}}$; условие $(C2)$ – из-за зависимости $A$ от ${{\gamma }_{{ii}}}h$ в силу (1.3). Потребуем, чтобы диагональные члены ${{\gamma }_{{ii}}}$ были равными. Это требование упростит дальнейшее рассмотрение. Тогда $(C2)$ будет выполнено в силу условия $(A3)$ в (1.4). Поэтому далее будем считать что условие $(C2)$ заменено на $s - 1$ уравнений ${{\gamma }_{{11}}} = {{\gamma }_{{22}}} = \ldots = {{\gamma }_{{ss}}}$.

Трехстадийный, L-устойчивый (для скалярного линейного уравнения) метод 3-го порядка, удовлетворяющий (1.4), невозможен. Это будет показано в конце раздела. Поэтому далее будет рассматриваться случай $s = 4$.

При $s = 4$ число неизвестных (тех, которые могут быть отличны от нуля) ${{b}_{i}} - 4,\;{{\alpha }_{{ij}}} - 6,\;{{\gamma }_{{ij}}} - 10$. Всего 20 неизвестных и 12 уравнений. Исследование решений такой системы во всей полноте довольно трудоемко. Ниже будет указан (несложный) алгоритм, позволяющий легко получать частные, L-устойчивые решения системы (1.4).

Условия $(A1),\;(A2),\;(B1),\;(B2)$ в (1.4) совпадают с условиями на коэффициенты для явного метода Рунге–Кутты 3-го порядка. Поэтому величины ${{b}_{i}},\;{{\alpha }_{{ij}}}$ можно считать известными. При известных ${{b}_{i}},\;{{\alpha }_{{ij}}}$ задача сводится к нахождению ${{\gamma }_{{ij}}}$. Удобно ввести следующие величины:

Уравнения для ${{\gamma }_{i}}$ и ${{\sigma }_{k}}$ в силу (1.4), (1.5) примут следующий вид:

По определению (1.5) ${{\sigma }_{4}} = {{b}_{4}}{{\gamma }_{{44}}} = {{b}_{4}}{{\gamma }_{1}}$, поэтому можно считать, что неизвестных ${{\sigma }_{k}}$ только 3.

Заметим, что при известных ${{\gamma }_{i}},{{\sigma }_{k}}$( удовлетворяющих (1.6)) легко находится общее (однопараметрическое) решение для ${{\gamma }_{{ij}}}$. Приведем ответ:

где ${{\gamma }_{{32}}}$ – параметр. Метод предполагается 4-стадийным, поэтому ${{b}_{4}} \ne 0$.

Перейдем к ограничениям на коэффициенты, которые накладывает условие L-устойчивости. В общем случае исследование устойчивости методов типа Розенброка очень сложно. Рассмотрим скалярный линейный случай $f = \lambda y$. Напомним определения: метод называется A-устойчивым, если $\left| {R(z)} \right| \leqslant 1$ при $z \in {{\mathbb{C}}_{ - }}$, где $z = \lambda h$, ${{\mathbb{C}}_{ - }}$ – левая полуплоскость, $R(z)$ – функция устойчивости, ${{y}_{1}} = R(z){{y}_{0}}$; метод называется L-устойчивым, если он А-устойчив и $R(\infty ) = 0$. При $A = \lambda $ метод (1.1) сводится к однократно диагонально неявному методу Рунге–Кутты (SDIRK) с матрицей ${{a}_{{ij}}} = {{\alpha }_{{ij}}} + {{\gamma }_{{ij}}}$. Такие методы хорошо изучены. Условия А-устойчивости для таких методов при $R(\infty ) = 0$ приведены в [8, c. 117]. В нашем случае ${{\gamma }_{1}}$ должна удовлетворять неравенству

Далее, для $R(\infty )$ справедлива формула (см. [8, c. 59]): Условие $R(\infty ) = 0$, вычисленное по формуле (1.9), в нашем случае принимает вид В (1.10) коэффициенты ${{c}_{i}}$ явно зависят только от ${{b}_{i}},\;{{\alpha }_{{ij}}},\;{{\gamma }_{i}},\;{{\sigma }_{j}}$, но не от ${{\gamma }_{{ij}}}$.

Приведем алгоритм нахождения частных, L-устойчивых решений (1.4).

АЛГОРИТМ

Шаг 1. Находятся/задаются ${{b}_{j}},\;{{\alpha }_{{ij}}}$, удовлетворяющие в (1.4) уравнениям $(A1),\;(A2),\;(B1),\;(B2).$

Шаг 2. Из первых трех линейных уравнений (1.6) для ${{\gamma }_{i}}$ находятся ${{\gamma }_{2}},\;{{\gamma }_{3}},\;{{\gamma }_{4}}$ как линейные функции ${{\gamma }_{1}}$. При известных ${{\gamma }_{i}}$ из последних трех линейных уравнений (1.6) для ${{\sigma }_{i}}({{\sigma }_{4}} = {{b}_{4}}{{\gamma }_{1}})$ находятся ${{\sigma }_{1}},\;{{\sigma }_{2}},\;{{\sigma }_{3}}$ как рациональные функции ${{\gamma }_{1}}$.

Шаг 3. Вычисляются коэффициенты ${{c}_{3}}$ и ${{c}_{4}}$ в (1.10) и находится ${{\gamma }_{{32}}}$ как (рациональная) функция ${{\gamma }_{1}}$.

Шаг 4. По (1.7) вычисляются ${{\gamma }_{{ij}}}$. Задается величина ${{\gamma }_{1}}$, удовлетворяющая (1.8), и определяются числовые значения коэффициентов.

Пример 1. Пусть выполнено ${{\beta }_{i}} = {{b}_{i}}(1 - {{\alpha }_{i}}),\;i = 2,3,\;{{\alpha }_{4}} = 0$. Tогда из $(A1),\;(A2),\;(B1),\;(B2)$ в (1.4) находим 4-х параметрическое (параметры ${{\alpha }_{2}},\;{{\alpha }_{3}},\;{{\alpha }_{{32}}},\;{{b}_{4}}$) решение:

Далее детально рассмотрим случай ${{\alpha }_{2}} = 1{\text{/}}3,\;{{\alpha }_{3}} = 1{\text{/}}2,\;{{b}_{4}} = 1{\text{/}}2$. Применяя шаг 1 алгоритма, по формулам (1.11) получаем таблицу Бутчера для ${{b}_{i}}$ и ${{\alpha }_{{ij}}}$ (${{\alpha }_{{32}}}$ – параметр):

На шаге 2 из (1.6) получаем (${{\beta }^{{\text{т}}}} = (1{\text{/}}2, - 1,1,0)$):

На шаге 3, используя (1.12), (1.13) для коэффициентов в (1.10), получаем ${{c}_{3}} = 1{\text{/}}6 - {{\gamma }_{1}} + \gamma _{1}^{2},$ ${{c}_{4}} = (2{{\gamma }_{1}} - 1){\text{/}}3$. Условие (1.10) перепишется в виде На шаге 4 из (1.7), (1.12) для ${{\gamma }_{{ij}}}$ получаем (1.12), (1.14), (1.15) задают двухпараметрическое (${{\alpha }_{{32}}},\;{{\gamma }_{1}}$ – параметры) семейство решений.

Пример 2. Условиям $(A1),\;(A2),\;(B1),\;(B2)$ в (1.4) будут удовлетворять коэффициенты любого явного метода Рунге–Кутты 4-го порядка точности. На шаге 1 возьмем “правило 3/8”:

На шаге 2 система в (1.6) для определения ${{\gamma }_{2}},\;{{\gamma }_{3}},\;{{\gamma }_{4}}$ в данном случае оказывается вырожденной, но разрешимой. Для простоты выкладок положим ${{\gamma }_{4}} = {{\gamma }_{1}}$. Тогда из (1.6) получаем (${{\beta }^{{\text{т}}}} = (1{\text{/}}8,1{\text{/}}4,1{\text{/}}8,0)$):

Подсчет коэффициентов в (1.10) на шаге 3 по (1.16), (1.17) показывает, что ${{c}_{i}},i \leqslant 3$, совпадают с ${{c}_{i}}$ примера 1, а для ${{c}_{4}}$ и ${{\gamma }_{{32}}}$ имеем формулы На шаге 4 положим ${{\gamma }_{1}} = 1{\text{/}}2$, тогда из (1.18) ${{\gamma }_{{32}}} = - 3{\text{/}}4$ и ${{\gamma }_{{ij}}}$, вычисленные по (1.7), имеют вид

Покажем, что не существует трехстадийного L-устойчивого метода 3-го порядка, т.е. метода, для которого выполнено (1.4) (без условия $(C2)$ и предположения о равенстве ${{\gamma }_{{ii}}}$) при $s = 3$. Будем считать, что ${{\gamma }_{{11}}}{{\gamma }_{{22}}}{{\gamma }_{{33}}} \ne 0$. В противном случае функция устойчивости $R(z)$ будет иметь вид $R(z) = P(z){\text{/}}Q(z),$ $\deg Q(z) = \deg (1 - {{\gamma }_{{11}}}z)(1 - {{\gamma }_{{22}}}z)(1 - {{\gamma }_{{33}}}z) = 2$. (Случай $\deg Q \leqslant 1$, очевидно, невозможен из-за условия порядка.) Из условия L-устойчивости ($R(\infty ) = 0$) следует, что $\deg P = 1$. Из условия порядка ($p = 3$) вытекает, что $R(z)$ – (2.1) аппроксимация Паде для ${{e}^{z}}$. Такая аппроксимация единственна; из явного вида для $Q(z)$ следует, что $Q(z)$ имеет только комплексные корни. Получили противоречие (${{\gamma }_{{ii}}}$ – действительны). Из условий ${{\gamma }_{{ii}}} \ne 0$ вытекает, что введенные в (1.5) векторы $\gamma $ и $\sigma $ отличны от нуля. Далее, общее решение уравнений $(A1),\;(A2),\;(B1),\;(B2)$ в (1.4) при $s = 3$ дается формулами (1.11) при ${{b}_{4}} = 0,\;{{\alpha }_{{32}}} = 1{\text{/}}(6{{b}_{3}}{{\alpha }_{2}})$. Три первых однородных уравнения в (1.6) должны иметь нетривиальное рещение $\gamma $. Приравняв соответствующий детерминант нулю и воспользовавшись формулами (1.11), получим ${{\alpha }_{3}} = 1$. Разрешим эту систему относительно ${{\gamma }_{1}}$, тогда ${{\gamma }_{2}}{\text{/}}{{\gamma }_{1}}$ и ${{\gamma }_{3}}{\text{/}}{{\gamma }_{1}}$ будут зависеть только от ${{\alpha }_{2}}$. Поэтому коэффициенты однородной системы (1.6) для $\sigma $ можно считать зависящими только от ${{\alpha }_{2}}$. Приравняв нулю соответствующий определитель, получим ${{\alpha }_{2}} = 1$. Получили противоречие: ${{\alpha }_{2}} = {{\alpha }_{3}}$. Случаи, неохваченные формулами (1.11) при ${{\alpha }_{2}} = {{\alpha }_{3}}$ или ${{\alpha }_{3}} = 0$, рассматриваются аналогично. Случай ${{\alpha }_{2}} = 0$ невозможен.

2. ФАКТОРИЗАЦИЯ НЕЯВНОГО ОПЕРАТОРА И АППРОКСИМАЦИЯ ПРОСТРАНСТВЕННЫХ ПРОИЗВОДНЫХ

Факторизация (1.3) осуществляется так же, как для линейно неявного метода Эйлера. Такой подход хорошо известен, см., например, [9]. Приведем основные шаги его построения для уравнений Навье–Стокса, записанных в криволинейной системе координат $\xi ,\eta ,\zeta (\xi = \xi (x,y,z),$ $\eta = \eta (x,y,z),$ $\zeta = \zeta (x,y,z),$ $x,\;y,\;z$ – декартовые координаты):

где индексом $v$ обозначены диффузионные члены уравнений Навье–Стокса.

Достаточно описать вид первого множителя в правой части (1.3); остальные множители строятся аналогично с очевидной заменой. Пусть

где $\Lambda ({\text{|}}\Lambda {\text{|}})$ – диагональная матрица, составленная из (модуля) собственных значений матрицы Якоби $J,\;{{S}_{\xi }}$ – матрица (размером $5 \times 5$) , приводящая $J$ к диагональному виду. Пусть тройка целых чисел n, m, l – номер центра ячейки разностной сетки; шаги сетки считаем равномерными, равными 1. Далее, считаем индексы m, l фиксированными, поэтому их опускаем. Тогда действие первого множителя в (1.3) на cеточную функцию ${{f}_{n}}$ задается в виде где $\nu $ – величина, пропорциональная коэффициенту кинематической вязкости. Обращение матрицы в (2.1) на каждом этапе сводится к умножениям матриц ${{S}_{\xi }},\;S_{\xi }^{{ - 1}}$ на вектор и решению трехдиагональной системы уравнений. Такое решение легко находится трехточечной скалярной прогонкой.

Применение направленных (в зависимости от знака собственных значений) разностей в (2.1) позволяет получить свойство диагонального преобладания (что способствует повышению устойчивости метода прогонки), однако, зависимость ${{A}_{\xi }}$ от ${{g}_{i}}$ становится, вообще говоря, негладкой. Гладкость нарушается в точках, в которых собственное значение обращается в нуль. По-видимому, можно надеяться, что в большинстве случаев таких точек относительно немного и это обстоятельство не оказывает решающего влияния на точность расчета.

Способ аппроксимации конвективных членов проще всего описать для случая скалярного уравнения; аппроксимация диффузионных (вязких) членов осуществлялась центральными разностями. Введем по оси $x$ равномерную сетку с шагом $h$ и проинтегрируем уравнение ${{u}_{t}} + f{{(u)}_{x}} = 0$ по ячейке с номером с $n$. Получим

Для линейного уравнения $f = au$ схема типа Годунова имеет вид где ${{U}_{{ \pm ,n + 1/2}}}$ – предельные справа (слева) значения $u(x,t)$ на границе ячейки $n + 1{\text{/}}2$. Величины ${{U}_{{ \pm ,n \mp 1/2}}}$ вычисляются по рекуррентным формулам: где $C_{n}^{l}$ – биномиальный коэффициент, число сочетаний из $n$ по $l;{{\Delta }^{l}}{{u}_{n}}$ – конечная разность $l$-го порядка: ${{\Delta }^{l}}{{u}_{n}} = {{\Delta }^{{l - 1}}}{{u}_{{n + 1}}} - {{\Delta }^{{l - 1}}}{{u}_{n}}$; $\phi _{l}^{k}(x,y)$ – некоторые однородные степени нуль функции двух переменных. Число $k$ определяет точность схемы – аппроксимация производной по $x$ будет точной на многочленах степени $k - 1$; случай $k = 2$ соответствует TVD-схеме.

Применительно к уравнениям Навье–Стокса/Эйлера схема типа Годунова строится следующим, хорошо известным способом. Формулы (2.2) применяются для “характеристических переменных” $W = S_{\xi }^{{ - 1}}U$. Далее, с помощью процедуры решения задачи о произвольном разрыве [10] вычисляется поток на границе ячейки $\xi = n + 1{\text{/}}2$. Для направлений $\eta ,\;\zeta $ построение аналогично.

3. ПОСТАНОВКА ЗАДАЧИ

Задача решалась в рамках трехмерных уравнений Навье–Стокса для сжимаемого газа с источниковыми членами $F,\;{{F}_{E}}$:

где $\rho $ – плотность, $v$ – скорость, $p$ – давление, $T$ – температура, $\delta $ – единичный тензор, $\gamma $ – показатель адиабаты, ${\text{Re,}}\;{\text{Pr}}$ – числа Рейнольдса и Прандтля. Система (3.1) записана в безразмерной форме: линейные размеры отнесены к некоторой характерной величине ${{L}_{\infty }}$, $\rho $ отнесено к ${{\rho }_{\infty }}$, $v$ – к ${{V}_{\infty }}$ , $p$ – к ${{\rho }_{\infty }}V_{\infty }^{2}$, $T$ – к $V_{\infty }^{2}{\text{/}}{{c}_{p}},({{c}_{p}}$ – удельная теплоемкость при постоянном объеме), $\mu $ – к ${{\mu }_{\infty }}$.

Источниковые члены $F,{{F}_{E}}$ задавались так, чтобы стационарное решение системы (3.1) являлось ABC-потоком:

Явный вид $F,{{F}_{E}}$ легко получить, подставляя (3.2) в (3.1):

Для системы (3.1), (3.3) ставилась задача Коши на трехмерном торе ${{\mathbb{T}}^{3}}$: где М – число Маха. Предполагая, что при $t \to \infty $ решение задачи (3.1) стремится к стационарному решению (3.2), легко определить константу $H$ в (3.2). Проинтегрируем уравнение энергии по области ${{\mathbb{T}}^{3}} \times [0,t]$. С учетом (3.4) получим (интеграл от источникового члена ${{F}_{E}}$ равен нулю) Угловые скобки в (3.5) – осреднение по тору ${{\mathbb{T}}^{3}}$. Устремляя $t \to \infty $, из (3.2), (3.5) находим константу Бернулли:

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Численное решение задачи (3.1), (3.3), (3.4) находилось при значениях параметров $\gamma = 1.4,$ $\mu = 1,$ ${\text{Pr}} = 0.72,$ ${\text{Re = 10}}0,$ ${\text{M}} = 0.5,$ $A = 1,$ $B = 0.5,$ $C = 0,$ $\lambda = 2\pi $.

Шаги разностной сетки $h$ выбирались равномерными по всем направлениям $h = 1{\text{/}}N$, где $N$ – число узлов в каком-либо направлении. Число Куранта $K$ вычислялось по формуле ${\text{K}} = 3\Delta t({\text{M}} + 1){\text{/}}({\text{M}}h)$, где $\Delta t$ – величина шага по времени.

Аппроксимация пространственных производных осуществлялась с 3-м порядком точности по формулам (2.2) с $k = 4,$ $\phi _{l}^{n} = \phi (x,y) = \left| x \right|{\text{/}}(\left| x \right| + \left| y \right|)$; вязкие члены аппроксимировались с 4-м порядком центральными разностями.

Интегрирование по времени с 3-м порядком осуществлялось методом (1.16), (1.19). Также проводились расчеты методом (1.12), (1.15) с ${{\alpha }_{{32}}} = 0,$ ${{\gamma }_{{32}}} = - 1{\text{/}}9,$ ${{\gamma }_{{11}}} = 1{\text{/}}3$. Распределения, полученные этим методом, оказались близки к распределениям, рассчитанным по (1.16), (1.19) и поэтому далее не приводятся. Для сравнения выполнялись расчеты линейно неявным методом Эйлера (${{\alpha }_{{11}}} = 0,\;{{b}_{1}} = {{\gamma }_{1}} = 1$) и методом 2-го порядка [11]:

На границах вычислительной области $(x = 0,1;\;y = 0,1;\;z = 0,1)$ выставлялись условия периодичности. Обращение матрицы неявного оператора осуществлялось скалярной циклической, трехточечной прогонкой (см., например, [12, c. 86]).

Временные затраты расчета каким-либо методом с большой точностью пропорциональны числу его стадий: время выполнения 1-го шага по времени для методов 2-го и 3-го порядков в 2 и, соответственно, 4 раза больше по сравнению с линейно неявным методом Эйлера.

Точное нестационарное решение задачи (3.1), (3.3), (3.4) неизвестно, поэтому точность метода интегрирования по времени оценивалась сравнением со стационарным (пределом при $t \to \infty $) решением (3.2). Заметим, что численное “стационарное” решение данной задачи зависит от метода интегрирования по времени. Действительно, если сделать 1 шаг по времени двумя различными методами и, приняв полученные (и восполненные соответствующим образом) распределения в качестве начальных данных, получить точное решение, то, как видно из (3.5), эти (точные) решения будут, вообще говоря, различаться. Этим, в частности, задача на торе отличается, например, от задачи расчета стационарного обтекания профиля, – для такой задачи стационарное решение определяется аппроксимацией пространственных производных/граничных условий и не зависит от метода интегрирования по времени и начального приближения.

Сопоставление численного решения с решением (3.2) осуществлялось следующим образом.

Рассчитанные распределения в сечении $x = 0.5,\;y = 0.5$ записывались на каждом временном шаге и затем осреднялись на промежутке $80 \leqslant t \leqslant 100$. Как показывают расчеты, при $t \geqslant 50$ численное решение можно считать стационарным с довольно большой точностью: величина невязки решения в норме ${{l}_{\infty }}$ равнялась $ \approx {\kern 1pt} {{10}^{{ - 7}}}{\kern 1pt} - {\kern 1pt} {{10}^{{ - 9}}}$.

Хорошо известно (см., например, [13]), что течения, задаваемые ABC-потоком, могут оказаться крайне неустойчивыми. При выбранных в настоящей работе параметрах численное решение, полученное методом 3-го порядка при $N = 33,K = 1$, оставалось устойчивым до $t \leqslant 500$. Что произойдет с численным решением при $t$, равным, например, ${{10}^{6}}$ – неизвестно.

Из (3.2), (3.6) следует, что при выбранных параметрах $A,\;B,\;C$ при $x = 0.5,$ $y = 0.5$ давление $p$ не зависит от $z$ и равно ${{p}_{ * }} = 1{\text{/}}(\gamma {{M}^{2}}) - 0.5(\gamma - 1)({{A}^{2}} + {{B}^{2}} + {{C}^{2}}) = 2.60714285$.

Рассчитанные распределения $p(z)$ при $K = 4$ и $N = 9,15,33,63$ приведены на фиг. 1а, б, в, г соответственно. На фиг. 1 и далее сплошная линия с маркером $ \circ $ соответствует расчету методом Эйлера, линия с $\square $ – расчету методом 2-го порядка, линия с $ + $ – расчету методом 3-го порядка. Горизонтальная штриховая линия – $p = {{p}_{ * }}$. На фиг. 1 по оси ординат выбран различный масштаб. Распределение ошибки ${{\varepsilon }_{p}}$ для давления в норме ${{l}_{\infty }}$ в зависимости от величины шага $h$ при $K = 4$ приведено на фиг. 2 (обе шкалы логарифмические); распределения в нормах ${{l}_{1}},\;{{l}_{2}}$ аналогичны и поэтому не приводятся. Штриховые прямые на фиг. 2 – зависимости, пропорциональные $h,\;{{h}^{2}}$, ${{h}^{3}}$ и ${{h}^{4}}$. Аналогичное распределение для ошибки скорости ${{\varepsilon }_{v}}$ приведено на фиг. 3, на этой фигуре штриховые прямые – зависимости, пропорциональные ${{h}^{3}}$ и ${{h}^{4}}$. Как и следовало ожидать, при фиксированном $h$ точность расчета возрастает с увеличением порядка аппроксимации по времени. При этом распределения давления, полученные различными методами, с очень большой точностью отличаются сдвигом на постоянную по оси ординат (фиг. 1). Как видно из фиг. 2, для метода Эйлера ошибка ${{\varepsilon }_{p}}$ падает чуть быстрее, чем $O(h)$, для метода 2-го порядка точность метода близка к теоретической, для метода 3-го порядка убывание ошибки (примерно) соответствует порядку аппроксимации пространственных производных. Ошибка ${{\varepsilon }_{v}}$ для компонент скорости слабо зависит от метода интегрирования по времени и порядок убывания ошибки также соответствует 3–4 порядку аппроксимации пространственных производных (фиг. 3). При этом отличие в компонентах скорости, полученных различными методами, примерно на 2 порядка меньше отличия в давлении. Распределения температуры вполне аналогичны распределениям давления и поэтому не приводятся. Величина сдвига в распределениях давления с точностью $ \approx 10\% $ равна разности констант $H$: для $N = 9$ рассчитанные значения констант равны 3.22570980, 3.23797792, 3.24142873, для $N = 15$ – 3.22619314, 3.23271241, 3.23427132, для $N = 33$ – 3.22898588, 3.23158921, 3.23220390, для $N = 63$ – 3.23052568, 3.23189454, 3.23214559 для методов 1‑го–3-го порядков соответственно. Значение, рассчитанное по (3.6), $H = 3.23214285$.

Результаты расчетов для $N = 33$ и $K = 0.5,1,2,4$ приведены на фиг. 4, 5. На фиг. 4 (обе шкалы логарифмические) приведено распределение ошибки ${{\varepsilon }_{p}}$ для давления в зависимости от величины шага по времени $\Delta t$. Штриховая прямая на фиг. 4 соответствует зависимости, пропорциональной $\Delta t$ . Как видно из графиков, для метода Эйлера уменьшение ошибки ${{\varepsilon }_{p}}$ близко к линейной зависимости. Для метода 3-го порядка ошибка ${{\varepsilon }_{p}}$ практически не меняется от величины шага по времени и, по-видимому, приблизительно равна ошибке аппроксимации по пространственным производным. Штриховая кривая с маркером $\diamondsuit $ на фиг. 4 – зависимость $5.13016 \times {{10}^{{ - 5}}} + 209.95\Delta {{t}^{2}}$, – в этом выражении первый член соответствует ошибке ${{\varepsilon }_{p}}$ для метода 3-го порядка, а коэффициент при $\Delta {{t}^{2}}$ определялся из условия равенства величине ошибки для метода 2-го порядка при $\Delta t = 1.3468 \times {{10}^{{ - 2}}}$. Аналогичное распределение для ошибки скорости ${{\varepsilon }_{v}}$ представлено на фиг. 5 (шкала по оси абсцисс – логарифмическая). Как видно из графиков, влияние величины шага по времени и порядка метода интегрирования по времени на компоненты скорости относительно невелико.

Таким образом, для данной задачи при $t \to \infty $ численные решения для компонент скорости, полученные различными методами, – близки и в основном определяются аппроксимацией пространственных производных. Распределения давления и температуры с точностью до аддитивной постоянной также хорошо совпадают. Порядок аппроксимации метода по времени главным образом влияет на рассчитанную величину константы Бернулли $H$.