Improvement of the meshless method for numerical simulation of supersonic viscous gas flows

Full Text

Введение

Настоящая статья продолжает цикл работ авторов, посвященных численному исследованию газодинамического взаимодействия крупнодисперсных частиц с ударным слоем [1, 2]. Данное направление является важной составляющей решения актуальной задачи расчета обтекания затупленного тела высокоскоростным потоком газа с дисперсной примесью [3–6]. В работе развивается бессеточный метод численного интегрирования систем уравнений газовой динамики [2, 7–10].

В работе [1] показано, что интерференция ударных волн, вызванная выходом крупной частицы за пределы ударного слоя, приводит к возникновению локальных зон кратного усиления конвективного теплового потока от газа к поверхности тела [11–13], что согласуется с данными стендовых испытаний [14, 15]. Для детального исследования подобных явлений каждая моделируемая частица должна рассматриваться как макроскопический объект, вследствие чего задача приобретает комплексный и многомасштабный характер. При использовании традиционных методов решения систем уравнений газовой динамики возникает необходимость построения детальной расчетной сетки как у поверхности тела, так и в окрестности движущихся частиц.

Метод адаптивных скользящих декартовых сеток в сочетании с методом погруженной границы с фиктивными ячейками [16] позволяет решить задачу моделирования движения крупных частиц в ударном слое в плоской двумерной и осесимметричной постановках [1]. Для исследования коллективных эффектов газодинамического взаимодействия нескольких частиц, а также для расчета движения частиц вдоль сложных пространственных траекторий требуется решение системы уравнений газовой динамики в трехмерной постановке. Применение декартовых сеток с кубическими ячейками нецелесообразно ввиду высоких вычислительных затрат. Методы конечных элементов [17] и конечных объемов на неструктурированных сетках требуют либо модификации единой вычислительной сетки по мере перемещения частиц, либо использования движущихся сеток [18].

В работах настоящего цикла для решения трехмерных задач был выбран иной подход — бессеточный метод [7, 8], в котором в качестве альтернативы ячейкам вычислительной сетки используются распределенные в пространстве узлы. Основу метода составляет аппроксимация частных производных газодинамических величин и их функций по пространственным переменным через параметры газа в вычислительных узлах. Решение системы уравнений газовой динамики осуществляется в основной системе координат, связанной с обтекаемым телом, а также в локальных системах координат, перемещающихся вместе с частицами [9].

При практическом применении бессеточного метода к решению нестационарных пространственных задач вычислительной аэродинамики авторы выявили ряд особенностей, недостаточно освещенных в современной научной литературе и потребовавших дополнительной проработки. Одной из таких проблем, потребовавших решения, явилась асимметрия полей газодинамических величин, наблюдаемая при трехмерном расчете обтекания осесимметричного затупленного тела набегающим вдоль оси симметрии потоком газа. Проведенный анализ результатов расчетов показал, что данная асимметрия была обусловлена использованием единой системы координат для расчета потоков между узлами в системах уравнений Эйлера и Навье—Стокса. Введение для расчета конвективных потоков между парой вычислительных узлов локальной системы координат, ориентация которой определяется пространственной конфигурацией соседних узлов, позволяет существенно повысить точность решения задачи. Описание модифицированного бессеточного метода для решения системы уравнений Эйлера представлено авторами в работе [10]; в настоящей статье данный подход применяется для моделирования вязких течений путем решения системы уравнений Навье—Стокса.

В качестве модельной задачи рассматривалось обтекание сферы сверхзвуковым потоком вязкого газа. Наилучшие результаты получены при использовании комбинации полиномиальной (Polynomial Least Squares, PLS) [19] МНК-аппроксимации частных производных для конвективных потоков и аппроксимации рядом Тейлора (Taylor Least Squares, TLS) для вязких потоков [8, 20, 21] при дискретизации системы уравнений Навье—Стокса. Детальное описание разработанного варианта бессеточного метода приводится в настоящей статье.

1. Математическая модель обтекания затупленного тела потоком вязкого газа

Система нестационарных уравнений Навье—Стокса в консервативных переменных вместе с уравнением состояния идеального газа описывает течение вязкого газа в трехмерном пространстве в декартовой системе координат [22]:
\[\begin{gather*}
\frac{\partial \boldsymbol{q}}{\partial t} +\frac{\partial \boldsymbol{F}(\boldsymbol{q})}{\partial x} +\frac{\partial \boldsymbol{G}(\boldsymbol{q})}{\partial y} +\frac{\partial \boldsymbol{H}(\boldsymbol{q})}{\partial z} =\frac{\partial \boldsymbol{F}^{v}(\boldsymbol{q})}{\partial x} +\frac{\partial \boldsymbol{G}^{v}(\boldsymbol{q})}{\partial y} +\frac{\partial \boldsymbol{H}^{v} (\boldsymbol{q})}{\partial z} ,
\\
\boldsymbol{q}=\begin{pmatrix} {\rho } \\ {\rho u} \\ {\rho v} \\ {\rho w} \\ {\rho e} \end{pmatrix}, \quad
\boldsymbol{F}=\begin{pmatrix} {\rho u} \\ {\rho u^{2} +p} \\ {\rho uv} \\ {\rho uw} \\ {\rho uH} \end{pmatrix}, \quad
\boldsymbol{G}=\begin{pmatrix}{\rho v} \\ {\rho uv} \\ {\rho v^{2} +p} \\ {\rho vw} \\ {\rho vH} \end{pmatrix}, \quad
\boldsymbol{H}=\begin{pmatrix} {\rho w} \\ {\rho uw} \\ {\rho vw} \\ {\rho w^{2} +p} \\ {\rho wH} \end{pmatrix},
\\
\boldsymbol{F}^{v} =\begin{pmatrix} {0} \\ {\tau _{xx} } \\ {\tau _{xy} } \\ {\tau _{xz} } \\ {\tau _{xx} u+\tau _{xy} v+\tau _{xz} w-q_{x} } \end{pmatrix}, \quad
\boldsymbol{G}^{v} =\begin{pmatrix} {0} \\ {\tau _{yx} } \\ {\tau _{yy} } \\ {\tau _{yz} } \\ {\tau _{yx} u+\tau _{yy} v+\tau _{yz} w-q_{y} } \end{pmatrix}, 
\\
\boldsymbol{H}^{v} =\begin{pmatrix} {0} \\ {\tau _{zx} } \\ {\tau _{zy} } \\ {\tau _{zz} } \\ {\tau _{zx} u+\tau _{zy} v+\tau _{zz} w-q_{z} } \end{pmatrix},
\\
e=\frac{p}{\rho(\gamma-1)} +\frac{1}{2}(u^{2} +v^{2} +w^{2}), \quad H=e+\frac{p}{\rho} , \quad p=\rho RT,
\end{gather*}\]
где $t$ — время; $\rho $ — плотность; $p$ — давление; $T$ — температура; $u$, $v$, $w$ — проекции вектора скорости $\boldsymbol{v}$ на оси координат $x$, $y$ и $z$ соответственно; $\gamma $ — показатель адиабаты; $R$ — газовая постоянная; $H$ — полная энтальпия газа; $\boldsymbol{F}$, $\boldsymbol{G}$, $\boldsymbol{H}$ — векторы конвективных потоков, а $\boldsymbol{F}^{v} $, $\boldsymbol{G}^{v} $, $\boldsymbol{H}^{v} $ — векторы вязких потоков вдоль соответствующих координатных осей.

Компоненты тензора вязких напряжений определяются соотношениями:
\[\begin{gather*}
\tau _{xx} =\frac{2}{3} \mu \Bigl(2\frac{\partial u}{\partial x} -\frac{\partial v}{\partial y} -\frac{\partial w}{\partial z} \Bigr), \quad
\tau _{yy} =\frac{2}{3} \mu \Bigl(2\frac{\partial v}{\partial y} -\frac{\partial u}{\partial x} -\frac{\partial w}{\partial z} \Bigr),
\\
\tau _{zz} =\frac{2}{3} \mu \Bigl(2\frac{\partial w}{\partial z} -\frac{\partial u}{\partial x} -\frac{\partial v}{\partial y} \Bigr),
\quad
\tau _{xy} =\tau _{yx} =\mu \Bigl(\frac{\partial u}{\partial y} +\frac{\partial v}{\partial x} \Bigr), 
\\
\tau _{xz} =\tau _{zx} =\mu \Bigl(\frac{\partial u}{\partial z} +\frac{\partial w}{\partial x} \Bigr), \quad
\tau _{yz} =\tau _{zy} =\mu \Bigl(\frac{\partial v}{\partial z} +\frac{\partial w}{\partial y} \Bigr).
\end{gather*}\]
Плотность теплового потока определяется вектором с компонентами 
\[
q_{x} =-\lambda \frac{\partial T}{\partial x} ,\quad 
q_{y} =-\lambda \frac{\partial T}{\partial y} ,\quad 
q_{z} =-\lambda \frac{\partial T}{\partial z} .
\] 
Величина динамической вязкости вычисляется по формуле Сазерленда: 
\[
\mu =\mu ^{*} \Bigl(\frac{T}{T^{*} } \Bigr)^{{3}/{2} } \frac{T^{*} +C}{T+C} ,
\] 
где $\mu ^{*} =0.0000178$ (H${}\cdot{}$с)/м$^2$, $T^{*} =273.15$ K, $C =110.4$ K. 
Коэффициент теплопроводности пропорционален динамической вязкости: $\lambda {C_{p} \mu }/\textsf{Pr} $, где $C_{p}= {\gamma R}/({\gamma -1}) $ — удельная теплоемкость газа при постоянном давлении; число Прандтля для воздуха $\textsf{Pr} =0.72$.

Для дискретного представления поля газодинамических величин в расчетной области формируется конечное множество вычислительных узлов с фиксированными координатами. Положение узлов определяется геометрией области и особенностями течения; распределение узлов сгущается вблизи обтекаемых поверхностей в направлении внешней нормали для детального разрешения особенностей течения (рис. 1). В модельном примере обтекания сферы расчетные узлы располагаются на поверхностях концентрических сфер равномерно, так что на единицу площади сферы в среднем приходится равное количество узлов (рис. 2) [7].

 

Рис. 1. Распределение вычислительных узлов в центральном сечении расчетной области
[Figure 1. Distribution of computational nodes in the central cross-section of the computational domain]

 

Рис. 2. Распределение вычислительных узлов на поверхности сферы
[Figure 2. Distribution of computational nodes on the sphere surface]

 

Для каждого вычислительного узла $i$ определяется облако окружающих его соседних узлов $j\in C_{i} $ (рис. 3, 4). Соседние узлы выбираются в пределах заданного расстояния с фильтрацией узлов, имеющих близкое угловое расположение.

 

Рис. 3. Схема, иллюстрирующая вычислительный узел и его окрестность (облако соседних узлов)
[Figure 3. Schematic diagram of a computational node and its neighborhood (cloud of neighboring nodes)]

 

Рис. 4. Распределение облаков соседних узлов в расчетной области
[Figure 4. Distribution of neighbor node clouds in the computational domain]

 

2. Аппроксимация частных пространственных производных скалярной функции методом наименьших квадратов

При дискретизации векторов невязких потоков в системе уравнений Навье—Стокса используется метод наименьших квадратов с полиномиальной аппроксимацией [19], который определяет выражения для частных производных газодинамических параметров и содержащих их скалярных функций по пространственным координатам в виде линейной комбинации разностей значений функции в соседних узлах $j\in C_{i} $ и расчетном узле $i$. Вводится локальная нумерация узлов в облаке $C_{i}^{+} =\{i\}\cup C_{i} $, где узлу $i$ присваивается номер $0$, а его соседям — номера с $1$ по $n$; начало локальной системы координат помещается в узел $i$. Полиномиальное приближение первого порядка для скалярной функции $\varphi =\varphi(\boldsymbol{x})$, $\boldsymbol{x}=\begin{pmatrix} {x} & {y} & {z} \end{pmatrix}^{\top} $, с известными значениями в узлах облака $\varphi _{k} =\varphi(\boldsymbol{x}_{k})$, $k=0,1,\ldots, n$, имеет вид
\[
\hat{\varphi }(\boldsymbol{x})=\boldsymbol{p}^{\top}(\boldsymbol{x})\boldsymbol{a}, \quad 
\boldsymbol{p}^{\top}(\boldsymbol{x})=\begin{pmatrix} {1} & {x} & {y} & {z} \end{pmatrix}, \quad 
\boldsymbol{a}^{\top}=\begin{pmatrix} {a_{0} } &{a_{1} } & {a_{2} } & {a_{3} } \end{pmatrix}.
\] 

Вектор весовых коэффициентов $\boldsymbol{a}$ вычисляется методом наименьших квадратов по значениям в узлах $\hat{\varphi }(\boldsymbol{x}_{k})=\boldsymbol{p}^{\top}(\boldsymbol{x}_{k})\boldsymbol{a}$, $k=0,1,\ldots, n$, путем минимизации функционала
\[
f(\boldsymbol{a})=\sum _{k=0}^{n}\omega _{k}(\boldsymbol{p}^{\top}(\boldsymbol{x}_{k})\boldsymbol{a}-\varphi _{k})^{2}  \to \min ,
\] 
где в центральном узле $\omega _{0} =1$, а весовые коэффициенты соседних узлов обратно пропорциональны расстоянию между соседним и центральным узлами:
\[
\omega _{k} = {1}/{d_{0k}^{r} } , \quad 
d_{0k} =\sqrt{(x_{k} -x_{0})^{2} +(y_{k} -y_{0})^{2} +(z_{k} -z_{0})^{2} } , \quad 
k=1,\ldots, n.
\] 
В настоящей работе использовался показатель $r=3$ для весовых коэффициентов $\omega _{k} $ при аппроксимации конвективных потоков и $r^{v} =2$ для весовых коэффициентов $\omega _{ij}^{v} = {1}/{d_{ij}^{r^{v} } } $ при аппроксимации вязких потоков. Рекомендуемые в литературе значения $r$ лежат в пределах от $1$ до $6$ [23].

Согласно методу наименьших квадратов, вектор $\boldsymbol{a}$ определяется решением системы линейных алгебраических уравнений
\[
\boldsymbol{P}^{\top} \boldsymbol{W} \boldsymbol{P}\boldsymbol{a}=\boldsymbol{P}^{\top} \boldsymbol{W} \boldsymbol{\varphi} ,
\] 
где
\[\begin{gather*}
\boldsymbol{P}=\begin{pmatrix} {\boldsymbol{p}^{\top} (\boldsymbol{x}_{0})} \\ {\boldsymbol{p}^{\top}(\boldsymbol{x}_{1})} \\ {\vdots } \\ {\boldsymbol{p}^{\top}(\boldsymbol{x}_{n})} \end{pmatrix}
=
\begin{pmatrix} {1} & {x_{0} } & {y_{0} } & {z_{0} } \\ {1} & {x_{1} } & {y_{1} } & {z_{1} } \\ {\vdots } & {\vdots } & {\vdots } & {\vdots } \\ {1} & {x_{n} } & {y_{n} } & {z_{n} } \end{pmatrix}, 
\quad
\boldsymbol{\varphi} =\begin{pmatrix} {\varphi _{0} } \\ {\varphi _{1} } \\ {\vdots } \\ {\varphi _{n} } \end{pmatrix},
\\ 
\boldsymbol{W}=\operatorname{diag}\begin{pmatrix} {\omega _{0} } & {\omega _{1} } & {\ldots } & {\omega _{n} } \end{pmatrix}=\begin{pmatrix} {\omega _{0} } & {0} & {\cdots } & {0} \\ {0} & {\omega _{1} } & {\cdots } & {0} \\ {\vdots } & {\vdots } & {\ddots } & {\vdots } \\ {0} & {0} & {\cdots } & {\omega _{n} } \end{pmatrix}.
\end{gather*}\]

Искомый вектор коэффициентов:
\[
\boldsymbol{a}=(\boldsymbol{P}^{\top} \boldsymbol{W} \boldsymbol{P})^{-1} \boldsymbol{P}^{\top} \boldsymbol{W} \boldsymbol{\varphi} = \boldsymbol{C} \boldsymbol{\varphi} ,
\] 
где матрица $\boldsymbol{C}$ является постоянной и не зависит от $\varphi _{k} $:
\[
\boldsymbol{C}_{4\times (n+1)} =(\boldsymbol{P}^{\top} \boldsymbol{W} \boldsymbol{P})^{-1} \boldsymbol{P}^{\top} \boldsymbol{W}=\begin{pmatrix} {c_{00} } & {c_{01} } & {\cdots } & {c_{0n} } \\ {c_{10} } & {c_{11} } & {\cdots } & {c_{1n} } \\ {c_{20} } & {c_{21} } & {\cdots } & {c_{2n} } \\ {c_{30} } & {c_{31} } & {\cdots } & {c_{3n} } \end{pmatrix}.
\]  
Аппроксимация частных производных функции $\varphi (\boldsymbol{x})$ в узле $i(0)$:
\[\begin{gather*}
\frac{\partial \varphi }{\partial x} \approx \frac{\partial \hat{\varphi }}{\partial x} =a_{1} =\sum _{k=0}^{n}c_{1k} \varphi _{k}  , 
\\
\frac{\partial \varphi }{\partial y} \approx \frac{\partial \hat{\varphi }}{\partial y} =a_{2} =\sum _{k=0}^{n}c_{2k} \varphi _{k}  ,
\\ 
\frac{\partial \varphi }{\partial z} \approx \frac{\partial \hat{\varphi }}{\partial z} =a_{3} =\sum _{k=0}^{n}c_{3k} \varphi _{k}  .
\end{gather*} \]
В работе [24] показано, что $\displaystyle \sum _{k=0}^{n}c_{1k}  =\sum _{k=0}^{n}c_{2k}  =\sum  _{k=0}^{n}c_{3k}  =0$. Поэтому выражения для частных производных могут быть записаны в форме  
\[
\frac{\partial \hat{\varphi }}{\partial x} =\sum _{k=1}^{n}c_{1k}(\varphi _{k} -\varphi _{0}) , \quad 
\frac{\partial \hat{\varphi }}{\partial y} =\sum _{k=1}^{n}c_{2k}(\varphi _{k} -\varphi _{0}) , \quad
\frac{\partial \hat{\varphi }}{\partial z} =\sum _{k=1}^{n}c_{3k}(\varphi _{k} -\varphi _{0}) .
\] 
Возвращаясь к исходной нумерации узлов и обозначая $\alpha _{ij} =c_{1k} $, $\beta _{ij} =c_{2k} $, $\gamma _{ij} =c_{3k} $, получаем
\[
\frac{\partial \varphi }{\partial x} \Big|_{i} =\sum _{j\in C_{i} }\alpha _{ij}  \Delta \varphi _{ij} , \quad 
\frac{\partial \varphi }{\partial y} \Big|_{i} =\sum _{j\in C_{i} }\beta _{ij}  \Delta \varphi _{ij} , 
\] 
\[
\frac{\partial \varphi }{\partial z} \Big|_{i} =\sum _{j\in C_{i} }\gamma _{ij}  \Delta \varphi _{ij} ,  \quad 
\Delta \varphi _{ij} =\varphi _{j} -\varphi _{i} .
\] 

При дискретизации векторов вязких потоков применяется рассмотренный авторами в предыдущих работах [8, 9] метод аппроксимации, основанный на разложении функции $\varphi(\boldsymbol{x})$ в ряд Тейлора. В отличие от полиномиальной аппроксимации, он накладывает строгое требование на равенство в точке $\boldsymbol{x}_{0} $ значений исходной и приближенной функций: $\hat{\varphi }(\boldsymbol{x}_{0})=\varphi _{0} =\varphi(\boldsymbol{x}_{0})$. Для каждого расчетного узла $i$ и облака окружающих его соседних узлов $j\in C_{i} $ записываются приближенные равенства:
\[\begin{gather*}
\varphi _{j} =\varphi _{i} +\Delta x_{ij}\frac{\partial \varphi }{\partial x} \Big|_{i}^{v} +\Delta y_{ij} \frac{\partial \varphi }{\partial y} \Big|_{i}^{v} +\Delta z_{ij}\frac{\partial \varphi }{\partial z} \Big|_{i}^{v} +O(h^{2}),
\\
\Delta x_{ij} =x_{j} -x_{i} , \quad \Delta y_{ij} =y_{j} -y_{i} , \quad \Delta z_{ij} =z_{j} -z_{i} , \quad \Delta \varphi _{ij} =\varphi _{j} -\varphi _{i} .
\end{gather*}\]
Оптимальное приближение искомых частных производных $\dfrac{\partial \varphi }{\partial x} \Big|_{i}^{v} $, $\dfrac{\partial \varphi }{\partial y} \Big|_{i}^{v} $, $\dfrac{\partial \varphi }{\partial z} \Big|_{i}^{v} $ достигается минимизацией функционала
\[
\sum _{j\in C_{i} }\omega _{ij}^{v}  \Bigl(\Delta \varphi _{ij} -\Delta x_{ij} \frac{\partial \varphi }{\partial x} \Big|_{i}^{v} -\Delta y_{ij} \frac{\partial \varphi }{\partial y} \Big|_{i}^{v} -\Delta z_{ij} \frac{\partial \varphi }{\partial z} \Big|_{i}^{v} \Bigr)^{2} \to \min .
\] 
Значения коэффициентов линейной комбинации для представления пространственных производных функции $\varphi $ могут быть получены решением системы линейных уравнений
\[\begin{gather*}
\frac{\partial \varphi }{\partial x} \Big|_{i}^{v} =\sum _{j\in C_{i} }\alpha _{ij}^{v}  \Delta \varphi _{ij} , \quad 
\frac{\partial \varphi }{\partial y} \Big|_{i}^{v} =\sum _{j\in C_{i} }\beta _{ij}^{v}  \Delta \varphi _{ij} , \quad
\frac{\partial \varphi }{\partial z} \Big|_{i}^{v} =\sum _{j\in C_{i} }\gamma _{ij}^{v}  \Delta \varphi _{ij} ,
\\
\begin{bmatrix} 
{\alpha _{ij}^{v} } \rule{0mm}{16pt}% 
\\ 
{\beta _{ij}^{v} } \rule{0mm}{16pt}% 
\\ 
{\gamma _{ij}^{v} } \rule{0mm}{16pt}% 
\\
\end{bmatrix}=
\left[\!\!\!
\begin{array}{lll}
{\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta x_{ik}^{2}  } & {\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta x_{ik} \Delta y_{ik}  } & {\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta x_{ik} \Delta z_{ik}  } \\ 
{\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta x_{ik} \Delta y_{ik}  } & {\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta y_{ik}^{2}  } & {\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta y_{ik} \Delta z_{ik}  } \\ 
{\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta x_{ik} \Delta z_{ik}  } & {\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta y_{ik} \Delta z_{ik}  } & {\sum\limits _{k\in C_{i} }\omega _{ik}^{v} \Delta z_{ik}^{2}  } \end{array}\!\!\! 
\right]
^{-1} 
\!\!
\begin{bmatrix} 
{\omega _{ij}^{v} \Delta x_{ij} } \rule{0mm}{16pt}% 
\\ 
{\omega _{ij}^{v} \Delta y_{ij} } \rule{0mm}{16pt}% 
\\ 
{\omega _{ij}^{v} \Delta z_{ij} } \rule{0mm}{16pt}% 
\\[8pt]
\end{bmatrix}.
\end{gather*}\]

Коэффициенты линейной комбинации $\alpha _{ij} $, $\beta _{ij} $, $\gamma _{ij} $, а также $\alpha _{ij}^{v} $, $\beta _{ij}^{v} $, $\gamma _{ij}^{v} $ для расчета частных пространственных производных функции $\varphi $ определяются исключительно взаимным геометрическим расположением соседних узлов и являются константами для каждой пары узлов. 

3. Дискретизация системы уравнений газовой динамики и преобразование систем координат

При численном интегрировании системы уравнений Эйлера в роли функции $\varphi $ выступают газодинамические переменные (плотность, давление, температура), проекции вектора скорости газа на координатные оси, компоненты конвективных потоков и т. д.

Система уравнений газовой динамики с использованием полученной аппроксимации частных производных может быть записана в полудискретной форме:
\[\begin{multline*}
\frac{\partial \boldsymbol{q}_{i} }{\partial t} +2\sum _{j\in C_{i} }
\left[\alpha _{ij} (\boldsymbol{F}_{ij} -\boldsymbol{F}_{i} )+\beta _{ij} (\boldsymbol{G}_{ij} -\boldsymbol{G}_{i} )+\gamma _{ij} (\boldsymbol{H}_{ij} -\boldsymbol{H}_{i} )\right] =
\\
=2\sum _{j\in C_{i} }\left[\alpha _{ij}^{v} (\boldsymbol{F}_{ij}^{v} -\boldsymbol{F}_{i}^{v} )+\beta _{ij}^{v} (\boldsymbol{G}_{ij}^{v} -\boldsymbol{G}_{i}^{v} )+\gamma _{ij}^{v} (\boldsymbol{H}_{ij}^{v} -\boldsymbol{H}_{i}^{v} )\right] .
\end{multline*} \]
Здесь $\boldsymbol{F}_{ij} $, $\boldsymbol{G}_{ij} $, $\boldsymbol{H}_{ij} $ — векторы конвективных потоков, а $\boldsymbol{F}_{ij}^{v} $, $\boldsymbol{G}_{ij}^{v} $, $\boldsymbol{H}_{ij}^{v} $ — векторы вязких потоков в точке $M_{ij} $, являющейся серединой отрезка, соединяющего узлы $i$ и $j$ (рис. 4); $\boldsymbol{F}_{i} =\boldsymbol{F}(\boldsymbol{q}_{i})$; $\boldsymbol{G}_{i} =\boldsymbol{G}(\boldsymbol{q}_{i})$; $\boldsymbol{H}_{i} =\boldsymbol{H}(\boldsymbol{q}_{i})$, $\boldsymbol{F}_{i}^{v} =\boldsymbol{F}^{v} (\boldsymbol{q}_{i})$, $\boldsymbol{G}_{i}^{v} =\boldsymbol{G}^{v} (\boldsymbol{q}_{i})$, $\boldsymbol{H}_{i}^{v} =\boldsymbol{H}^{v} (\boldsymbol{q}_{i})$.

Поскольку в используемой для повышения порядка точности метода процедуре реконструкции векторов состояния газа ограничитель применяется отдельно к каждой компоненте вектора состояния, результат реконструкции, а в дальнейшем и решение задачи в целом зависят от выбора системы координат. В целях обеспечения инвариантности решения задачи относительно ориентации в пространстве облака узлов-соседей $j\in C_{i} $ узла $i$ вводится система координат $C'(0x'y'z')$, базис которой состоит из ортонормированных собственных векторов матрицы $A_{i} $ [23]:
\[
A_{i} =
\left[
\begin{array}{lll} 
{\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta x_{ij}^{2} } & {\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta x_{ij} \Delta y_{ij} } & {\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta x_{ij} \Delta z_{ij} } 
\\ 
{\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta x_{ij} \Delta y_{ij} } & {\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta y_{ij}^{2} } & {\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta y_{ij} \Delta z_{ij} } 
\\ 
{\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta x_{ij} \Delta z_{ij} } & {\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta y_{ij} \Delta z_{ij} } & {\sum\limits _{j\in C_{i} }\omega _{ij}  \Delta z_{ij}^{2} } 
\end{array}
\right]
,
\] 
где $\Delta x_{ij} =x_{j} -x_{i} $; $\Delta y_{ij} =y_{j} -y_{i} $; $\Delta z_{ij} =z_{j} -z_{i} $; $\omega _{ij} = {1}/{d_{ij}^{r} } $;  $d_{ij}  =\sqrt{\Delta x_{ij}^{2} +\Delta y_{ij}^{2} +\Delta z_{ij}^{2} } $. Матрица МНК $A_{i} $ является симметрической и положительно определенной, ее собственные значения положительны, а собственные векторы линейно независимы [23].

Авторами предлагается модификация данного подхода с адаптацией ориентации системы $C'(0x'y'z')$ для каждой пары соседних узлов. Ось $0x'$ определяется коэффициентами линейной комбинации для вычисления частных производных $\alpha _{ij} $, $\beta _{ij} $, $\gamma _{ij} $; направление $0y'$ — неколлинеарным оси $0x'$ собственным вектором $\boldsymbol{b}$ матрицы $A_{i} $ с наибольшим собственным значением, а ось $0z'$ дополняет систему координат $C'$ до правой:
\[\begin{gather*}
\boldsymbol{i}'=\begin{pmatrix} {x'_{x} } \\ {x'_{y} } \\ {x'_{z} } \end{pmatrix}=\frac{\boldsymbol{g}_{ij} }{\|\boldsymbol{g}_{ij}\|} , \quad 
\boldsymbol{g}_{ij} =\begin{pmatrix} {\alpha _{ij} } \\ {\beta _{ij} } \\ {\gamma _{ij} } \end{pmatrix}, 
\\
\boldsymbol{j}'=\begin{pmatrix} {y'_{x} } \\ {y'_{y} } \\ {y'_{z} } \end{pmatrix}=\frac{\boldsymbol{b}-(\boldsymbol{b}\cdot \boldsymbol{i}')\boldsymbol{i}'}{\|\boldsymbol{b}-(\boldsymbol{b}\cdot \boldsymbol{i}')\boldsymbol{i}'\|} ,
\quad 
\boldsymbol{b}\nparallel \boldsymbol{i}', \quad 
\boldsymbol{k}'=\begin{pmatrix} {z'_{x} } \\ {z'_{y} } \\ {z'_{z} } \end{pmatrix}=\boldsymbol{i}'\times \boldsymbol{j}'.
\end{gather*}\]
Как исходная система координат $C(0xyz)$, так и $C'(0x'y'z')$ являются декартовыми с ортонормированным базисом:
\[
\boldsymbol{j}'\bot \boldsymbol{i}',\quad \boldsymbol{k}'\bot \boldsymbol{i}',\quad \boldsymbol{k}'\bot \boldsymbol{j}';\]
\[
\boldsymbol{i}'=x'_{x} \boldsymbol{i}+x'_{y} \boldsymbol{j}+x'_{z} \boldsymbol{k},\quad 
\boldsymbol{j}'=y'_{x} \boldsymbol{i}+y'_{y} \boldsymbol{j}+y'_{z} \boldsymbol{k}, \quad 
\boldsymbol{k}'=z'_{x} \boldsymbol{i}+z'_{y} \boldsymbol{j}+z'_{z} \boldsymbol{k}.
\] 
Обратное преобразование орт системы координат:
\[
\boldsymbol{i}=x'_{x} \boldsymbol{i}'+y'_{x} \boldsymbol{j}'+z'_{x} \boldsymbol{k}', \quad 
\boldsymbol{j}=x'_{y} \boldsymbol{i}'+y'_{y} \boldsymbol{j}'+z'_{y} \boldsymbol{k}', \quad 
\boldsymbol{k}=x'_{z} \boldsymbol{i}'+y'_{z} \boldsymbol{j}'+z'_{z} \boldsymbol{k}'.
\] 
Прямое и обратное преобразование векторов при переходе между системами координат сводится к проецированию, то есть вычислению скалярного произведения вектора с ортом соответствующей оси:
\[\begin{gather*}
\boldsymbol{v}'=\begin{pmatrix} {u'} \\ {v'} \\ {w'} \end{pmatrix}=\begin{pmatrix} {\boldsymbol{v}\cdot \boldsymbol{i}'} \\ {\boldsymbol{v}\cdot \boldsymbol{j}'} \\ {\boldsymbol{v}\cdot \boldsymbol{k}'} \end{pmatrix}=\begin{pmatrix}{c} {ux'_{x} +vx'_{y} +wx'_{z} } \\ {uy'_{x} +vy'_{y} +wy'_{z} } \\ {uz'_{x} +vz'_{y} +wz'_{z} } \end{pmatrix},
\\
\boldsymbol{v}=\begin{pmatrix} {u} \\ {v} \\ {w} \end{pmatrix}=\begin{pmatrix} {\boldsymbol{v}'\cdot \boldsymbol{i}} \\ {\boldsymbol{v}'\cdot \boldsymbol{j}} \\ {\boldsymbol{v}'\cdot \boldsymbol{k}} \end{pmatrix}=\begin{pmatrix} {u'x'_{x} +v'y'_{x} +w'z'_{x} } \\ {u'x'_{y} +v'y'_{y} +w'z'_{y} } \\ {u'x'_{z} +v'y'_{z} +w'z'_{z} } \end{pmatrix}.
\end{gather*}\]
Скалярные величины при переходе между системами координат остаются неизменными, для них символ <<штрих>> далее опускается:
\[
\rho '=\rho , \quad p'=p, \quad e'=e, \quad T'=T.
\] 

Несмотря на то, что сами скалярные величины остаются неизменными при повороте системы координат, их градиент, используемый в процедуре MUSCL-реконструкции (Monotonic Upstream-centered Scheme for Conservation Laws) [21], требует пересчета. Данная процедура применяется для повышения порядка точности аппроксимации конвективных потоков. Рассмотрим вычисление градиента скалярной величины $\varphi \in\{\rho ,p,e,T\}$, компоненты которого являются частными производными по направлению, в системе координат $C'$ посредством преобразования вектора градиента этой величины, вычисленного в исходной системе координат $C$. Преобразование сводится к проецированию вектора градиента на оси новой системы координат:
\[\begin{gather*}
\frac{\partial \varphi }{\partial x'} \Big|_{i} =\sum _{k\in C_{i} }\alpha '_{ik}(\varphi _{k} -\varphi _{i}) =\boldsymbol{i}'\cdot \boldsymbol{\nabla} \varphi =x'_{x} \frac{\partial \varphi }{\partial x} \Big|_{i} +x'_{y} \frac{\partial \varphi }{\partial y} \Big|_{i} +x'_{z} \frac{\partial \varphi }{\partial z} \Big|_{i} ,
\\
\frac{\partial \varphi }{\partial y'} \Big|_{i} =\sum _{k\in C_{i} }\beta '_{ik} (\varphi _{k} -\varphi _{i})=\boldsymbol{j}'\cdot \boldsymbol{\nabla} \varphi =y'_{x} \frac{\partial \varphi }{\partial x} \Big|_{i} +y'_{y} \frac{\partial \varphi }{\partial y} \Big|_{i} +y'_{z} \frac{\partial \varphi }{\partial z} \Big|_{i} ,
\\
\frac{\partial \varphi }{\partial z'} \Big|_{i} =\sum _{k\in C_{i} }\gamma '_{ik} (\varphi _{k} -\varphi _{i})=\boldsymbol{k}'\cdot \boldsymbol{\nabla} \varphi =z'_{x} \frac{\partial \varphi }{\partial x} \Big|_{i} +z'_{y} \frac{\partial \varphi }{\partial y} \Big|_{i} +z'_{z} \frac{\partial \varphi }{\partial z} \Big|_{i} ;
\end{gather*} \]
вектор $\boldsymbol{g}'_{ik} $ определяет коэффициенты для вычисления градиента в $C'$:
\[
\boldsymbol{g}'_{ik} =
\begin{pmatrix} {\alpha '_{ik} } \\ {\beta '_{ik} } \\ {\gamma '_{ik} } \end{pmatrix}=\begin{pmatrix} {\boldsymbol{g}_{ik} \cdot \boldsymbol{i}'} \\ {\boldsymbol{g}_{ik} \cdot \boldsymbol{j}'} \\ {\boldsymbol{g}_{ik} \cdot \boldsymbol{k}'} \end{pmatrix}=\begin{pmatrix} {\alpha _{ik} x'_{x} +\beta _{ik} x'_{y} +\gamma _{ik} x'_{z} } \\ {\alpha _{ik} y'_{x} +\beta _{ik} y'_{y} +\gamma _{ik} y'_{z} } \\ {\alpha _{ik} z'_{x} +\beta _{ik} z'_{y} +\gamma _{ik} z'_{z} } 
\end{pmatrix}.
\] 
Рассмотрим вычисление градиента компонент вектора скорости $\boldsymbol{v}'$ в новом базисе путем преобразования градиента компонент скорости, вычисленных в старом базисе:
\[\begin{multline*}
\frac{\partial u'}{\partial x'} \Big|_{i} =\frac{\partial u'}{\partial x} \Big|_{i} \frac{\partial x}{\partial x'} \Big|_{i} +\frac{\partial u'}{\partial y} \Big|_{i} \frac{\partial y}{\partial x'} \Big|_{i} +\frac{\partial u'}{\partial z} \Big|_{i} \frac{\partial z}{\partial x'} \Big|_{i} = 
\\ =x'_{x} \frac{\partial u'}{\partial x} \Big|_{i} +x'_{y} \frac{\partial u'}{\partial y} \Big|_{i} +x'_{z} \frac{\partial u'}{\partial z} \Big|_{i} =
\\
=x'_{x} \Big(x'_{x} \frac{\partial u}{\partial x} \Big|_{i} +x'_{y} \frac{\partial v}{\partial x} \Big|_{i} +x'_{z} \frac{\partial w}{\partial x} \Big|_{i} \Big)+x'_{y} \Big(x'_{x} \frac{\partial u}{\partial y} \Big|_{i} +x'_{y} \frac{\partial v}{\partial y} \Big|_{i} +x'_{z} \frac{\partial w}{\partial y} \Big|_{i} \Big)+
\\
+x'_{z} \Big(x'_{x} \frac{\partial u}{\partial z} \Big|_{i} +x'_{y} \frac{\partial v}{\partial z} \Big|_{i} +x'_{z} \frac{\partial w}{\partial z} \Big|_{i} \Big)=\boldsymbol{i}'\cdot 
\begin{pmatrix} {\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } & {\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } & {\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } \end{pmatrix}^{\top} ,
\end{multline*} \]
\[\begin{multline*}
\frac{\partial u'}{\partial y'} \Big|_{i} =\frac{\partial u'}{\partial x} \Big|_{i} \frac{\partial x}{\partial y'} \Big|_{i} +\frac{\partial u'}{\partial y} \Big|_{i} \frac{\partial y}{\partial y'} \Big|_{i} +\frac{\partial u'}{\partial z} \Big|_{i} \frac{\partial z}{\partial y'} \Big|_{i} = \\ 
=y'_{x} \frac{\partial u'}{\partial x} \Big|_{i} +y'_{y} \frac{\partial u'}{\partial y} \Big|_{i} +y'_{z} \frac{\partial u'}{\partial z} \Big|_{i} =
\\
=y'_{x} \Big(x'_{x} \frac{\partial u}{\partial x} \Big|_{i} +x'_{y} \frac{\partial v}{\partial x} \Big|_{i} +x'_{z} \frac{\partial w}{\partial x} \Big|_{i} \Big)+y'_{y} \Big(x'_{x} \frac{\partial u}{\partial y} \Big|_{i} +x'_{y} \frac{\partial v}{\partial y} \Big|_{i} +x'_{z} \frac{\partial w}{\partial y} \Big|_{i} \Big)+
\\
+y'_{z} \Big(x'_{x} \frac{\partial u}{\partial z} \Big|_{i} +x'_{y} \frac{\partial v}{\partial z} \Big|_{i} +x'_{z} \frac{\partial w}{\partial z} \Big|_{i} \Big)=\boldsymbol{j}'\cdot 
\begin{pmatrix} {\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } & {\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } & {\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } \end{pmatrix}^{\top} ,
\end{multline*}\]
и т. д. Здесь введены следующий обозначения:  
\[\begin{gather*}
\frac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} =
\begin{pmatrix} {\dfrac{\partial u}{\partial x} \Big|_{i} } & { \dfrac{\partial v}{\partial x} \Big|_{i} } & {\dfrac{\partial w}{\partial x} \Big|_{i} } \end{pmatrix}^{\top}, 
\quad
\frac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} =
\begin{pmatrix} {\dfrac{\partial u}{\partial y} \Big|_{i} } & {\dfrac{\partial v}{\partial y} \Big|_{i} } & {\dfrac{\partial w}{\partial y} \Big|_{i} } \end{pmatrix}^{\top}, 
\\
\frac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} =
\begin{pmatrix} {\dfrac{\partial u}{\partial z} \Big|_{i} } & {\dfrac{\partial v}{\partial z} \Big|_{i} } & {\dfrac{\partial w}{\partial z} \Big|_{i} } \end{pmatrix}^{\top}.
\end{gather*} \]
В результате получим
\[\begin{gather*}
\frac{\partial u'}{\partial x'} \Big|_{i} =\boldsymbol{i}'\cdot 
\begin{pmatrix} 
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\quad 
\frac{\partial u'}{\partial y'} \Big|_{i} =\boldsymbol{j}'\cdot 
\begin{pmatrix} 
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\quad 
\frac{\partial u'}{\partial z'} \Big|_{i} =\boldsymbol{k}'\cdot 
\begin{pmatrix} 
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\ 
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{i}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\\
\frac{\partial v'}{\partial x'} \Big|_{i} =\boldsymbol{i}'\cdot 
\begin{pmatrix} 
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\quad 
\frac{\partial v'}{\partial y'} \Big|_{i} =\boldsymbol{j}'\cdot 
\begin{pmatrix} 
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\ 
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\quad 
\frac{\partial v'}{\partial z'} \Big|_{i} =\boldsymbol{k}'\cdot 
\begin{pmatrix} 
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{j}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\\
\frac{\partial w'}{\partial x'} \Big|_{i} =\boldsymbol{i}'\cdot 
\begin{pmatrix} 
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \big|_{i} } \\
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\ 
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\;\;
\frac{\partial w'}{\partial y'} \Big|_{i} =\boldsymbol{j}'\cdot 
\begin{pmatrix} 
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix},
\; \;
\frac{\partial w'}{\partial z'} \Big|_{i} =\boldsymbol{k}'\cdot 
\begin{pmatrix} 
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial x} \Big|_{i} } \\
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial y} \Big|_{i} } \\
{\boldsymbol{k}'\cdot \dfrac{\partial \boldsymbol{v}}{\partial z} \Big|_{i} } 
\end{pmatrix}.
\end{gather*} \]

Преобразование градиентов позволяет однократно на каждом шаге расчета вычислить градиент скалярной величины или компоненты вектора в узле $i$ в основной системе координат $C$ и получить его компоненты в $C'$ для реконструкции векторов состояния и расчета векторов конвективных потоков для каждой пары, которую узел $i$ образует со своими соседями.

4. Вычисление векторов конвективных потоков

Покомпонентная процедура MUSCL-реконструкции  с ограничителем Ван Албады 2 (van Albada 2 limiter) [21] применяется к вектору состояния в физических переменных $\boldsymbol{\Psi} '=\begin{pmatrix} \rho & u' & v' & w'& p \end{pmatrix}^{\top} $ в системе координат $C'$:
\[\begin{gather*}
{\psi '}_{ij}^{+} =\psi '_{i} +\frac{s_{i} }{4}\bigl[(1-ks_{i})\Delta _{ij}^{-} +(1+ks_{i})(\psi '_{j} -\psi '_{i} )\bigr], \;
\Delta _{ij}^{-} =2\boldsymbol{r}'_{ij} \cdot \boldsymbol{\nabla }\psi '_{i} -(\psi '_{j} -\psi '_{i} ),
\\
{\psi '}_{ij}^{-} =\psi '_{j} -\frac{s_{j} }{4}\bigl[(1-ks_{j})\Delta _{ij}^{+} +(1+ks_{j})(\psi '_{j} -\psi '_{i})\bigr], \; 
\Delta _{ij}^{+} =2\boldsymbol{r}'_{ij} \cdot \boldsymbol{\nabla }\psi '_{j} -(\psi '_{j} -\psi '_{i} ),
\\ 
s_{i} =\max \Bigl(0,\frac{2\Delta _{ij}^{-} (\psi '_{j} -\psi '_{i})+\varepsilon }{\Delta _{ij}^{-2} +(\psi '_{j} -\psi '_{i} )^{2} +\varepsilon } \Bigr), 
\quad
s_{j} =\max \Bigl(0, \frac{2\Delta _{ij}^{+} (\psi '_{j} -\psi '_{i} )+\varepsilon }{\Delta _{ij}^{+2} +(\psi '_{j} -\psi '_{i} )^{2} +\varepsilon } \Bigr),
\\
\boldsymbol{r}'_{ij} =
\begin{pmatrix} {x'_{j} -x'_{i} } & {y'_{j} -y'_{i} } & {z'_{j} -z'_{i} } \end{pmatrix}^{\top}, 
\quad 
\boldsymbol{\nabla} \psi '_{n} =
\begin{pmatrix} {\dfrac{\partial \psi '}{\partial x'} \Big|_{n} } & 
{\dfrac{\partial \psi '}{\partial y'} \Big|_{n} } & {\dfrac{\partial \psi '}{\partial z'} \Big|_{n} } \end{pmatrix}^{\top},
\end{gather*}\]
где коэффициент $k$ определяет порядок аппроксимации схемы: при $k= {1}/{3} $ — третий порядок, при $k=0$ или $k=\pm 1$ — второй. 
В настоящей работе применяются $k={1}/{3} $, $\varepsilon =10^{-13} $. Здесь следует обратить внимание на необходимость вычисления градиента $\boldsymbol{\nabla} \psi '_{n} $ в системе координат $C'$.

Производится восстановление векторов консервативных переменных ${\boldsymbol{q}'}^{-} =\boldsymbol{q}({\boldsymbol{\Psi} '}_{ij}^{-} )$ и ${\boldsymbol{q}'}^{+} =\boldsymbol{q}({\boldsymbol{\Psi} '}_{ij}^{+} )$ в системе координат $C'$.

Для расчета потоков на границе разрыва используются усредненные по методу Роу (Roe approximate Riemann solver) физические величины, которые обеспечивают выполнение свойств сохранения и позволяют избежать появления нефизических решений:
\[\begin{gather*}
\hat{\rho }=\sqrt{\rho ^{-} \rho ^{+} } , \; 
\hat{H}=\frac{H^{-} \sqrt{\rho ^{-} } +H^{+} \sqrt{\rho ^{+} } }{\sqrt{\rho ^{-} } +\sqrt{\rho ^{+} } } , \; 
\hat{c}'=\sqrt{\left\{\hat{H}-\frac{\hat{u}'^{2} +\hat{v}'^{2} +\hat{w}'^{2} }{2} \right\}(\gamma -1)} ,
\\
\hat{u}'=\frac{u'^{-} \sqrt{\rho ^{-} } +u'^{+} \sqrt{\rho ^{+} } }{\sqrt{\rho ^{-} } +\sqrt{\rho ^{+} } } , \; 
\hat{v}'=\frac{v'^{-} \sqrt{\rho ^{-} } +v'^{+} \sqrt{\rho ^{+} } }{\sqrt{\rho ^{-} } +\sqrt{\rho ^{+} } } , \;
\hat{w}'=\frac{w'^{-} \sqrt{\rho ^{-} } +w'^{+} \sqrt{\rho ^{+} } }{\sqrt{\rho ^{-} } +\sqrt{\rho ^{+} } } ,
\\
\hat{\boldsymbol{v}}'=\begin{pmatrix} {\hat{u}'} & {\hat{v}'} & {\hat{w}'} \end{pmatrix}^{\top} .
\end{gather*} \]

Схема расчета векторов конвективных потоков SLAU (Simple Low-dissipation AUSM) [25], относящаяся к семейству схем AUSM (Advection Upstream Splitting Method), адаптирована для бессеточного метода по аналогии со схемой HLLC (Harten–Lax–van Leer Contact)  в работе [26] и применяется к вектору $\boldsymbol{K}$, представляющему собой линейную комбинацию векторов потоков $\boldsymbol{F}$, $\boldsymbol{G}$ и $\boldsymbol{H}$ вдоль координатных осей, а не к каждому из них по отдельности:
\[
\boldsymbol{K}=\eta _{x} \boldsymbol{F}+\eta _{y} \boldsymbol{G}+\eta _{z} \boldsymbol{H}=\begin{pmatrix} {\rho \vartheta } & {\rho u\vartheta +p\eta _{x} } & {\rho v\vartheta +p\eta _{y} } & {\rho w\vartheta +p\eta _{z} } & {\rho \vartheta H} \end{pmatrix}^{\top}, 
\]
\[
\boldsymbol{\eta} =\frac{\boldsymbol{g}_{ij} }{\lambda _{ij} } , \quad
\lambda _{ij} =\|\boldsymbol{g}_{ij}\|=\sqrt{\alpha _{ij}^{2} +\beta _{ij}^{2} +\gamma _{ij}^{2} } ,
\] 
где $\vartheta =\boldsymbol{v}\cdot \boldsymbol{\eta} $ — проекция вектора скорости газа $\boldsymbol{v}$ на направление $\boldsymbol{\eta} $.

Расчет конвективного потока $\boldsymbol{K}'_{ij} $ в системе координат $C'$ в точке $M_{ij} $ в соответствии со схемой SLAU производится следующим образом:
\[\begin{gather*}
\boldsymbol{q}'_{L} =\boldsymbol{q}'^{+} , \quad \boldsymbol{q}'_{R} =\boldsymbol{q}'^{-} , \quad \boldsymbol{\eta} '=\frac{\boldsymbol{g}'_{ij} }{\lambda '_{ij} } , \quad \lambda '_{ij} =\|\boldsymbol{g}'_{ij} \|=\|\boldsymbol{g}_{ij} \|=\lambda _{ij} ,
\\
\boldsymbol{K}'_{ij} =\frac{m_{ij}^{*} +|m_{ij}^{*}|}{2} \boldsymbol{h}'_{L} +\frac{m_{ij}^{*} -|m_{ij}^{*}|}{2} \boldsymbol{h}'_{R} +p_{ij}^{*} \boldsymbol{N}', 
\\
\boldsymbol{h}'=\begin{pmatrix} {0} & {u'} & {v'} & {w'} & {e+\frac{p}{\rho } } \end{pmatrix}^{\top}, 
\quad 
\boldsymbol{N}'=\begin{pmatrix} {0} & {\eta '_{x} } & {\eta '_{y} } & {\eta '_{z} } & {0} \end{pmatrix}^{\top},
\\
\vartheta _{L} =\boldsymbol{v}'_{L} \cdot \boldsymbol{\eta} ', \quad 
\vartheta _{R} =\boldsymbol{v}'_{R} \cdot \boldsymbol{\eta} ', \quad
\bar{c}=\frac{c_{L} +c_{R} }{2} , \quad
M_{L} =\frac{\vartheta _{L} }{\bar{c}} , \quad
M_{R} =\frac{\vartheta _{R} }{\bar{c}} ,
\\
g=-\max(\min(M_{L} ,0),-1)\cdot \min(\max(M_{R} ,0),1)\in[{0,1}],
\\
\tilde{M}=\min \Bigl(1,\frac{1}{\bar{c}}  \sqrt{\frac{\|\boldsymbol{v}'_{L}\|^{2} +\|\boldsymbol{v}'_{R}\|^{2} }{2} } \Bigr), \quad
f_{p} =(1-\tilde{M})^{2} , \quad
|\bar{\vartheta }|=\frac{\rho _{L}|\vartheta _{L}|+\rho _{R}|\vartheta _{R}|}{\rho _{L} +\rho _{R} } ,
\\
m_{ij}^{*} =\frac{1-g}{2} \cdot \bigl(\rho _{L} \vartheta _{L} +\rho _{R} \vartheta _{R} -|\bar{\vartheta }|(\rho _{R} -\rho _{L} ) \bigr)-f_{p} \cdot \frac{p_{R} -p_{L} }{2\bar{c}} ,
\\
\tilde{p}_{L} \big|_{\alpha = {3}/{16} } =\begin{cases}
\frac{1}{4}(M_{L} +1)^{2}(2-M_{L})+\alpha M_{L}(M_{L}^{2} -1)^{2} , & |M_{L}|\le 1 ;
\\ 
\frac{1}{2} (1+\operatorname{sign}(M_{L} )), & |M_{L}|>1 ,
\end{cases}
\\
\tilde{p}_{R} \big|_{\alpha = {3}/{16} } =\begin{cases}
\frac{1}{4}(M_{R} -1)^{2}(2+M_{R})-\alpha M_{R}(M_{R}^{2} -1)^{2} , & |M_{R}|\le 1 ;
\\ 
\frac{1}{2}(1-\operatorname{sign}(M_{R})), & |M_{R}|>1 ,
\end{cases}
\\
p_{ij}^{*} =\frac{\tilde{p}_{L} +\tilde{p}_{R} }{2} \cdot (p_{L} -p_{R})+ 
\bigl(1+(1-f_{p})(\tilde{p}_{L} +\tilde{p}_{R} -1) \bigr)\cdot \frac{p_{L} +p_{R} }{2} .
\end{gather*}\]

Вектор конвективного потока $\boldsymbol{K}_{ij} $ вычисляется посредством обратного преобразования физических компонент вектора $\boldsymbol{K}'_{ij} $ из системы координат $C'$ в исходную систему координат $C$:
\[
\rho =\rho ', \quad e=e', \quad 
\boldsymbol{v}=\begin{pmatrix} {u} & {v} & {w} \end{pmatrix}^{\top}=
\begin{pmatrix} {\boldsymbol{v}'\cdot \boldsymbol{i}} & {\boldsymbol{v}'\cdot \boldsymbol{j}} & {\boldsymbol{v}'\cdot \boldsymbol{k}} \end{pmatrix}^{\top}.
\]

5. Вычисление векторов вязких потоков

Расчет векторов вязких потоков ввиду линейности и отсутствия ограничителей в процедуре реконструкции производится в системе координат $C$, преобразование координат при этом не требуется. Частные пространственные производные температуры и компонент скорости
\[\begin{gather*}
\frac{\partial u}{\partial x} \Big|_{i}^{v} =\sum _{j\in C_{i} }\alpha _{ij}^{v}  (u_{j} -u_{i} ), \; 
\frac{\partial u}{\partial y} \Big|_{i}^{v} =\sum _{j\in C_{i} }\beta _{ij}^{v}  (u_{j} -u_{i}), \; \dots  , \; 
\frac{\partial T}{\partial z} \Big|_{i}^{v} =\sum _{j\in C_{i} }\gamma _{ij}^{v} (T_{j} -T_{i}),
\\
\tau _{xx} \big|_{i} =\frac{2}{3} \mu _{i} \Bigl(2\frac{\partial u}{\partial x} \Big|_{i}^{v} -\frac{\partial v}{\partial y} \Big|_{i}^{v} -\frac{\partial w}{\partial z} \Big|_{i}^{v} \Bigr), \; \dots  , \; 
\tau _{xy} \big|_{i} =\mu _{i} \Bigl(\frac{\partial u}{\partial y} \Big|_{i}^{v} +\frac{\partial v}{\partial x} \Big|_{i}^{v} \Bigr), \; \dots  , 
\\
 q_{z} \big|_{i} =-\lambda _{i} \frac{\partial T}{\partial z} \Big|_{i}^{v} 
\end{gather*}\]
непосредственно используются для расчета векторов вязких потоков $\boldsymbol{F}_{i}^{v}  =\boldsymbol{F}^{v} (\boldsymbol{q}_{i})$, $\boldsymbol{G}_{i}^{v} =\boldsymbol{G}^{v} (\boldsymbol{q}_{i})$, $\boldsymbol{H}_{i}^{v} =\boldsymbol{H}^{v} (\boldsymbol{q}_{i})$ в узле $i$ согласно соотношениям, приведенным в первом разделе настоящей статьи.

Реконструкция векторов градиента физических переменных $u$, $v$, $w$, $T$, необходимых для вычисления вязких потоков $\boldsymbol{F}_{ij}^{v} $, $\boldsymbol{G}_{ij}^{v} $, $\boldsymbol{H}_{ij}^{v} $ в середине отрезка, соединяющего узлы $i$ и $j$, выполняется согласно [27]:
\[
\boldsymbol{\nabla }\varphi \big|_{ij}^{v} =
\overline{\boldsymbol{\nabla }\varphi \big|_{ij}^{v} }-
\Bigl(
\overline{\boldsymbol{\nabla }\varphi \big|_{ij}^{v} }\cdot \frac{\boldsymbol{r}_{ij} }{\|\boldsymbol{r}_{ij}\|} -\frac{\varphi _{j} -\varphi _{i} }{\|\boldsymbol{r}_{ij}\|} 
\Bigr)
\frac{\boldsymbol{r}_{ij} }{\|\boldsymbol{r}_{ij}\|} , \quad 
\overline{\boldsymbol{\nabla }\varphi \big|_{ij}^{v} }=\frac{\boldsymbol{\nabla }\varphi \big|_{i}^{v} +\boldsymbol{\nabla }\varphi \big|_{j}^{v} }{2} .
\] 
Например,
\[\begin{multline*}
\frac{\partial u}{\partial y} \Big|_{ij}^{v} =
\overline{\frac{\partial u}{\partial y} \Big|_{ij}^{v} }-
\Bigl(\overline{\frac{\partial u}{\partial x} \Big|_{ij}^{v} }(x_{j} -x_{i})+\overline{\frac{\partial u}{\partial y} \Big|_{ij}^{v} }(y_{j} -y_{i}) +
\\
 +\overline{\frac{\partial u}{\partial z} \Big|_{ij}^{v} }(z_{j} -z_{i})
-\frac{u_{j} -u_{i} }{\|\boldsymbol{r}_{ij}\|} \Bigr)\frac{y_{j} -y_{i} }{\|\boldsymbol{r}_{ij}\|} ,
\end{multline*} \]
\[
\overline{\frac{\partial u}{\partial x} \Big|_{ij}^{v} }=\frac{1}{2} \Bigl(\frac{\partial u}{\partial x} \Big|_{i}^{v} +\frac{\partial u}{\partial x} \Big|_{j}^{v} \Bigr),
\quad 
\overline{\frac{\partial u}{\partial y} \Big|_{ij}^{v} }=\frac{1}{2} \Bigl(\frac{\partial u}{\partial y} \Big|_{i}^{v} +\frac{\partial u}{\partial y} \Big|_{j}^{v} \Bigr),
\quad 
\overline{\frac{\partial u}{\partial z} \Big|_{ij}^{v} }=\frac{1}{2} \Bigl(\frac{\partial u}{\partial z} \Big|_{i}^{v} +\frac{\partial u}{\partial z} \Big|_{j}^{v} \Bigr),
\]
\[\begin{multline*}
\frac{\partial v}{\partial x} \Big|_{ij}^{v} =
\overline{\frac{\partial v}{\partial x} \Big|_{ij}^{v} }-
\Bigl(\overline{\frac{\partial v}{\partial x} \Big|_{ij}^{v} }(x_{j} -x_{i})+\overline{\frac{\partial v}{\partial y} \Big|_{ij}^{v} }(y_{j} -y_{i})+
\\
+\overline{\frac{\partial v}{\partial z} \Big|_{ij}^{v} }(z_{j} -z_{i})-\frac{v_{j} -v_{i} }{\|\boldsymbol{r}_{ij}\|} \Bigr) \frac{x_{j} -x_{i} }{\|\boldsymbol{r}_{ij}\|} ,
\end{multline*} \]
\[
\overline{\frac{\partial v}{\partial x} \Big|_{ij}^{v} }=\frac{1}{2} \Bigl(\frac{\partial v}{\partial x} \Big|_{i}^{v} +\frac{\partial v}{\partial x} \Big|_{j}^{v} \Bigr),
\quad
\overline{\frac{\partial v}{\partial y} \Big|_{ij}^{v} }=\frac{1}{2} \Bigl(\frac{\partial v}{\partial y} \Big|_{i}^{v} +\frac{\partial v}{\partial y} \Big|_{j}^{v} \Bigr),
\quad 
\overline{\frac{\partial v}{\partial z} \Big|_{ij}^{v} }=\frac{1}{2} \Bigl(\frac{\partial v}{\partial z} \Big|_{i}^{v} +\frac{\partial v}{\partial z} \Big|_{j}^{v} \Bigr),
\]
вязкость определяется усреднением значений в узлах $\mu _{ij} =\dfrac{\mu _{i} +\mu _{j} }{2} $, тогда $\tau _{xy} \big|_{ij} =\mu _{ij} \Bigl(\dfrac{\partial u}{\partial y} \Big|_{ij}^{v} +\dfrac{\partial v}{\partial x} \Big|_{ij}^{v} \Bigr)$, остальные компоненты тензора вязких напряжений вычисляются аналогичным образом.

6. Численное интегрирование системы уравнений газовой динамики

Численное интегрирование системы уравнений газовой динамики по времени выполняется явным методом Рунге—Кутты третьего порядка точности [28]:
\[\begin{gather*}
\frac{\partial \boldsymbol{q}_{i} }{\partial t} +R(\boldsymbol{q}_{i})=R^{v}(\boldsymbol{q}_{i}), \quad 
R(\boldsymbol{q}_{i})=2\sum _{j\in C_{i} }\lambda _{ij}[\boldsymbol{K}_{ij} -\boldsymbol{K}(\boldsymbol{q}_{i})] ,
\\
R^{v}(\boldsymbol{q}_{i})=2\sum _{j\in C_{i} }\left[\alpha _{ij}^{v}(\boldsymbol{F}_{ij}^{v} -\boldsymbol{F}_{i}^{v})+\beta _{ij}^{v}(\boldsymbol{G}_{ij}^{v} -\boldsymbol{G}_{i}^{v} )+\gamma _{ij}^{v} (\boldsymbol{H}_{ij}^{v} -\boldsymbol{H}_{i}^{v})\right] ,
\\
\boldsymbol{q}_{i}^{(1)} =\boldsymbol{q}_{i}^{n} -\Delta tR(\boldsymbol{q}_{i}^{n}), \quad 
\boldsymbol{q}_{i}^{(2)} =\frac{3}{4} \boldsymbol{q}_{i}^{n} +\frac{1}{4} \boldsymbol{q}_{i}^{(1)} -\frac{1}{4} \Delta tR(\boldsymbol{q}_{i}^{(1)}),
\\
\boldsymbol{q}_{i}^{n+1} =\frac{1}{3} \boldsymbol{q}_{i}^{n} +\frac{2}{3} \boldsymbol{q}_{i}^{(2)} -\frac{2}{3} \Delta tR(\boldsymbol{q}_{i}^{(2)}).
\end{gather*} \]

Величина временного шага интегрирования определяется критерием Куранта по скоростям, полученным с помощью усредненных по Роу (Roe-averaging) физических величин, которые обеспечивают корректное определение характеристических скоростей на границе разрыва:
\[\begin{gather*}
\Delta t=\min \Bigl(\frac{\mathrm{CFL}}{\max\limits_{i}|\nu _{i}|} ,\frac{\mathrm{VNN}}{\max\limits_{i}|\nu _{i}^{v}|}\Bigr),
\end{gather*}\]
где $\mathrm{CFL}=0.5$ — число Куранта–Фридрихса–Леви (Courant–Friedrichs–Lewy condition) для конвективных членов, $\mathrm{VNN}=0.4$ — параметр устойчивости для вязких членов (Viscous Neumann Number),
\[\begin{gather*}
\nu _{i} =\sum _{j\in C_{i} } \bigl(\alpha '_{ij} \hat{u}'_{ij} +\beta '_{ij} \hat{v}'_{ij} +\gamma '_{ij} \hat{w}'_{ij} +\lambda '_{ij} \hat{c}'_{ij} \bigr) ,
\\
\nu _{i}^{v} =\max \Bigl(\frac{4}{3} ,\frac{\gamma ^{ {3}/{2} } }{\mathsf{Pr} } \Bigr)
\sum _{j\in C_{i} }\frac{2\bar{\mu }_{ij} ({\alpha _{ij}^{v}}^{2} +{\beta _{ij}^{v}}^{2} +{\gamma _{ij}^{v}}^{2} )}{\rho _{i} +\rho _{j} }  ,
\quad 
\bar{\mu }_{ij} =\mu \Bigl(\frac{T_{i} +T_{j} }{2} \Bigr).
\end{gather*}\]

7. Вычислительный эксперимент

Программная реализация представленных выше моделей и алгоритмов выполнена авторами статьи на языке программирования C++ в сочетании с технологией гетерогенных параллельных вычислений OpenCL. Программа предназначена для работы под управлением операционных систем семейств Linux и Windows и для ускорения вычислений позволяет использовать графические процессоры различных производителей.

В качестве модельной решается задача сверхзвукового обтекания сферы радиуса $R_{s} =55$ мм воздухом с показателем адиабаты $\gamma =1.4$. На входе в расчетную область задаются граничные условия первого рода, определяющие невозмущенный сверхзвуковой поток с фиксированной температурой $T_{\infty } =300$ K, числами Маха $\mathsf{M}_{\infty } =3$ и Рейнольдса $\mathsf{Re}_{\infty } =  {\rho _{\infty } u_{\infty } L}/{\mu } =10^{5} $, где $L$ — характерный размер преграды (диаметр сферы); на выходной границе — условия Неймана $\frac{\partial \boldsymbol{q}}{\partial \boldsymbol{n}} =0$, где $\boldsymbol{n}$ — вектор внешней нормали к границе. Поверхность сферы рассматривается как твердая стенка с известной температурой $T_{{w}} =500$ K, условиями прилипания и непротекания газа $\boldsymbol{v}=0$. На каждом вычислительном шаге параметры газа в узлах, расположенных на поверхности тела, рассчитываются согласно общему алгоритму для внутренних узлов расчетной области, вычитается вектор скорости $\rho u=\rho v=\rho w=0$, $\tilde{e}=e-\frac{1}{2} \|\boldsymbol{v}\|^{2} $, затем выполняется корректировка плотности для выравнивания температуры газа с температурой стенки $\tilde{\rho }=\frac{\rho \tilde{e}(\gamma -1)}{RT_{{w}} } $. Количество узлов в вычислительных экспериментах при решении задачи в трехмерной постановке варьировалось в пределах от $N_{1} \approx 2\cdot 10^{5} $ до $5^{3} \cdot N_{1} \approx 2.5\cdot 10^{7} $.

На рис. 5 показана теневая картина течения газа [29]; штриховой кривой отмечено приближенно-аналитическое положение головной ударной волны [30]. Представлены поля давления, плотности и числа Маха в сечении области расчета полуплоскостью, проходящей через ось симметрии (рис. 6–8).

 

Рис. 5. Теневая картина сверхзвукового обтекания сферы
[Figure 5. Schlieren image of supersonic flow past a sphere]

 

Рис. 6. Распределение плотности газа в центральном сечении области расчета при сверхзвуковом обтекании сферы
[Figure 6. Gas density distribution in the central cross-section of the computational domain for supersonic flow past a sphere]

 

Рис. 7. Распределение давления газа в центральном сечении области расчета при сверхзвуковом обтекании сферы
[Figure 7. Gas pressure distribution in the central cross-section of the computational domain for supersonic flow past a sphere]

 

Рис. 8. Распределение числа Маха в центральном сечении области расчета при сверхзвукового обтекании сферы
[Figure 8. Mach number distribution in the central cross-section of the computational domain for supersonic flow past a sphere]

 

Даже при относительно небольшом количестве вычислительных узлов наблюдается хорошее согласование давления газа на поверхности сферы, полученного расчетом бессеточным методом, с эталонными данными [31] (рис. 9).

 

Рис. 9. Сравнение расчетного распределения давления на поверхности сферы с эталонными данными
[Figure 9. Comparison of the calculated surface pressure distribution on the sphere with benchmark data]

 

На рис. 10 показан график конвективного теплового потока от газа к поверхности тела $Q=-\lambda _{w} \frac{\partial T}{\partial \boldsymbol{n}_{w} } $, где $\lambda _{w} =\lambda (T_{w})$ — коэффициент теплопроводности газа у поверхности тела; $\boldsymbol{n}_{w} $ — вектор внешней нормали к границе; частная производная температуры газа по нормали к поверхности тела в вычислительном узле $i$ на поверхности тела 
\[
\frac{\partial T}{\partial \boldsymbol{n}_{w} } \Big|_{i} =\sum _{j\in C_{i} }[(n_{wx} \alpha _{ij}^{v} +n_{wy} \beta _{ij}^{v} +n_{wz} \gamma _{ij}^{v})(T_{j} -T_{i})] .
\]
В качестве эталонной используется приближенно-аналитическая кривая зависимости конвективного теплового потока от угла отклонения от критической точки $\alpha $: $Q(\alpha)=Q_{FR} (0.55+0.45\cos 2\alpha)$ [32]. Значение $Q_{FR} $ в критической точке рассчитывается по формуле Фэя и Ридделла [33].

 

Рис. 10. Сравнение расчетного конвективного теплового потока к поверхности сферы с приближенно-аналитическим решением
[Figure 10. Comparison of the calculated convective heat flux to the sphere surface with an approximate analytical solution]

 

По мере увеличения количества вычислительных узлов отклонение расчетной кривой конвективного теплового потока от приближенно-аналитической сокращается и в расчете с максимальным числом узлов не превышает 2%. Поведение представленных графиков в зависимости от количества узлов свидетельствует о сходимости рассматриваемого численного метода.

Заключение

В результате серии работ по совершенствованию бессеточного метода численного решения уравнений газовой динамики [21, 27, 34] разработан комплекс моделей и алгоритмов для расчета сверхзвукового обтекания в трехмерной постановке. Отличительными чертами настоящей работы, позволившими повысить точность относительно предыдущих реализаций метода [2, 7–10] и снизить погрешности, обусловленные асимметрией в распределении вычислительных узлов [10], являются:

1. адаптивное преобразование системы координат для каждой пары узлов при расчете конвективных потоков;
2. применение схемы расчета векторов конвективных потоков к их линейной комбинации вдоль координатных осей;
3. использование различных типов аппроксимации пространственных производных при расчете векторов конвективных и вязких потоков.

Предложенный модифицированный бессеточный метод апробирован на задаче сверхзвукового обтекания сферы. Показано, что даже при относительно небольшом количестве вычислительных узлов ($\approx 2\cdot 10^{5}$) достигается хорошее согласование с эталонными данными по распределению давления на поверхности тела. При увеличении количества узлов до $2.5\cdot 10^{7}$ отклонение расчетной кривой конвективного теплового потока от приближенно-аналитической не превышает 2%, что свидетельствует о сходимости метода.

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

Конкурирующие интересы. У нас нет конфликта интересов в отношении авторства и публикации этой статьи.
Авторский вклад и ответственность. Все  авторы  принимали  участие  в  разработке  концепции статьи и в написании рукописи. Авторы несут  полную  ответственность  за  предоставление  окончательной рукописи в печать. Окончательная  версия рукописи была одобрена всеми авторами.
Финансирование. Работа выполнена при поддержке Министерства науки и высшего образования Российской Федерации (тема № FSFF-2023-0008) в рамках государственного задания  Московского авиационного института (национального исследовательского университета).

About the authors

Dmitry L. Reviznikov

Moscow Aviation Institute (National Research University)

Author for correspondence.
Email: reviznikov@inbox.ru
ORCID iD: 0000-0003-0998-7975
SPIN-code: 9763-9898
Scopus Author ID: 6602701797
ResearcherId: T-4571-2018
https://www.mathnet.ru/rus/person26170

Dr. Phys. & Math. Sci., Professor; Professor; Dept. of Computational Mathematics and Programming

Russian Federation, 125993, Moscow, Volokolamskoe Shosse, 4

Andrey V. Sposobin

Moscow Aviation Institute (National Research University)

Email: spise@inbox.ru
ORCID iD: 0009-0004-7720-2556
SPIN-code: 3028-0859
Scopus Author ID: 42462268200
ResearcherId: AAV-1002-2020
https://www.mathnet.ru/rus/person34344

Dr. Phys. & Math. Sci.; Senior Researcher; R&D Dept. 806

Russian Federation, 125993, Moscow, Volokolamskoe Shosse, 4

References

Supplementary files