Журнал вычислительной математики и математической физики, 2020, T. 60, № 12, стр. 2050-2054

Одномерное уравнение Бюргерса представлено в виде

где u – вычисляемый параметр, изменяющийся по времени t и в направлении х, Re – критерий Рейнольдса.

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

Разностные методы первого порядка аппроксимации по временному ∆t и пространственному ∆х шагам интегрирования испытывают сильное влияние схемной вязкости, но обладают свойством монотонности, сохраняющим знак первой производной [1].

Для численного решения некоторых задач бывает достаточно получить результат со вторым порядком аппроксимации (∆t², ∆х²). Например, при использовании схемы Лакса–Вендроффа необходимо определять параметры течения в полуцелых пространственных и временных узлах разностной сетки, а для учета второй производной в (1) надо брать значения параметра не в трех, а четырех соседних пространственных узлах сетки.

В предлагаемой схеме на первом этапе, предикторе, для разностной схемы Мак-Кормака вводится параметр θ, как и в схеме Кранка–Никольсона. Для j-го узла сетки по координате х с n-го временного слоя значение скорости ${{\bar {u}}_{j}}$ определяется из уравнения

где введено обозначение $\gamma = \frac{{\Delta t}}{{\operatorname{Re} \Delta x}}$. Вид члена с первой производной в (2) периодически меняется с представленного в (2) на $\frac{{{{{(u_{j}^{n})}}^{2}} - {{{(u_{{j - 1}}^{n})}}^{2}}}}{{2\Delta x}}\Delta t$. Одна из причин состоит в устранении влияния направления, в котором она вычисляется. При решении трехмерной задачи в [2] дана таблица с возможными периодическими изменениями в определении первых пространственных производных. Целесообразно нелинейный член в (1) представлять в (2) в дивергентной форме.

Второй этап, корректор, основан на процедуре перешагивания (“чехарда”, “leap-frog”) [3]. Уравнение (1) решается для j-го пространственного узла между n-м и n + 1-м временными слоями. Нелинейный член определен как среднее значение, взятое на n-м временном слое и на этапе предиктор, но с противоположным направлением в дифференцировании, как и в методе Мак-Кормака. При вычислении члена со второй производной в (1) значение скорости в j-м узле берется на n + 1-м временном слое, а остальные на n-м слое [4]. В результате получаем

Данная схема, как будет показано, имеет второй порядок аппроксимации и, как следствие, не обладает свойством монотонности. Поэтому, если модуль величины $\frac{{\partial u}}{{\partial x}}$ в (1) имеет большое значение, то знак производной может отличаться от того, который есть на самом деле. Для предотвращения этого явления рекомендуется после каждого этапа корректора проводить “градиентное” сглаживание, решая уравнение

Для каждой узловой точки с координатой ${{x}_{j}}$ величина

где $\left( {\frac{{\partial u}}{{\partial x}}_{j}^{{n + 1}}} \right)_{{{\text{max}}}}^{2}$ – максимальное среди всех значений ${{\left( {\frac{{\partial u}}{{\partial x}}_{j}^{{n + 1}}} \right)}^{2}}$. Значение коэффициента сглаживания α принимается чуть больше нуля. Обычно $0 < \alpha \leqslant 0.2$. Начальные и граничные значения ${{u}_{s}}$ берутся равными соответствующим величинам $u_{j}^{{n + 1}}$.

Полученные при решении (4) значения ${{u}_{s}}$ принимаются равными $u_{j}^{{n + 1}}$ на n + 1 временном слое.

Порядок аппроксимации первой производной по времени после этапа корректора равен ∆t² из-за центральной разности. Периодическое чередование направления численного дифференцирования при $n \to \infty $ приведет в среднем к центральной разности второго порядка (∆х²). Правый член в (1) имеет второй порядок аппроксимации по пространственной координате и такой же порядок по времени. Это можно получить, проведя в правой части уравнения (3) разложение в ряд Тейлора в j-м узле после этапа коррекции.

Анализ устойчивости проводится для линеаризованного уравнения (1) заменой нелинейного члена на ${{u}_{0}}\frac{{\partial u}}{{\partial x}}$, где ${{u}_{0}} = {\text{const}}$. Вводится число Куранта $c = \frac{{\left| {{{u}_{0}}} \right|\Delta t}}{{\Delta x}}$. Для комплексного фурье-разложения пространственных возмущений с частотой ω вводится угол φ = ω∆х. Модуль коэффициента перехода $q = \left| {\frac{{\delta u_{j}^{{n + 1}}}}{{\delta u_{j}^{n}}}} \right|$ равен модулю отношения ошибки $\delta u_{j}^{{n + 1}}$ на n + 1-м временном слое к ошибке $\delta u_{j}^{n}$. Чередование направления численного дифференцирования дает право заменить $1 - {{e}^{{ - i\varphi }}}$ на $\frac{{{{e}^{{i\varphi }}} - {{e}^{{ - i\varphi }}}}}{2} = i\sin {\varphi }$. Анализ проведен без учета влияния сглаживания.

Тогда из (2) и (3) получим $Z\delta u_{j}^{{n + 1}} = (A - iB)\delta u_{j}^{n}$, где

Имеем После преобразований получим

Из (4) следует, что достаточно выполнения условий $\theta \geqslant 0.5$ и ${{с}^{2}} \leqslant 1$, чтобы ${{q}^{2}} \leqslant 1$, и предложенная схема будет устойчива для любых значений Re в уравнении (1).Тогда шаг интегрирования по времени получается из условия $\Delta {{t}^{2}} \leqslant {{\left( {\frac{{\Delta x}}{{{{u}_{{{\text{max}}}}}}}} \right)}^{2}}$. Чтобы учесть сглаживание, необходимо умножить полученное значение q² на квадрат коэффициента перехода уравнения (5), который ввиду малости значения α не превысит единицу.

Во всех приводимых примерах расчета θ = 0.55 значение числа Куранта не превышало 0.95 для всех узлов сетки. Расчет заканчивался, когда $t \leqslant {{t}_{k}}$, после чего вычислялась среднеквадратичная ошибка

относительно N значений точного решения задачи u_Т (см. [5]). Во всех примерах рассматривается процесс распространения возмущений в положительном направлении.

Первый пример расчета для уравнения (1) взят для $ - 1 \leqslant х \leqslant 1$, Re = 10³, t_к = 0.92 при начальном условии для t = 0

На границе расчетной области u(x = –1) = 1, u(x = 1) = 0. Для учащения сетки [5] вблизи х = 0 введена замена ξ = arctg(x). Значения σ по результатам расчета приведены в табл. 1 для различных значений α и количества отрезков интегрирования k при $х \geqslant 0$. Наилучшее совпадение с точным решением, как следует из табл. 2, наблюдается при k = 10 и α = 0.1. В этом случае при резком, скачкообразном, изменении параметра течения, наблюдается наибольшая длина отрезка, где $c \approx 1$ и получены наиболее точные результаты счета. Но при этом сильнее выражены осцилляции на фронте волны (фиг. 1).

Для второго примера $0 \leqslant x \leqslant 1$, Re = 10⁴, t_k = 0.5, начальные условия $u = \sin \left( {\pi x} \right)$ при t = 0. На границе u(x = 0) = u(x = 1) = 0. Разбиение сетки было равномерным. Резкое изменение значения u наблюдалось вблизи х = 1. При увеличении k точность повышается до определенного значения k, так как в дальнейшем наступает неустойчивость дифференцирования. На фиг. 2 показаны результаты расчета.