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

1.1. Система уравнений

Рассмотрим безграничный объем ньютоновской несжимаемой жидкости. Ее изотермическое течение подчиняется уравнениям Навье–Стокса:

где $x$, $y$, $z$ – декартовы пространственные координаты, $t$ – время, ${\mathbf{v}} = ({{{v}}_{x}},{{{v}}_{y}},{{{v}}_{z}})$ – вектор скорости жидкости, $p$ – давление жидкости, ${\text{Re}} = {\text{const}} > 0$ – число Рейнольдса, символ ${{\partial }_{t}} \equiv \partial {\text{/}}\partial t$, символ $ \otimes $ – знак тензорного произведения. Система уравнений (1) уже записана в безразмерной форме, при этом конвективный член $({\mathbf{v}} \times \nabla ){\mathbf{v}}$ в левой части уравнений движения представлен с помощью уравнения неразрывности в дивергентном виде как ${\text{div}}({\mathbf{v}} \otimes {\mathbf{v}})$.

1.2. Расщепление по физическим процессам

Построим для уравнений Навье–Стокса (1) бикомпактную схему четвертого порядка аппроксимации по пространству. Воспользуемся методом расщепления Марчука–Стрэнга [25[, [26] по физическим процессам. Систему (1) можно разбить на три части – конвективную:

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

Ясно, что любая машинная реализация упомянутого алгоритма потребует конечности расчетной области и постановки каких-то условий на ее границах. Поскольку в настоящей работе проводятся расчеты задачи о распаде вихря Тейлора–Грина, имеющей периодическое решение, ограничимся для определенности областью $D = {{(0,2\pi )}^{3}}$ и периодическими условиями на ее границе $\partial D$.

1.3. Сетка и представление решения на ней

Введем в замкнутой области $\bar {D}$ сетку

Условимся нумеровать узлы по направлениям $x$, $y$, $z$ индексами $j$, $k$, $\ell $ соответственно. Введем также временные слои $0 = {{t}^{0}} < {{t}^{1}} < {{t}^{2}} < \; \ldots $, а через $\tau $ обозначим (быть может, переменный) шаг по времени. Определим несколько сеточных функций:

1) ${\mathbf{V}}_{{j,k,\ell }}^{n}$ – задана во всех узлах сетки $\Omega $ (индексы $j$, $k$, $\ell $ принимают как целые, так и полуцелые значения) и аппроксимирует в них скорость ${\mathbf{v}}$ в момент времени $t = {{t}^{n}}$, $n = 0,\;1, \ldots $;

2) ${\mathbf{v}}_{{j,k,\ell }}^{n}$, $({\mathbf{v}}_{x}^{'})_{{j,k,\ell }}^{n}$, $({\mathbf{v}}_{y}^{'})_{{j,k,\ell }}^{n}$, $({\mathbf{v}}_{z}^{'})_{{j,k,\ell }}^{n}$, $({\mathbf{v}}_{{xy}}^{{''}})_{{j,k,\ell }}^{n}$, $({\mathbf{v}}_{{yz}}^{{''}})_{{j,k,\ell }}^{n}$, $({\mathbf{v}}_{{xz}}^{{''}})_{{j,k,\ell }}^{n}$, $({\mathbf{v}}_{{xyz}}^{{'''}})_{{j,k,\ell }}^{n}$ – заданы только в целых узлах сетки $\Omega $ (индексы $j$, $k$, $\ell $ принимают только целые значения) и аппроксимируют в них при $t = {{t}^{n}}$ значения скорости ${\mathbf{v}}$ и семи ее производных по $x$, $y$, $z$, $xy$, $yz$, $xz$, $xyz$ соответственно. На совокупность этих сеточных функций можно смотреть как на одну сеточную функцию с большим числом компонент на один узел (24 компоненты – 3 компоненты скорости и по 7 производных для каждой из них). Отметим, что такая функция “хранит” примерно столько же чисел, сколько и функция ${\mathbf{V}}_{{j,k,\ell }}^{n}$, это количество порядка $24{{N}_{x}}{{N}_{y}}{{N}_{z}}$ при ${{N}_{x}},{{N}_{y}},{{N}_{z}} \gg 1$;

3) $p_{{j,k,\ell }}^{n}$, $(p_{x}^{'})_{{j,k,\ell }}^{n}$, $(p_{y}^{'})_{{j,k,\ell }}^{n}$, $(p_{z}^{'})_{{j,k,\ell }}^{n}$, $(p_{{xy}}^{{''}})_{{j,k,\ell }}^{n}$, $(p_{{yz}}^{{''}})_{{j,k,\ell }}^{n}$, $(p_{{xz}}^{{''}})_{{j,k,\ell }}^{n}$, $(p_{{xyz}}^{{'''}})_{{j,k,\ell }}^{n}$ – то же, что и в п. 2), только для давления $p$.

1.4. Аппроксимация конвективной части уравнений Навье–Стокса

Представим систему (2) в виде

где векторы ${{{\mathbf{F}}}^{d}}({\mathbf{v}}) = {{{v}}_{d}}{\mathbf{v}}$, $d = x,y,z$. Так как мы уже начали пользоваться методом расщепления для численного решения системы (1), разумно продолжить пользоваться этим методом и для дискретизации системы (2), чтобы добиться наибольшей экономичности результирующей численной схемы для полной системы (1). Применим к системе многомерных уравнений (5) локально-одномерное расщепление по пространству [25], [27], [28] (см. также [16]) и глобальное потоковое расщепление Лакса–Фридрихса (см. [12]), что приведет к ее разбиению на шесть одномерных систем вида: где расщепленные потоки ${{{\mathbf{F}}}^{{d, \pm }}}({\mathbf{v}})$ задаются следующими выражениями: Параметры $C_{2}^{d} > 0$ подбираются так, чтобы на некотором множестве ${\mathbf{v}}$, например, множестве значений сеточной функции ${\mathbf{V}}_{{j,k,\ell }}^{n}$, матрицы Якоби потоков (7) были соответственно положительно- и отрицательно-определенными. В расчетах настоящей работы параметры $C_{2}^{d}$ остаются фиксированными в течение каждого перехода со слоя ${{t}^{n}}$ на слой ${{t}^{{n + 1}}}$, но между этими переходами пересчитываются по формуле: где $\delta > 0$ – константа, отвечающая за “запас” положительной/отрицательной определенности матриц ${{{\mathbf{A}}}^{ \pm }}({\mathbf{v}})$.

Наконец, каждая из шести систем (6) решается при помощи одномерной бикомпактной схемы с аппроксимацией по времени методом трапеций (как в схеме Кранка–Николсона). В силу аналогичности этих систем друг другу мы остановимся подробно только на системе с потоком ${{{\mathbf{F}}}^{{x, + }}}({\mathbf{v}})$. Для нее указанная бикомпактная схема запишется следующим образом [12]:

Поясним обозначения: ${{\hat {V}}_{{j,k,\ell }}}$ и – это сеточные функции типа ${\mathbf{V}}_{{j,k,\ell }}^{n}$; первая означает искомое решение на неком промежуточном слое $t = \hat {t}$, а вторая – известное решение на промежуточном слое ; в данном случае . Далее, в описании общего алгоритма счета объясняется, что конкретно фигурирует в схеме (9) вместо ${\mathbf{\hat {V}}}$ и . Символами $A_{0}^{x}$, $\Lambda _{1}^{x}$, $\Lambda _{2}^{x}$ обозначены следующие одномерные сеточные операторы: где ${{U}_{{j,k,\ell }}}$ – произвольная сеточная функция, определенная во всех узлах $\Omega $. Потоки $\hat {F}_{{j,k,\ell }}^{{x, + }}$ и для всех $j$, $k$, $\ell $ (и целых, и полуцелых) определяются формулами: Схема (9) дополняется периодическим граничным условием: Искомое решение ${\mathbf{\hat {V}}}$ вычисляется из нелинейных сеточных уравнений (9), (10) посредством метода Ньютона (c относительной погрешностью ${\text{rtol}} > 0$) и одномерных циклических двухточечных прогонок вдоль оси $Ox$. Устойчивость этих прогонок гарантируется знакоопределенностью матриц ${{{\mathbf{A}}}^{{x, + }}}({{{\mathbf{\hat {V}}}}_{{j,k,\ell }}})$. Обозначим нахождение решения ${\mathbf{\hat {V}}}$ по известному решению символически как действие оператора послойного перехода $S_{{{\text{conv}}}}^{{x, + }}(\tau )$: Поток с противоположным знаком ${{{\mathbf{F}}}^{{x, - }}}({\mathbf{v}})$, а также направления $y$, $z$ (т.е. все остальные одномерные системы (6)) рассматриваются аналогично; операторы перехода соответствующих им одномерных бикомпактных схем обозначим через $S_{{{\text{conv}}}}^{{x, - }}(\tau )$ и $S_{{{\text{conv}}}}^{{d, \pm }}(\tau )$, $d = y,z$.

Итоговая симметризованная по $x$, $y$ локально-одномерная бикомпактная схема для конвективной части (2) определяется своим оператором перехода ${{S}_{{{\text{conv}}}}}(\tau )$:

Данная схема аппроксимирует периодическую краевую задачу для уравнений (2) с погрешностью $O({{\tau }^{2}},{{h}^{4}})$ при $\tau ,h \to 0$, где $h$ – максимальный пространственный шаг.

Заметим, что для аппроксимации по времени совсем не обязательно использовать метод трапеций, вместо него, к примеру, можно выбрать любой другой $A$- или $L$-устойчивый диагонально-неявный метод Рунге–Кутты высокого порядка. Наше предпочтение методу трапеций обусловлено двумя факторами.

1. Порядок аппроксимации метода трапеций равен двум, что совпадает с порядком метода расщепления Марчука–Стрэнга.

2. Этот второй порядок экономично достигается лишь за одну неявную стадию.

1.5. Аппроксимация вязкой части уравнений Навье–Стокса

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

Пусть $u(x)$ – произвольная достаточно гладкая функция. Введем тогда сеточные функции ${{u}_{j}} = u({{x}_{j}})$ и ${{(u_{x}^{'})}_{j}} = du({{x}_{j}}){\text{/}}dx$, где индекс $j$ пробегает только целые значения. Обратимся к ячейке $x \in [{{x}_{j}},{{x}_{{j + 1}}}]$: имея на ее границах четыре значения: ${{u}_{j}},\;{{u}_{{j + 1}}},\;{{(u_{x}^{'})}_{j}},\;{{(u_{x}^{'})}_{{j + 1}}}$, мы можем построить по ним кубический интерполяционный полином Эрмита ${{u}_{h}}(x)$, который будет приближать функцию $u(x)$ внутри ячейки с погрешностью $O(h_{{x,j + 1/2}}^{4})$ при ${{h}_{{x,j + 1/2}}} \to 0$. Нетрудно показать, что для полинома ${{u}_{h}}(x)$ имеют место такие формулы (ради компактности записи опускаем индекс $j + 1{\text{/}}2$ у шага ${{h}_{{x,j + 1/2}}}$):

(член ${{h}_{x}}\tilde {A}_{0}^{x}{{u}_{{j + 1/2}}}$ совпадает с известной квадратурой Эйлера–Маклорена). Формулы (12) определяют сеточные операторы $\tilde {A}_{0}^{x}$, $\tilde {\Lambda }_{1}^{x}$, $\tilde {\Lambda }_{2}^{x}$, $\tilde {\Lambda }_{3}^{x}$, действующие на сеточные функции, заданные только в целых узлах и дополненные своими сеточными функциями-производными. Эти новые сеточные операторы представляются в виде комбинаций трех более простых операторов $M_{0}^{x}$, $\Delta _{0}^{x}$, ${{P}^{x}}$: Очевидно, оператор ${{P}^{x}}$ при действии на сеточную функцию ${{u}_{j}}$ заменяет ее значение на значение соответствующей ей сеточной функции-производной ${{(u_{x}^{'})}_{j}}$ в том же узле.

Вернемся к системе (3). Возьмем ее саму и семь ее дифференциальных следствий – производных по $x$, $y$, $z$, $xy$, $yz$, $xz$, $xyz$, осредним эти восемь векторных уравнений по ячейке $[{{x}_{j}},{{x}_{{j + 1}}}] \times [{{y}_{k}},{{y}_{{k + 1}}}] \times [{{z}_{\ell }},{{z}_{{\ell + 1}}}]$ (т.е. проинтегрируем и поделим на объем ячейки), сведем кратные интегралы к повторным, одномерные интегралы приблизим по формулам (12), для аппроксимации по времени применим метод трапеций. В результате перечисленных действий получится следующая бикомпактная схема для системы (3):

Сеточные функции ${{{\mathbf{\hat {v}}}}_{{j,k,\ell }}}$ и – это сеточные функции типа ${\mathbf{v}}_{{j,k,\ell }}^{n}$; первая означает некое промежуточное искомое решение, а вторая – некое известное решение (сказанное распространяется также на отвечающие им всем их сеточные функции-производные). Ради ясности перечислим все возможные варианты действия операторов ${{P}^{x}}$, ${{P}^{y}}$, ${{P}^{z}}$: причем все они полагаются коммутирующими между собой: Дискретные уравнения (13) дополняются периодическими граничными условиями, по сути представляющими из себя проекцию граничных условий для системы дифференциальных уравнений (3): Соотношения (14) являются достаточными, поскольку граничные условия для остальных искомых сеточных функций ${{({\mathbf{\hat {v}}}{{_{{}}^{{''}}}_{{xy}}})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}{{_{{}}^{{''}}}_{{yz}}})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}{{_{{}}^{{''}}}_{{xz}}})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}{{_{{}}^{{'''}}}_{{xyz}}})}_{{j,k,\ell }}}$ получаются путем подходящего применения операторов ${{P}^{x}}$, ${{P}^{y}}$, ${{P}^{z}}$ к равенствам (14) (например, к первым двум из них необходимо применять операторы ${{P}^{y}}$, ${{P}^{z}}$ и ${{P}^{y}}{{P}^{z}}$, что соответствует дифференцированию вдоль границы условий для системы (3)).

Схема (13), (14) аппроксимирует периодическую краевую задачу для системы (3) с погрешностью $O({{\tau }^{2}},{{h}^{4}})$. Важно отметить, что описанный выше метод получения бикомпактной пространственной аппроксимации четвертого порядка для многомерного уравнения теплопроводности обобщает метод работы [11], в которой был рассмотрен только одномерный случай.

Обозначим оператор послойного перехода бикомпактной схемы (13), (14) как ${{S}_{{{\text{visc}}}}}(\tau )$:

Здесь и далее, записывая “${\mathbf{\hat {v}}}$” без индексов, мы будем иметь в виду всю восьмерку сеточных функций ${{{\mathbf{\hat {v}}}}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}{{_{x}^{'}}_{x}})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}_{y}^{'})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}_{z}^{'})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}_{{xy}}^{{''}})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}_{{yz}}^{{''}})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}_{{xz}}^{{''}})}_{{j,k,\ell }}}$, ${{({\mathbf{\hat {v}}}_{{xyz}}^{{'''}})}_{{j,k,\ell }}}$ (аналогично для , ${{{\mathbf{v}}}^{n}}$).

Поясним теперь, как именно из уравнений (13), (14) вычисляется искомое решение ${\mathbf{\hat {v}}}$. Предлагается использовать метод итерируемой приближенной факторизации (ИПФ). Пусть ${{{\mathbf{\hat {v}}}}^{{(s)}}}$ – приближение с номером $s$ к искомому решению ${\mathbf{\hat {v}}}$ ($s = 0,\;1,\; \ldots $), начальное приближение ; обозначим через $\delta {{{\mathbf{\hat {v}}}}^{{(s + 1)}}} = {{{\mathbf{\hat {v}}}}^{{(s + 1)}}} - {{{\mathbf{\hat {v}}}}^{{(s)}}}$. Введем сеточные операторы:

где $\alpha > 0$ – параметр. Тогда уравнения метода ИПФ, записанные посредством операторов (16), имеют вид: где $\tau * = \tau {\text{/}}(2{\text{Re}})$, Граничные условия для $($17) получаются путем добавления “$\delta $” и верхнего индекса “$(s + 1)$” в соотношения (14). Искомое приращение $\delta {{{\mathbf{\hat {v}}}}^{{(s + 1)}}}$ находится из уравнений (17) обращением пар одномерных операторов ($B_{1}^{d}(\tau {\kern 1pt} *{\kern 1pt} ),B_{2}^{d}(\tau {\kern 1pt} *{\kern 1pt} ))$ в такой последовательности: по $z$, по $y$, по $x$. В свою очередь, обращение каждой такой пары на всей сетке сводится к совокупности независимых одномерных циклических двухточечных прогонок, устойчивых коль скоро $\tau {\kern 1pt} * > 0$, т.е. всегда. Конечно, тот факт, что векторное уравнение (3) распадается на три независимых скалярных уравнения для компонент скорости, отражается и в схеме (13) , и в методе ИПФ (17), благодаря чему все вычисления можно выполнять отдельно для каждой компоненты скорости независимо от других компонент (т.e., например, сначала находится ${{{\hat {v}}}_{x}}$, далее ${{{\hat {v}}}_{y}}$, потом ${{{\hat {v}}}_{z}}$). Итерирование кончается, если выполнено условие сходимости

1.6. Аппроксимация части уравнений Навье–Стокса, отвечающей за силы давления и несжимаемость жидкости

Здесь мы будем действовать иначе – сперва дискретизируем уравнения (4) по времени при помощи неявного методом Эйлера:

Затем, взяв дивергенцию от обеих частей первого уравнения и воспользовавшись вторым уравнением, получим два независимых уравнения: Таким образом, численное решение (4) разбивается на две отдельные подзадачи: решение уравнения Пуассона для давления (18) и пересчет скорости жидкости по формуле (19).

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

Уравнение Пуассона (18), как и полная система уравнений Навье–Стокса, решается при периодических граничных условиях, проекция которых на сетку дает дискретные граничные условия для схемы (20) :

Бикомпактная схема (20), (21) аппроксимирует периодическую краевую задачу для уравнения Пуассона (18) с погрешностью $O({{h}^{4}})$. Заметим, что в погрешности нет члена $O({{h}^{4}}{\text{/}}\tau )$: хотя дивергенция скорости делится на $\tau $, она сама обязана иметь асимптотику $O(\tau )$, так как при сгущении сетки численное распределение давления сходится к точному и составляет конечную величину $O(1)$.

Для реализации схемы (20) применяется метод ИПФ. Уравнения этого метода для нее совершенно аналогичны тем, что были приведены выше для схемы (13) . Разница состоит лишь в том, что в качестве аргумента у операторов $B_{{1,2}}^{d}({\kern 1pt} \cdot {\kern 1pt} )$, $d = x,y,z$, выступает не модифицированный временной шаг $\tau {\kern 1pt} {\text{*}}$, а параметр $\theta > 0$ – настраиваемый параметр итерационного метода.

Обозначим нахождение давления $\hat {p}$ по схеме (20), (21) символически как действие оператора ${{S}_{{{\text{ell}}}}}(\tau )$ на :

Обратимся к пространственной дискретизации соотношения (19). Для компоненты скорости ${{{v}}_{x}}$ имеем

Проекция на сетку соотношения (23) и его дифференциальных следствий – производных по $y$, $z$, $yz$ дает следующие формулы: Легко видеть, что формулы (24) позволяют найти лишь четыре из восьми искомых сеточных функций. Оставшиеся четыре функции должны вычисляться из еще не рассмотренных дифференциальных следствий (23) – его производных по $x$, $xy$, $xz$, $xyz$. Ясно, что эти следствия нельзя спроецировать на сетку, так как тогда мы, например, столкнемся со второй производной от давления по $x$, которая не содержится непосредственно в решении для $\hat {p}$. Вспомним пример с функцией $u(x)$, в рамках которого выведены формулы (12). Аппроксимируем значение $u_{{xx}}^{{''}}({{x}_{j}})$ с погрешностью $O(h_{x}^{3})$ на трехточечном шаблоне $\{ {{x}_{{j - 1}}},{{x}_{j}},{{x}_{{j + 1}}}\} $ по значениям сеточных функций ${{u}_{j}}$ и ${{(u_{x}^{'})}_{j}}$ и обозначим эту аппроксимацию через $D_{2}^{x}{{u}_{j}}$: На равномерной сетке оператор $D_{2}^{x}$ выражается в виде а формула (25) повышает свою точность до $O(h_{x}^{4})$. Используя трехточечную аппроксимацию (25), напишем расчетные формулы для оставшейся половины $\mathop {{\hat {v}}}\nolimits_x $: Заметим, что диапазоны изменения индексов $j$, $k$, $\ell $ в формулах (26) приведены корректно: например, значение ${{\hat {p}}_{{ - 1,k,\ell }}}$ заменяется на ${{\hat {p}}_{{{{N}_{x}} - 1,k,\ell }}}$ в силу периодичности вычисляемого решения задачи.

Что касается дискретизации соотношения (19) применительно к ${{{v}}_{y}}$ и ${{{v}}_{z}}$, то она выводится в полной аналогии с рассуждениями предыдущего абзаца, поэтому выписывать варианты формул (24), (26) для ${{{v}}_{y}}$ и ${{{v}}_{z}}$ мы не будем. Обозначим нахождение ${\mathbf{\hat {v}}}$ по формулам (24), (26) и их аналогам для ${{{v}}_{y}}$ и ${{{v}}_{z}}$ символически как действие оператора перехода ${{S}_{{{\text{pr}}}}}(\tau ,\hat {p})$:

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

Для составления полного алгоритма счета осталось рассмотреть и описать две операции: пересчет численного решения ${{{\mathbf{v}}}^{n}}$ в решение ${{{\mathbf{V}}}^{n}}$ и наоборот.

Переход от ${{{\mathbf{v}}}^{n}}$ к ${{{\mathbf{V}}}^{n}}$ осуществляется достаточно просто: решение ${{{\mathbf{v}}}^{n}}$ содержит всю информацию, необходимую для построения в каждой ячейке сетки трикубического (то есть трехмерного аналога бикубического) интерполяционного полинома Эрмита. Отметим, что совокупность этих полиномов представляет из себя трикубический сплайн гладкости 1. Этот сплайн позволяет посчитать скорость жидкости в любой точке $\bar {D}$ с четвертым порядком точности, в частности, в полуцелых узлах сетки $\Omega $. Следуя нашим обозначениям, обозначим пересчет ${{{\mathbf{v}}}^{n}}$ в ${{{\mathbf{V}}}^{n}}$ через действие оператора ${{S}_{{{\mathbf{v}} \to {\mathbf{V}}}}}$:

Переход от ${{{\mathbf{V}}}^{n}}$ к ${{{\mathbf{v}}}^{n}}$ выполняется при помощи пятиточечной конечно-разностной аппроксимации для первой производной с третьим порядком точности. На равномерной сетке формулы пересчета ${{{\mathbf{V}}}^{n}}$ в ${{{\mathbf{v}}}^{n}}$ имеют четвертый порядок точности и записываются в виде

Поставим в соответствие этим формулам пересчета оператор ${{S}_{{{\mathbf{V}} \to {\mathbf{v}}}}}$:

1.8. Результирующая бикомпактная схема и общий алгоритм счета

Введенные выше операторы перехода (11), (15), (22), (27), (28) и (30) позволяют записать результирующую бикомпактную схему в максимально лаконичном виде:

Бикомпактная схема (31) аппроксимирует систему уравнений Навье–Стокса (1) с погрешностью $O({{\tau }^{2}},{{h}^{3}})$ на неравномерных пространственных сетках и с погрешностью $O({{\tau }^{2}},{{h}^{4}})$ на равномерных сетках. В любом случае, при $\tau = O(h)$ погрешность аппроксимации составляет $O({{\tau }^{2}})$. Отметим, что при такой связи между шагами $\tau $ и $h$ максимальное конвективное (гиперболическое) число Куранта остается постоянным при сгущении сетки, что вполне уместно и желательно при расчетах течений с доминированием процесса конвекции, т.е. при больших числах Рейнольдса.

Далее схему (31) мы будем сокращенно называть T2B4. Тем самым мы подчеркиваем, что для счета по времени применяется метод трапеций второго порядка, а для дискретизации пространственных производных в уравнениях (2), (3), (18) – бикомпактная аппроксимация четвертого порядка.

Дадим также словесное описание общего алгоритма счета по схеме T2B4. Переход со слоя ${{t}^{n}}$ на слой ${{t}^{{n + 1}}}$ делится на 9 этапов.

Шаг 1. Вычислить параметры потокового расщепления $C_{2}^{d}$, $d = x,y,z$, по формуле (8).

Шаг 2. Добавить к скорости ${{{\mathbf{v}}}^{n}}$ слагаемое $( - (\tau {\text{/}}2)\nabla {{p}^{n}})$ по формулам (24), (26) и их вариантам для ${{{v}}_{y}}$, ${{{v}}_{z}}$.

Шаг 3. Сделать полшага ($\tau {\text{/}}2$) по бикомпактной схеме (13), (14) для вязкой части уравнений Навье–Стокса, взяв в качестве скорость, найденную на шаге 2.

Шаг 4. Построить в каждой ячейке сетки трикубический интерполяционный полином Эрмита для скорости, найденной на шаге 3. При помощи этого полинома посчитать значения скорости в полуцелых узлах сетки $\Omega $, отбросив после этого значения производных скорости в целых узлах.

Шаг 5. Сделать полный шаг ($\tau $) по локально-одномерной бикомпактной схеме (11) для конвективной части уравнений Навье–Стокса, взяв вместо скорость, полученную на шаге 4.

Шаг 6. Посчитать по формулам (29) значения производных скорости, найденной на шаге 5, в целых узлах сетки $\Omega $, отбросив после этого значения скорости в полуцелых узлах.

Шаг 7. Сделать полшага ($\tau {\text{/}}2$) по бикомпактной схеме (13), (14) для вязкой части уравнений Навье–Стокса, взяв в качестве скорость, найденную на шаге 6.

Шаг 8. Найти по бикомпактной схеме (20), (21) искомое давление на верхнем слое ${{p}^{{n + 1}}}$, заменив при этом шаг $\tau $ на $\tau {\text{/}}2$, а скорость на скорость, найденную на шаге 7.

Шаг 9. Получить искомую скорость на верхнем слое ${{{\mathbf{v}}}^{{n + 1}}}$, добавив к скорости, найденной на шаге 7, слагаемое $( - (\tau {\text{/}}2)\nabla {{p}^{{n + 1}}})$ по формулам (24), (26) и их вариантам для ${{{v}}_{y}}$, ${{{v}}_{z}}$.

3.2. Результаты для Re = 400

Возьмем сетку по пространству ${{64}^{3}}$ ячеек, положим временной шаг $\tau = 0.02$ (конвективное число Куранта равно 0.5). Посчитанное поле модуля завихренности $\omega = \left| {\mathbf{\omega }} \right|$ изображено в моменты времени $t = 0$, 5, 10, 15, 17.5, 20 на фиг. 1. Видно, что распад начального вихря действительно происходит в ламинарном режиме: активного перемешивания слоев жидкости не наблюдается, завихренность по своей структуре остается упорядоченно ячеистой, с течением времени (при $t > 10$) вытягиваясь в длинные филаменты (волокна) и затухая. На фиг. 2 показаны две кривые скорости диссипации кинетической энергии $\varepsilon = - d{{E}_{K}}{\text{/}}dt$: одна получена по схеме T2B4, другая взята из работы [24] и полагается эталонной. Хотя в окрестности $t = 8$ схема T2B4 несколько занижает диссипацию ${{E}_{K}}$, в целом можно сказать, что на достаточно грубой сетке ${{64}^{3}}$ ячеек схема T2B4 воспроизводит эту зависимость близко к эталону.

3.3. Результаты для Re = 800

Возьмем более подробную сетку по пространству ${{128}^{3}}$ ячеек, временной шаг пропорционально уменьшим до $\tau = 0.01$. Визуализация найденного поля $\omega $ в разные моменты времени приведена на фиг. 3. Из этой визуализации следует, что течение сохраняет более-менее упорядоченную ячеистую структуру вплоть до момента времени $t = 17.5$; однако, ко времени $t = 20$ движение жидкости становится полностью турбулентным, как и ожидалось. Из данных с фиг. 4 можно заключить, что при ${\text{Re}} = 800$ схема T2B4, как и в предыдущем варианте задачи при ${\text{Re}} = 400$, дает верную скорость диссипации кинетической энергии.

3.4. Результаты для Re = 1600

Параметры сетки выберем такими же, как в предыдущем варианте: по пространству ${{128}^{3}}$ ячеек, по времени шаг $\tau = 0.01$. Отметим, что в данном варианте задачи на данной сетке “чистая” (без монотонизации) бездиссипативная бикомпактная схема T2B4 “разваливается” в момент времени $t \approx 7.4$, сразу после начала интенсивного вихреобразования. Эта проблема устраняется, если оператор перехода ${{S}_{{{\text{conv}}}}}(\tau )$ (11) лимитирован при помощи метода консервативной монотонизации [21]. Положим параметр этого метода ${{C}_{1}} = 0.5$, что на порядок меньше тех значений, которые обычно употребляются при расчетах задач с сильными разрывами (см. [21]).

На фиг. 5 приведена визуализация поля $\omega $ в разные моменты времени. Что разумно, течение при ${\text{Re}} = 1600$ по сравнению с течением при ${\text{Re}} = 800$ характеризуется более сильным производством завихренности; уже к моменту времени времени $t = 15$ всякая связь с начальной структурой движения потеряна, а при $t = 20$ вычислительная область целиком заполнена практически однородной изотропной турбулентностью.

Перейдем к анализу кривых диссипации кинетической энергии с фиг. 6. Из этой фигуры ясно, что при ${\text{Re}} = 1600$ монотонизированная схема T2B4 воспроизводит зависимость $\varepsilon (t)$ хуже, чем ее чистая версия при ${\text{Re}} = 400$, 800. Тем не менее нельзя не заметить две особенности. Во-первых, монотонизированная схема T2B4 дает кривую диссипации, примерно правильную по форме, но смещенную влево на $\Delta t \approx 0.3$. Другими словами, процесс превращения кинетической энергии в тепло идет по траектории верной формы, но с небольшим опережением. Во-вторых, если судить по данным для “чистой” схемы T2B4 (очевидно, доступным только до момента ее развала), то метод консервативной монотонизации не дает существенной численной прибавки к реальной диссипации. Из второго замечания следует, что опережение в диссипации кинетической энергии связано с подробностью сетки и погрешностями чистой бикомпактной аппроксимации, а не с методом консервативной монотонизации.

Важнейшим требованием к любой численной схеме для счета турбулентных течений является то, что она должна как можно точнее воспроизводить спектр кинетической энергии ${{\sigma }_{K}}(k,t)$, где $k$ – волновое число, пробегающее значения $0,\;1,\; \ldots ,\;{{N}_{{min}}}{\text{/}}2$, ${{N}_{{min}}} = min({{N}_{x}},{{N}_{y}},{{N}_{z}})$. Уточним, что ${{\sigma }_{K}}(k,t)$ вычисляется по формуле (например, см. [24]):

где ${\mathbf{k}} = ({{k}_{x}},{{k}_{y}},{{k}_{z}})$ – вектор волновых чисел, ${{k}_{d}} = 0,\; \pm 1,\; \ldots ,\; \pm {{N}_{d}}{\text{/}}2$, $d = x,y,z$, ${\mathbf{A}}({\mathbf{k}},t)$ – коэффициент ряда Фурье скорости ${\mathbf{v}}({\mathbf{r}},t)$: Коэффициент ${\mathbf{A}}({\mathbf{k}},t)$ вычисляется по значениям только сеточной функции ${{{\mathbf{v}}}_{{j,k,\ell }}}$ при помощи квадратурной формулы трапеций и алгоритма быстрого преобразования Фурье.

Фигура 7 позволяет выяснить, насколько хорошо этому требованию удовлетворяет монотонизированная бикомпактная схема T2B4. Сравнение с эталонным спектром, полученным в работе [29] по псевдоспектральной схеме с очень высоким разрешением в ${{512}^{3}}$ степеней свободы (фактически это спектр сошедшегося решения), показывает, что схема T2B4 безукоризненно разрешает инерционный интервал, но несколько преувеличивает затухание мелких вихрей в диссипативном интервале. Впрочем, у схемы T2B4 диссипативный участок спектра не сваливается резко вниз, а идет почти параллельно эталону, что свидетельствует об умеренности численной диссипации, вносимой методом консервативной монотонизации [21].

Интересно сделать более серьезную проверку бикомпактной схемы T2B4 на разрешение спектра: выполним расчеты на менее подробной сетке из ${{64}^{3}}$ ячеек и на совсем редкой сетке из ${{32}^{3}}$ ячеек ($\tau = 0.02$ и 0.04 соответственно). Консервативную монотонизацию при этом выключим (${{C}_{1}} = 0$), так как схема T2B4 на этих сетках досчитывает задачу до конца без монотонизации. На фиг. 8 изображены полученные спектры. Очевидно, схема T2B4 корректно, без “загибов вниз”, воспроизводит инерционный участок спектра даже на сетке из ${{32}^{3}}$ ячеек, не разрешающей все мелкие вихри и часть вихрей среднего масштаба. Из этого можно сделать вывод, что бикомпактная схема T2B4 позволяет правильно рассчитывать динамику крупных вихрей, что очень важно при LES-моделировании.