Журнал вычислительной математики и математической физики, 2022, T. 62, № 1, стр. 175-192

2.1. Уравнение неразрывности

Скажем несколько слов об используемых сетках. В расчетной области строятся узлы, лежащие внутри области или на ее границе. Затем область разбивается на ячейки. Каждый узел принадлежит своей ячейке. Ячейками могут служить ячейки Дирихле.

Идею метода рассмотрим на примере уравнения неразрывности

Проинтегрируем уравнение по k-й ячейке и заменим интегралы от пространственных производных на поверхностные интегралы, когда граница ячейки разбита на участки, соответствующие границе между k-й и соседней m-й ячейками:

где ${{n}_{x}}$, ${{n}_{y}}$, ${{n}_{z}}$ – проекции нормали к границе между ячейками. Для аппроксимации функций, входящих в подынтегральные выражения, будем использовать полиномы от пространственных переменных. Плотность в первом слагаемом будем аппроксимировать полиномом второго порядка, а выражения $\rho u$, $\rho {v}$, $\rho w$ – полиномами третьего порядка.

2.2. Построение полиномов

Пусть у нас есть $I$ узлов, в которых нам известны значения функции $f(x,y,z)$. Обозначим эти значения ${{f}_{i}}$. Полином $F(x,y,z)$, аппроксимирующий функцию $f,$ имеет вид

где базисные функции ${{\psi }_{j}}(x,y,z) = {{X}^{p}}{{Y}^{q}}{{Z}^{s}}$, $X = x - {{x}_{c}}$, $Y = y - {{y}_{c}}$, $Z = z - {{z}_{c}}$, $p + q + s \leqslant r$. В случае построения полинома второго порядка для аппроксимации плотности $r = 2$, а ${{x}_{c}},~\;{{y}_{c}},~\;{{z}_{c}}$ совпадают с координатами k-го узла и ${{\psi }_{j}} = \{ 1,X,Y,Z,{{X}^{2}},{{Y}^{2}},{{Z}^{2}},XY,XZ,YZ\} $, $J = 10$. В дальнейшем будем обозначать функции ${{\psi }_{j}}$ как$~\psi _{j}^{k}$.

В случае построения полинома третьего порядка для аппроксимации $\rho u$, $\rho {v}$, $\rho w$ на границе между k-й и m-й ячейками $r = 3$, а ${{x}_{c}},~\;{{y}_{c}},\;~{{z}_{c}}$ – середина отрезка, соединяющего k-й и m-й узлы. Функции ${{\psi }_{j}}$ имеют вид ${{\psi }_{j}} = \{ 1,X,Y,Z,{{X}^{2}},{{Y}^{2}},{{Z}^{2}},XY,XZ,YZ,{{X}^{3}},{{Y}^{3}},{{Z}^{3}},{{X}^{2}}Y,{{X}^{2}}Z,{{Y}^{2}}X,{{Y}^{2}}Z,{{Z}^{2}}X,$ ${{Z}^{2}}Y,XYZ\} $, $J = 20.~$ В этом случае для ${{\psi }_{j}}$ будем использовать выражение $\psi _{j}^{{km}}$.

В дальнейшем мы считаем, что $I > J$. Полагаем, что вектор коэффициентов полиномов ${\mathbf{c}}$ с компонентами ${{c}_{j}}$ линейно связан со значениями функции $f$ в соответствующих узлах (${{f}_{i}}$): ${\mathbf{c}} = A{\mathbf{f}}$, $A$ – прямоугольная матрица из $J$ строк и $I$ столбцов. Чтобы найти матрицу $A$, потребуем выполнения следующего условия: если $f = {{\psi }_{j}}$, то коэффициенты ${{с}_{s}} = {{\delta }_{{sj}}}$ (${{с}_{j}} = 1$, остальные компоненты вектора $c$ равны нулю). Это условие можно записать в виде матричного уравнения

или где ${{\Psi }}$ – прямоугольная матрица, состоящая из $I$ строк и $J$ столбцов, а ее коэффициенты равны ${{{{\Psi }}}_{{{\text{ij}}}}} = {{\psi }_{j}}({{x}_{i}},{{y}_{i}},{{z}_{i}})$, а ${{E}_{J}}~$ – единичная матрица размером $J$. Одно из решений этого уравнения имеет вид

Легко убедиться, что такое решение удовлетворяет уравнению (2). Именно такое решение мы будем использовать при построении полиномов третьего порядка, аппроксимирующих функции $\rho u$, $\rho {v}$, $\rho w$.

Для построения квадратичного полинома, аппроксимирующего плотность, мы будем выбирать в некотором смысле наилучшее из решений уравнения (2). Отметим, что если ${{A}_{1}}$ и ${{A}_{2}}$ – решения уравнения (2), то их разность ${{A}_{d}} = {{A}_{1}} - {{A}_{2}}$ удовлетворяет уравнению ${{A}_{d}}{{\Psi }} = 0$. Тогда любая строка ${{z}^{{\text{T}}}}$ матрицы ${{A}_{d}}$ удовлетворяет уравнению ${{z}^{{\text{T}}}}{{\Psi }} = 0$ или линейной системе ${{{{\Psi }}}^{{\text{T}}}}z = 0$. Все решения системы можно выписать, используя метод Гаусса с выбором главного элемента. Сделаем это для простого случая. Обозначим первые $J$ столбцов матрицы ${{{{\Psi }}}^{{\text{T}}}}$, как матрицу ${{{{\Psi }}}_{A}}$. А последние $I - J$ столбцов этой матрицы, как матрицу ${{{{\Psi }}}_{B}}$. Тогда $({{{{\Psi }}}_{A}}\,{\text{|}}\,{{{{\Psi }}}_{B}})z = 0$. Здесь и далее вертикальная черта делит матрицу на две части. Левая и правая части задаются соответствующими выражениями. Предположим, что матрицу ${{{{\Psi }}}_{A}}$ можно обратить (этого можно добиться перестановкой столбцов в ${{{{\Psi }}}^{{\text{T}}}}$ или точек в шаблоне). Наше уравнение можно переписать в виде ${{{{\Psi }}}_{A}}({{{\text{E}}}_{J}}\,{\text{|}}\,{{\Psi }}_{A}^{{ - 1}}{{{{\Psi }}}_{B}})z = 0$ или $\left( {{{{\text{E}}}_{J}}\,{\text{|}}\,{\text{M}}} \right)z = 0$, где ${\text{M}} = {{\Psi }}_{A}^{{ - 1}}{{{{\Psi }}}_{B}}$.

Решение этого уравнения имеет вид $z = \left( {\begin{array}{*{20}{c}} { - M\lambda } \\ \lambda \end{array}} \right)$, где вектор $\lambda $ состоит из $I - J$ произвольных чисел. Это же решение можно записать в виде $z = Z\lambda $, где $Z = \left( {\begin{array}{*{20}{c}} { - M} \\ {{{E}_{{I - J}}}} \end{array}} \right)$ – матрица из $I$ строк и $I - J$ столбцов.

Тогда любое решение уравнения (2) можно представить в виде

где матрица ${{\Lambda }}$ имеет $I - J$ строк и $J$ столбцов и состоит из произвольных чисел.

2.3. Шаблоны

При построении полиномов мы будем использовать два вида шаблонов – множеств узлов, обозначим их ${{S}_{A}}(k)$ и ${{S}_{B}}(k)$. Индекс $k$ в этих обозначениях – номер центрального узла шаблона. Назовем два узла соседними по ребру, если соответствующие им ячейки имеют общее ребро. Шаблон ${{S}_{A}}(k)$ включает k-й узел и соседние с ним по ребру.

Назовем два узла соседними по грани, если соответствующие им ячейки имеют общую грань. Пусть $M$ – множество узлов, через $adj(M)$ обозначим множество узлов, включающее узлы множества $M$ и узлы, соседние с ними по грани. Используя это обозначение, можно определить шаблон ${{S}_{B}}(k) = adj(adj(adj(k))),$ т.е. шаблон ${{S}_{B}}(k)$ включает соседние узлы по грани до третьей степени соседства с k-м узлом.

Если в качестве узлов используются узлы декартовой сетки, то шаблон ${{S}_{A}}$ будет содержать 19 узлов, а шаблон ${{S}_{B}}$ – 63 узла. На фиг. 1 показаны шаблон ${{S}_{A}}$ и верхняя половина шаблона ${{S}_{B}}$.

2.4. Аппроксимация интегралов

Далее будем использовать следующие выражения. Выражение $i \in S$ (где $S$ – шаблон) означает, что индекс пробегает номера точек, принадлежащие шаблону $S$. Выражение $n(i)$ означает функцию, дающую номер узла в шаблоне. Значения функции – целые числа от 1 до числа точек в шаблоне.

Плотность под интегралом в первом слагаемом в уравнении (1) аппроксимируется полиномом второго порядка ${{R}_{k}}(x,y,z)$:

где матрица ${{A}_{k}}$ удовлетворяет уравнению (2) и имеет вид (4), $J1 = 10$. Тогда интеграл в первом слагаемом в уравнении (1) можно аппроксимировать следующим образом: где ${{B}_{{kn(i)}}}$ имеют вид Чтобы вычислить объемные интегралы в этом выражении, k-я ячейка разбивается на пирамиды  с  основанием  на  грани ячейки и вершиной пирамиды в k-м узле. Затем основание пирамиды разбивается на треугольники, и вычисление объемного интеграла по ячейке сводится к вычислению объемных интегралов по нескольким тетраэдрам от функций $\psi _{j}^{k}(x,y,z) = {{(x - {{x}_{k}})}^{p}}{{(y - {{y}_{k}})}^{q}}{{(z - {{z}_{k}})}^{s}}$, $p + q + s \leqslant 2$.

Произведение $\rho u$ во втором слагаемом уравнения (1) мы аппроксимируем полиномом третьего порядка ${{F}_{{km}}}(x,y,z)$:

где ${{S}_{B}}(k,m) = {{S}_{B}}(k) \cap {{S}_{B}}(m)$.

Матрица ${{A}_{{km}}}$ является решением уравнения (2) и имеет вид (3), $J2 = 20$. Поверхностный интеграл от величины, включающей это выражение, имеет вид

При интегрировании граница между двумя ячейками разбивается на треугольники. Мы получили поток из k-й ячейки в m-ю. Такой же поток по величине, но с обратным знаком, мы получим при построении аппроксимации в m-й ячейке. Это обеспечивает консервативность нашего метода.

Учитывая суммирование по $m$ в (1), можно записать

Аналогичные выражения можно записать, аппроксимируя $\rho {v}$ и $\rho w$ полиномами третьего порядка ${{G}_{{km}}}(x,y,z)$ и ${{H}_{{km}}}(x,y,z)$: и Отметим, что во всех формулах используется одна и та же матрица ${{A}_{{km}}}$.

Принимая во внимание суммирование по $m$, получаем следующие формулы:

При реализации метода коэффициенты ${{B}_{{kn(i)}}}$, $i \in {{S}_{A}}(k)$, и ${{c}_{{kn(i)}}}$, ${{d}_{{kn(i)}}}$, ${{e}_{{kn(i)}}}$, $i \in {{S}_{B}}(k)$, могут быть вычислены заранее и храниться в памяти компьютера.

Используем величину

для оценки точности полиномов в окрестности k-го узла $\left| {\rho - {{R}_{k}}} \right| < O(h_{k}^{3})$, $\left| {\rho u - {{F}_{{km}}}} \right| < O(h_{k}^{4})$, $\left| {\rho {v} - {{G}_{{km}}}} \right| < O(h_{k}^{4})$, $\left| {\rho w - {{H}_{{km}}}} \right| < O(h_{k}^{4})$.

Считая объем k-й ячейки порядка $h_{k}^{3}$, а площадь ее граней порядка $h_{k}^{2}$, можно оценить точность аппроксимации интегралов, входящих в (1):

где $N = n(i)$.

2.5. Схема для уравнения (1)

Используем полученные выше аппроксимации для записи явной схемы, аппроксимирующей уравнение (1) с первым порядком точности по времени:

где $\tau $ – шаг по времени, а ${{\hat {\rho }}_{i}}$ – значение плотности на следующем временном слое, $N = n(i)$.

Слагаемые в этом уравнении аппроксимируют интегралы с точностью $O(h_{k}^{6})$. Величины ${{B}_{{kN}}}$ имеют порядок, сопоставимый с объемом k-й ячейки, т.е. порядка $h_{k}^{3}$. Таким образом, можно говорить, что уравнение (1) аппроксимируется с точностью $O(h_{k}^{3} + \tau )$.

Слагаемое ${{D}_{4}}\rho $ положительное и обеспечивает устойчивость схемы. Оператор ${{D}_{4}}\rho = {{V}_{k}}{{D}_{2}}[{{\mu }_{4}}{{D}_{2}}[\rho ]]$, где ${{V}_{k}}$ – объем k-й ячейки, ${{\mu }_{4}}$ – константа. Оператор ${{D}_{2}}[\rho ]$ имеет вид

где ${{S}_{{km}}}$ – площадь грани между k-й и m-й ячейками, а ${{h}_{{km}}}$ – сумма расстояний от k-го и m-го узлов до плоскости, содержащей грань (границу) между k-й и m-й ячейками. Отметим, что если $f({{h}_{{km}}}) = 1$, то оператор ${{D}_{2}}[\rho ]$ аппроксимирует оператор Лапласа.

При практической реализации метода мы полагаем $f({{h}_{{km}}}) = h_{{km}}^{{3/2}}$. Тогда на равномерной декартовой сетке, имеющей шаг $h$, оператор ${{D}_{2}}[\rho ]$ аппроксимирует оператор ${{h}^{{3/2}}}\Delta $, а оператор ${{D}_{4}}\rho {\text{/}}{{V}_{k}}$ аппроксимирует оператор ${{\mu }_{4}}{{h}^{3}}\Delta \Delta $, который имеет порядок $O({{h}^{3}})$ и быстро убывает при измельчении сетки.

Для определения значений плотности на новом временном слое надо решить систему

правая часть которой определяется из уравнения (5). Чтобы эта система решалась эффективно, необходимо диагональное преобладание. Такое диагональное преобладание имеет место, когда коэффициент больше единицы. Этот коэффициент определяет скорость сходимости итерационного метода решения линейной системы относительно значений плотности на новом временном слое:

И учитывая уравнение (4), перепишем коэффициенты в виде

матрица ${{Z}_{k}}$ из уравнения (4), ${{\lambda }_{m}}$ являются произвольными числами, ${{I}_{k}}$ – количество узлов в шаблоне ${{S}_{A}}(k)$. Действительно,

${{\lambda }_{m}} = \sum\nolimits_{j = 1}^{J1} {{{{{\Lambda }}}_{{mj}}}} \int_{{{V}_{k}}}^{} {\psi _{j}^{k}(x,y,z)dV} $,

а элементы матрицы ${{\Lambda \;}}--$ произвольные числа. И всегда можно так выбрать элементы матрицы ${{\Lambda }}$, чтобы получить произвольно заданные значения ${{\lambda }_{m}}$.

Мы выбираем значения ${{\lambda }_{m}}$ так, чтобы улучшить диагональное преобладание (чтобы значение $\sigma $ было максимальным). Мы будем решать более простую задачу. Хорошими будут коэффициенты ${{B}_{{ki}}}$, если их часть будет равна нулю, исключая коэффициент ${{B}_{{kn(k)}}}$ на диагонали. Мы можем постараться сделать равными нулю $q = {{I}_{k}} - J1$ коэффициентов. Уравнения вида ${{B}_{{ki}}} = 0,~\;i = {{i}_{1}}, \ldots ,{{i}_{q}}$ ($1 \leqslant {{i}_{1}} < \ldots < {{i}_{q}} \leqslant {{I}_{k}},~\;{{i}_{s}} \ne n(k)$), задают линейную систему для определения множителей ${{\lambda }_{m}}$. Количество способов выбрать $q$ коэффициентов из ${{I}_{k}} - 1$ равно соответствующему числу сочетаний $C_{{Ik - 1}}^{q}$. Для каждой такой комбинации можно вычислить значение и определить лучшие значения ${{B}_{{ki}}}$. В случае декартовой трехмерной сетки число таких комбинаций $C_{{18}}^{9}$ (${{I}_{k}} = 19$, $J1$ = 10).

Так как эту задачу перебора надо решить для всех внутренних узлов сетки, такой подход может занимать слишком много времени. Еще сильнее упростим задачу нахождения коэффициентов ${{B}_{{ki}}}$, и для их поиска будем использовать методы решения задач линейного программирования. Запишем коэффициенты в виде ${{B}_{{ki}}} = B_{{ki}}^{*} + {{z}_{i}}$, где $B_{{ki}}^{*}$ было определено ранее, а вектор $z$ удовлетворяет уравнению ${{{{\Psi }}}^{{\text{T}}}}z = 0$. В дальнейшем для простоты будем использовать обозначение ${{b}_{i}} = {{B}_{{ki}}}$.

Сначала будем искать оптимальные значения коэффициентов ${{b}_{i}}$ такие, что ${{b}_{i}} \geqslant 0,~\;i = 1,2, \ldots ,{{I}_{k}}$. Коэффициенты ${{b}_{i}}$ удовлетворяют линейной системе ${{{{\Psi }}}^{{\text{T}}}}b = c$, где $c = {{{{\Psi }}}^{{\text{T}}}}{{B}_{k}}$. Решение будет в некотором смысле оптимальным, если знаменатель в выражении для диагонального преобладания $\sigma $ окажется минимальным, т.е. потребуем, чтобы сумма $\sum\nolimits_{i = 1,n(i) \ne k}^{Ik} {{{b}_{i}}} $ была минимальной. Задача

${{{{\Psi }}}^{{\text{T}}}}b = c$,   $\mathop {\min }\limits_b \left[ {\sum\nolimits_{i = 1,n(i) \ne k}^{Ik} {{{b}_{i}}} } \right]$,   ${{b}_{i}} \geqslant 0,$

является задачей линейного программирования и для ее решения можно использовать симплекс метод (см. [9]). Если минимум существует, мы считаем это решение оптимальным. Если минимум не достижим, мы считаем, что одно из значений ${{b}_{i}} < 0,~~n(i) \ne k$. Для каждого из этих случаев можно сформулировать и решить соответствующую задачу линейного программирования, а затем выбрать из найденных решений наилучшее. Опыт показывает, что таким образом достаточно приемлемое решение удавалось найти во всех рассматриваемых задачах.

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

Уравнения Навье–Стокса в декартовой системе координат в безразмерном виде имеют следующий вид:

$\rho $ – плотность, $u$, ${v}$, $w$ – компоненты скорости, $p$ – давление, $e$ – внутренняя энергия, $\mu $ – вязкость, ${\text{Re}}$, ${\text{Pr}} = 0.72$ – числа Рейнольдса и Прандтля.

Уравнения Навье–Стокса аппроксимируются по схеме, аналогичной схеме (5) для уравнения неразрывности

${{V}_{q}}$ – аппроксимация вязких членов, $q = 1,~2, \ldots ,5$ – номер уравнения.

При аппроксимации вязких членов необходимо уметь аппроксимировать вторые пространственные производные. Если в уравнении есть производная ${{\partial }^{2}}u{\text{/}}\partial x\partial y$, то после интегрирования по ячейке такое слагаемое будет аппроксимироваться следующим образом:

Такая аппроксимация содержит потоки через границы ячеек и не нарушает законов сохранения.

Рассмотрим аппроксимацию граничных условий. При постановке граничных условий или задается значение какой-либо физической величины в граничном узле, или граничное условие имеет вид $\frac{{\partial f}}{{\partial n}} = 0$, где $n$ – нормаль к границе области. Как правило, удается так расположить узлы вблизи границы расчетной области, что несколько узлов лежат на нормали к границе. Это дает возможность аппроксимировать производную по нормали с третьим порядком точности. При решении задач использовались следующие виды граничных условий.

1. На входе в область задаются полная энтальпия ${{h}_{0}} = \gamma e + {{u}^{2}}{\text{/}}2$, полное давление

где M – число Маха и ${v} = w = 0$. Для инварианта Римана ${{I}_{l}} = u - 2c{\text{/}}(\gamma - 1)$ задавалось условие $\frac{{\partial {{I}_{l}}}}{{\partial n}} = 0$.

2. На выходе из области задавалось давление и аппроксимировались условия

3. На стенке использовались условия $u = {v} = w = \partial p{\text{/}}\partial n = 0$, а значение внутренней энергии $e$ задавалось.

Отметим, что и в задаче об обтекании сферы, и в методической задаче об эволюции вихря использовались характерные величины: плотность ${{\rho }_{\infty }}$ и скорость звука ${{c}_{\infty }}$ в невозмущенном потоке (на бесконечности). Размерные величины $\rho {\kern 1pt} ',\;u{\kern 1pt} ',\;e{\kern 1pt} '$ связаны с безразмерными $\rho ,\;u,\;e$ соотношениями $\rho = \rho {\kern 1pt} {\text{'/}}{{\rho }_{\infty }}$, $u = u{\kern 1pt} {\text{'/}}{{c}_{\infty }}$, $e = e{\kern 1pt} {\text{'/}}c_{\infty }^{2}$.

2.7. Об аппроксимации вязких членов

При решении уравнений Навье–Стокса конвективные члены $\operatorname{conv} = \frac{{\partial {\mathbf{F}}}}{{\partial x}} + \frac{{\partial {\mathbf{G}}}}{{\partial y}} + \frac{{\partial {\mathbf{H}}}}{{\partial z}}$ аппроксимируются с точностью $O({{h}^{3}})$. А вязкие члены $\operatorname{vis} = - \left( {\frac{{\partial {{{\mathbf{F}}}_{V}}}}{{\partial x}} + \frac{{\partial {{{\mathbf{G}}}_{V}}}}{{\partial y}} + \frac{{\partial {{{\mathbf{H}}}_{V}}}}{{\partial z}}} \right)$ с точностью $O({{h}^{2}}){\text{/Re}}$. Это означает, что разность точных членов ${\text{conv}}$, ${\text{vis}}$ и их аппроксимаций $\operatorname{convh} $, $\operatorname{vish} $ удовлетворяет неравенствам

Рассмотрим случай $h > 1{\text{/Re}}$. Такое условие выполняется при решении задачи об обтекании сферы. Такое же условие обычно выполняется при численном решении других задач газовой динамики. В задачах об обтекании воздухом летательных аппаратов число Рейнольдса ${\text{Re}} = {{10}^{5}}{\kern 1pt} - {\kern 1pt} {{10}^{6}}$, и условие $h > 1{\text{/Re}}$ также выполняется, если не требуется экстремально повышенная точность.

В случае выполнения условия получаем

Если это условие не выполняется: $h < 1{\text{/Re}}~$ (это может потребоваться, если мы хотим получить решение точнее, чем $C{\text{/R}}{{{\text{e}}}^{3}}$). Тогда надо использовать другие методы, например, увеличивать шаблон для аппроксимации вязких членов.

3.1. Методические расчеты

Для проверки аппроксимационных свойств метода мы вычислим невязку (правую часть) на первом шаге по времени в задаче об эволюции вихря на расчетных сетках разной степени подробности. Областью, в которой будет решаться задача, является прямоугольный параллелепипед $0 \leqslant x \leqslant 10$, $0 \leqslant y \leqslant 10$, $0 \leqslant z \leqslant 25$. В расчетах мы используем сетки, узлы которых являются узлами равномерных декартовых сеток. Шаг по направлениям $x$ и $y$ равен $h$, а по направлению $z$ равен $10h$. В расчетах мы используем три сетки со следующим количеством узлов по трем координатам и соответствующим значением $h$: 1) $41 \times 41 \times 11$, $h = 0.25$; 2) $81 \times 81 \times 21$, $h = 0.125$; 3) $161 \times 161 \times 41$, $h = 0.0625$.

Прежде чем оценивать метод на разных сетках, проведем простой численный тест. Метод конечных объемов можно использовать для вычисления производных. Аппроксимируем выражение

$\int_{{{V}_{k}}}^{} {\partial u{\text{/}}\partial xdV} - \int_{{{S}_{k}}}^{} {u{{n}_{x}}dS} = 0$,

где интегрирование проводится по k-й ячейке и ее поверхности выражением

где ${{B}_{{kN}}}$ и ${{c}_{{kN}}}$ – коэффициенты операторов в уравнении (5), а $N = n(i)$. Первое слагаемое в этом уравнении примерно равно ${{V}_{k}}{{[\partial u{\text{/}}\partial x]}_{k}}$, ${{V}_{k}}$ – объем k-й ячейки. Поэтому выражение $\mathop {E(u) = \max }\limits_k \left| {{{Q}_{k}}} \right|{\text{/}}{{V}_{k}}$ определяет точность, с которой аппроксимируется производная $u$ по $x$. Мы вычисляли значение функции $E(u)$ на самой грубой расчетной сетке, которая была описана выше, максимум брался по всем внутренним узлам сетки, а функцией являлись полиномы от координат. Для функций $x,\;y,\;z$ значение $E(u)$ не превышало величины $1.22{{e}^{{ - 14}}}$, для квадратичных полиномов ${{x}^{2}},\;{{y}^{2}},\;{{z}^{2}},\;xy,\;xz,\;yz$ – величины $5.25{{e}^{{ - 13}}}$, для кубических полиномов ${{x}^{3}},\;{{y}^{3}},\;{{z}^{3}},~~\;{{x}^{2}}y,\;x{{y}^{2}},$ ${{x}^{2}}z,\;x{{z}^{2}},\;{{y}^{2}}z,\;y{{z}^{2}},\;~xyz$ – величины $2.07{{e}^{{ - 11}}}$. Похожие по порядку значения получались при дифференцировании тех же функций по $y$ и $z$. Таким образом, производные полиномов с порядком не выше третьего можно найти с точностью, не связанной с шагом сетки. Это подтверждает, что рассматриваемый метод имеет третий порядок аппроксимации по пространственным переменным.

В задаче об эволюции вихря газодинамические параметры связаны соотношениями

Это поле является точным решением уравнения Эйлера. В табл. 1 приведены значения нормы ${{\delta }_{{mq}}} = \sum\nolimits_k^{} {\left| {\delta {{f}_{{kq}}}} \right|{\text{/}}N} $ невязок для уравнений газовой динамики на первом шаге по времени, $m$ – номер сетки, $q$ – номер уравнения, максимум берется по всем внутренним узлам сетки, у которых 0.3 < z < 2.2. Мы рассматриваем только первый шаг, поскольку не удается поставить невозмущающих граничных условий на границах расчетной области, а постановка условий периодических по $z$ требует существенной доработки программы. В табл. 1 также приводится значение фактического порядка схемы, вычисляемого по формуле ${{\log }_{2}}{{\delta }_{{m - 1}}}{\text{/}}{{\delta }_{m}}$.

3.2. Задача о распространении гауссового импульса

Для оценки возможности метода, описывающего распространение акустических волн, рассмотрим тестовую задачу, предложенную в [10]. Решаются линеаризованные уравнения Эйлера, которые имеют вид

Задача рассматривается в области $ - D < x < D,$ $ - D < y < D$, $ - D < z < D$. В начальный момент времени $\rho = p = {\text{exp}}[ - {{r}^{2}}{\text{/}}2]$, ${{r}^{2}} = {{x}^{2}} + {{y}^{2}} + {{z}^{2}}$, $u = {v} = w = 0$. Точное аналитическое решение приведено в [11] и имеет вид

Рассматривается два варианта, для которых проводились расчеты на сетках разной степени подробности. В первом – размер области $~D = 10$, время расчета $t = 5$. Во втором $D = 20$, $t = 10$. Для аппроксимации по времени использовалась схема Рунге–Кутты четвертого порядка точности. При расчетах использовалось значение ${{\mu }_{4}} = 0.1$. Отличие численного давления от точного в нормах

${{s}_{{1k}}} = \sum\nolimits_i^{} {\left| {{{p}_{i}} - p_{i}^{{{\text{ex}}}}} \right|{\text{/}}L_{k}^{3}} $  и  ${{s}_{{2k}}} = {{\left( {\sum\nolimits_i^{} {{{{({{p}_{i}} - p_{i}^{{{\text{ex}}}})}}^{2}}} } \right)}^{{1/2}}}{\text{/}}L_{k}^{3}$

и соответствующий порядок при переходе на более подробную сетку

${{q}_{{1k}}} = \frac{{\log ({{s}_{{1k}}}{\text{/}}{{s}_{{1k + 1}}})}}{{\log ({{h}_{k}}{\text{/}}{{h}_{{k + 1}}})}}$  и  ${{q}_{{2k}}} = \frac{{\log ({{s}_{{2k}}}{\text{/}}{{s}_{{2k + 1}}})}}{{\log ({{h}_{k}}{\text{/}}{{h}_{{k + 1}}})}}$

приведены в табл. 2. Здесь L_k – количество узлов на ребре куба (расчетной области), h_k – шаг сетки. В некоторых случаях порядок меньше третьего. Но при переходе на самую подробную сетку порядок выше третьего достигается в обоих случаях.

3.3. Задача об обтекании пластины вязким газом

Рассмотрим задачу об обтекании плоской пластины вязким сжимаемым газом. Задача имеет простую геометрию. Известно приближенное аналитическое решение Блазиуса. Полученные численные решения можно использовать для оценки разных аппроксимаций вязких членов уравнений Навье–Стокса.

Расчетная область является прямоугольным параллелепипедом. Узлами используемой сетки являются узлы декартовой сетки с равными шагами по осям $x$ и $y$ и сгущением по оси $z$ к поверхности пластины. Ячейками являются ячейки Дирихле для этих узлов. Расчеты проводились на двух сетках, которые отличались количеством узлов в направлении $z.$ Количество узлов у грубой сетки $61 \times 21 \times 26$, минимальный шаг сетки по $z$ равен $0.01$. Количество узлов у подробной сетки $61 \times 21 \times 52$, минимальный шаг подробной сетки по $z$ равен $0.006$. Шаги в направлениях ${{h}_{x}} = {{h}_{y}} = 0.05$ у обеих сеток. Безразмерная длина области $3$ ($0 < x < 3$). Передняя кромка пластины расположена в точке с $x = 0.5$. Таким образом, точка с координатой $x = 1.5$ находится в середине расчетной области (по $x$), а соответствующая длина пластины равна единице. Для обезразмеривания плотности и скорости использовались плотность $\rho $ и скорость звука $c$ в набегающем потоке. Число Маха в набегающем потоке ${\text{M}} = 0.5$. Число Рейнольдса ${\text{Re}} = \frac{{\rho cl}}{\mu } = 250$.

При решении задачи использовались следующие граничные условия на гранях параллелепипеда:

– вход ($x = 0$): задаются параметры набегающего потока;

– выход ($x = 3$): скорость и внутренняя энергия экстраполируются, плотность находится по заданному выходному давлению;

– боковые грани: газодинамические параметры экстраполируются, ${v} = 0$;

– верхняя грань: горизонтальная компонента скорости экстраполируется, для плотности и внутренней энергии используются параметры набегающего потока;

– нижняя грань: поверхность пластины – скорость и нормальная производная давления равны нулю, задается температура пластины, равная температуре набегающего потока. Перед передней кромкой пластины: газодинамические параметры экстраполируются, $w = 0$.

На фиг. 2 приведено распределение коэффициента трения $~{{c}_{f}}$ вдоль линии равноудаленной от боковых краев пластины для численных решений, полученных на двух сетках (грубой и подробной):

На том же графике приведено значение коэффициента трения для решения Блазиуса, которое вычисляется по формуле Коэффициент трения полученных численных решений с высокой точностью согласуется с решением Блазиуса.

Используем полученные решения в качестве поля для оценки аппроксимаций вязких членов уравнений Навье–Стокса со вторым и четвертым порядками точности. Мы будем использовать уравнения Навье–Стокса, записанные в криволинейной системе координат. Узлы наших сеток (грубой и подробной) используются как узлы криволинейных сеток, на которых и будут вычисляться аппроксимации вязких членов. На фиг. 3 и 4 приведено сравнение аппроксимаций второго и четвертого порядков вязких членов в уравнении, содержащем $\left( {\frac{{\partial \rho u{\text{/}}J}}{{\partial t}}} \right)$, вдоль прямой, перпендикулярной пластине. Прямая проходит через точку на пластине, которая находится на расстоянии единицы (безразмерная длина) от ее передней кромки и на одинаковом расстоянии от боковых краев пластины. Отличия в вязких членах, полученных с разным порядком, малы на грубой сетке и практически незаметны на подробной.

3.4. Построение сетки в задаче об обтекании сферы

Построение сетки для данной задачи можно разбить на следующие три этапа: 1) построение квазиравномерной треугольной сетки на поверхности сферы, 2) перемещение точек по поверхности сферы со сгущением в области следа, 3) построение трехмерной сетки.

Для построения треугольной сетки на единичной сфере использовалась треугольная сетка, построенная на поверхности икосаэдра, вписанного в единичную сферу. Подобный подход описан в [12]. Вершины икосаэдра имеют следующие координаты: у точек ${{P}_{{ + / - }}}$ координаты $\left( {0,0, \pm 1} \right),$ координаты точек ${{Q}_{n}}$ – $\left( {r\cos \frac{{2\pi n}}{5},\;r\sin \frac{{2\pi n}}{5},\;h} \right)$, а точек ${{R}_{n}}$ – $\left( {r\cos \frac{{\pi (2n + 1)}}{5},\;r\sin \frac{{\pi (2n + 1)}}{5},\;h} \right)$, где $n = 1,2, \ldots ,5$. Учитывая, что длина всех ребер икосаэдра одинакова, можно записать $\left| {{{P}_{ + }}{{Q}_{1}}} \right| = \left| {{{Q}_{1}}{{Q}_{2}}} \right| = \left| {{{Q}_{1}}{{R}_{1}}} \right|$. Из этих двух уравнений, записанных в координатах, можно найти $r$ и $h$.

На каждой грани икосаэдра мы строим равномерную треугольную сетку. Если A, B, C – вершины некоторой грани икосаэдра, то радиус-вектор узлов этой сетки задается выражением ${{{\mathbf{r}}}_{p}} = {{{\mathbf{r}}}_{A}} + ({{{\mathbf{r}}}_{B}} - {{{\mathbf{r}}}_{A}})\left( {\frac{i}{m}} \right) + ({{{\mathbf{r}}}_{C}} - {{{\mathbf{r}}}_{A}})\left( {\frac{j}{m}} \right)$. Параметр $m$ – целое число, задающее подробность сетки, $i$ и $j$ – целые индексы, изменяющиеся от $0$ до $m$ и удовлетворяющие неравенству $i + j \leqslant m$.

Спроектируем эту сетку на поверхность единичной сферы. Точке P, являющейся узлом сетки на поверхности икосаэдра, соответствует точка S на поверхности сферы, так что ${{r}_{s}} = {{r}_{p}}{\text{/}}\left| {{{r}_{p}}} \right|$, ${{r}_{s}}$ и ${{r}_{p}}$ – радиус-векторы соответствующих точек. Квазиравномерная сетка на поверхности сферы построена.

Преобразуем эту сетку, чтобы она сгущалась вблизи точки с координатами $(1,0,0)$, т.е. трехмерная сетка будет сгущаться в области следа. Введем сферические координаты $\vartheta $, $\varphi $. Если точка на единичной сфере имеет координаты $x,\;y,\;z,$ то они будут связаны со сферическими координатами следующим образом: $\cos \varphi = y{\text{/}}{{r}_{{yz}}}$, $\sin \varphi = z{\text{/}}{{r}_{{yz}}}$, $\sin \vartheta = {{r}_{{yz}}} = \sqrt {{{y}^{2}} + {{z}^{2}}} $, $\cos \vartheta = - x.$ При преобразовании угол $\varphi $ не меняется, а зависимость нового угла $\theta $ от старого угла $\vartheta $ задается функцией $\theta (\vartheta ) = \sqrt {5{{\pi }^{2}} - {{{(\vartheta - 2\pi )}}^{2}}} - \pi $.

Для такой функции $\theta (0) = 0$, $\theta (\pi ) = \pi $, $\theta {\kern 1pt} '(0) = 2$, $\theta {\kern 1pt} '(\pi ) = 1{\text{/}}2$. Это означает, что сетка становится разреженной вблизи точки с координатами $( - 1,~0,~0)$ и сгущается вблизи точки $(1,~0,~0)$. Обе эти точки при преобразовании не перемещаются.

Трехмерная сетка включает узлы, расположенные на сферах с центром в начале координат. Координаты узлов на сфере с радиусом $r$ получаются умножением координат узлов на единичной сфере на этот радиус. Шаг, с которым увеличивается радиус сфер, экспоненциально увеличивается так, что величина $q = ({{r}_{{k + 1}}} - {{r}_{k}}){\text{/}}({{r}_{k}} - {{r}_{{k - 1}}})$ является постоянной. Проводя плоскости через середины отрезков, соединяющих близкие узлы сетки, мы получаем грани ячеек этих узлов. Ячейки нашей сетки являются пяти- или шестиугольными усеченными пирамидами. Прямые, проходящие через боковые ребра ячеек, пересекаются в начале координат.

Чтобы построить сетку по описанному выше методу, надо задать четыре параметра: количество точек на ребре икосаэдра, равное $m + 1$, минимальный шаг по радиусу ${{r}_{2}} - {{r}_{1}}$, коэффициент увеличения шага $q$, количество сферических слоев $l$.

3.5. Результаты численных расчетов в задаче об обтекании сферы

Задача об обтекании сферы жидкостью широко исследовалась во многих экспериментальных и численных работах (см. [13]–[19]). У задачи интересная геометрия, для которой нельзя построить одноблочную криволинейную (структурированную) сетку. В некоторых случаях задача решалась в сферической системе координат (см. [13], [14]) с расщеплением по времени, когда на втором шаге для давления решалось уравнение Пуассона. Для аппроксимации пространственных производных в [13] использовался метод коллокации на основе тригонометрических функций и полиномов Чебышёва, а в [14] использовались разностные аппроксимации второго порядка.

В [15] решалось уравнение для завихренности, а сама завихренность задавалась как линейная комбинация неких функций. В терминологии [15] лагранжевых частиц, характеризуемых местоположением (пространственными координатами) и весом. Такому полю завихренности соответствует определенное поле скорости, которое используется для вычисления местоположения частиц и их весов на следующем шаге по времени. В [16] задача об обтекании сферы решалась методом конечных объемов. Для интерполяции конвективных членов использовалась TVD-схема. Количество контрольных объемов в приведенных расчетах было 3.5–6 млн.

В случае метода погруженных границ (см. [17], [18]) данная задача использовалась в качестве тестовой. Узлы сетки в этом методе занимают всю область, включая обтекаемое тело. Однако в узлах, находящихся внутри тела или вблизи его поверхности, в уравнения добавляются корректирующие слагаемые.

Характер течения при обтекании сферы вязкой жидкостью существенно зависит от числа Рейнольдса, вычисленного по диаметру сферы $\operatorname{Re} = {{u}_{\infty }}d{\text{/}}\mu $. Перечислим некоторые диапазоны чисел Рейнольсда, приведенные в [16], при которых обтекание имеет определенный характер:

– $20 < {\text{Re}} < 212$ течение стационарное и осесимметричное;

– $212 < {\text{Re}} < 275$ течение стационарное, есть плоскость симметрии;

– $275 < {\text{Re}} < 375$ течение нестационарное, происходит регулярный сброс вихрей, есть плоскость симметрии;

– $375 < {\text{Re}} < 800$ азимутальный угол сброса вихрей меняется иррегулярно;

– $800 < {\text{Re}} < 3.7 \times {{10}^{5}}$ происходит отрыв ламинарного пограничного слоя.

Мы рассматриваем нестационарное обтекание сферы и приводим результаты расчетов для двух значений числа Рейнольдса ${\text{Re}} = 300$ и ${\text{Re}} = 350$. Результаты получены предлагаемым в статье методом. При решении линейной системы относительно значений в узлах ячеек на новом временном слое выполнялось 10 итераций. Значение константы ${{\mu }_{4}}$ равнялось 0.2 в расчете с ${\text{Re}} = 300$ и 0.3 в расчете с ${\text{Re}} = 350$. Для таких чисел Рейнольдса течение является нестационарным и характеризуется периодическим сходом вихрей с поверхности сферы.

Для расчета обтекания сферы с числом Рейнольдса ${\text{Re}} = 300$ использовалась сетка, которая характеризуется следующими параметрами. Количество точек на ребре икосаэдра 26, количество сферических слоев 51, минимальный шаг сетки по радиусу (у поверхности обтекаемой сферы) 0.02, коэффициент увеличения шага $q = 1.1$. Границы ячеек на поверхности сферы и граничные узлы этой сетки показаны на фиг. 5. Видно, что размеры ячеек в области натекания потока (в левой части рисунка) больше, чем в области следа (в правой части рисунка).

При рассмотрении подобных течений используют различные способы идентификации вихревых структур. Обычный подход, когда для идентификации вихрей используются значения модуля завихренности, плохо подходит для данной задачи из-за того, что модуль завихренности принимает большие значения в пограничном слое. Для определения области, где есть вихревые структуры, анализируют тензор скоростей деформации (см. [20]). В декартовой системе координат компоненты этого тензора – производные от компонент скорости $u,~\;{v},\;~w$ по $x,~\;y,\;~z$. Характер движения в окрестности точки определяется собственными значениями матрицы $A = \partial {{u}_{i}}{\text{/}}\partial {{x}_{j}}$, которые являются решениями уравнения ${{\lambda }^{3}} + P{{\lambda }^{2}} + Q\lambda + R = 0$, где величины $P,~\;Q,~\;R$ есть $P = - \operatorname{tr} [A]$, $Q = (\operatorname{tr} {{[A]}^{2}} - \operatorname{tr} [{{A}^{2}}]){\text{/}}2$, $R = - \det [A]$.

Если два собственных значения матрицы – комплексные сопряженные числа, то, как показано в [20], траектория в окрестности данной точки – это центр или фокус, что указывает на вихревой характер течения. В [16], где рассматривалось течение несжимаемой жидкости, и, значит, P $ = 0$, показано, что матрица имеет два комплексных собственных значения, если $D = {{(Q{\text{/}}3)}^{3}} + {{(R{\text{/}}2)}^{2}} > 0$. Наряду с критерием существования завихренности $D > 0$ в [13] используется более строгое условие $Q > 0$.

При рассматриваемом числе Рейнольдса ${\text{Re}} = 300$ течение нестационарное. В нашем расчете оно симметрично относительно плоскости XOZ и характеризуется квазипериодическим сходом вихрей с нижней части сферы. На фиг. 6 изображена изоповерхность $Q = {{10}^{{ - 5}}}$, соответствующая безразмерному времени $t = 500$. На фиг. 6 видно два последовательных вихря, сошедших с поверхности сферы.

На фиг. 7а показан график зависимости плотности от времени в точке за сферой с координатами (1, 0, –1/2). На фиг. 7б показан график спектральной плотности величины, показанной на левом рисунке, в зависимости от числа Струхаля ${\text{St}} = fd{\text{/}}{{u}_{\infty }}$. Второй максимум на этом графике соответствует частоте схода вихрей с поверхностей сферы. Это значение ${\text{St}} = 0.135$ хорошо соотносится с результатами других авторов, рассматривавших обтекание сферы несжимаемой жидкостью при числе Рейнольдса ${\text{Re}} = 300$: 0.137 из [14], 0.135 из [15], 0.133 из [16], 0.134 из [17], 0.134 из [18], 0.136 из [19].

Для расчета обтекания сферы с числом Рейнольдса ${\text{Re}} = 350$ использовалась более подробная сетка, которая характеризуется следующими параметрами: количество точек на ребре икосаэдра 34, количество сферических слоев 61, минимальный шаг сетки по радиусу (у поверхности обтекаемой сферы) 0.015, коэффициент увеличения шага $q = 1.075$.

Течение при этом числе Рейнольдса имеет более сложный характер и исследователи (см. [13]) отмечают две частоты в спектральной плоскости газодинамических параметров. На фиг. 8а приведен график зависимости плотности газа в точке за сферой с координатами (1, 0, 1/2). Видно, что примерно каждое третье колебание имеет большую амплитуду. На фиг. 8б приведена спектральная плотность приведенной на фиг. 8а величины. Максимумы достигаются в двух точках. Частота, соответствующая периоду схода вихрей, характеризуется числом Струхаля ${\text{St}} = 0.136$, что достаточно близко к значению $0.138$, приведенному в [13]. Значение второго числа Струхаля ${\text{St}} = 0.037$ согласуется со значением ${\text{St}} = 0.04$, которое было получено в [13]. Отметим также, что обтекание сферы при ${\text{Re}} = 350$ рассматривалось в [16], где максимумы соответствуют значениям ${\text{St}} = 0.0135$ и ${\text{St}} = 0.041$.