Research of methods for calculating gradient of gas-dynamic quantities in finite volume difference schemes in aerodynamic problems

Full Text

Введение

Точность численного моделирования напрямую зависит от точности выполнения различных операций, одной из которых является вычисление градиента. Значение, полученное в результате этой операции, применяется, например, при построении схем интерполяции второго порядка [1-2]. Градиенты чаще всего вычисляются либо методом наименьших квадратов (МНК), либо методом Грина-Гаусса. При работе с тонкоячеистыми сетками МНК может давать значения градиентов с большой погрешностью, что в конечном итоге снижает устойчивость всего численного решения [1]. В этом контексте метод Грина-Гаусса проявляет себя как более надежный инструмент, обеспечивающий повышенную точность вычислений на сетках данного типа. Если же сеточная модель содержит ячейки с неортогональными гранями, наоборот – наблюдается более высокая точность МНК.

Известно, что неструктурированные сетки боле востребованы на практике [3-6] ввиду более легкого процесса построения средствами сеточных генераторов и удобства заполнения расчетной области (тонкие вытянутые ячейки в призматическом слое, в основном ядре сетки – крупные многогранники). В настоящей работе разработан авторский гибридный алгоритм вычисления градиентов, предназначенный для повышения точности расчетов на неструктурированных сетках. Метод интегрирует преимущества двух классических подходов ‒ метода Грина-Гаусса и метода наименьших квадратов (МНК). Ключевыми аспектами разработки являются адаптивная комбинация методов в зависимости от локальных характеристик сетки и оптимизированный алгоритм выбора весовых коэффициентов. Практическая значимость предлагаемого решения подтверждена тестовыми расчетами и реализацией в программном комплексе.

Основные уравнения

Система уравнений Навье-Стокса [6-8], используемая для описания течения жидкости и газа, записывается следующим образом:

$\left\{\begin{array}{l}\frac{\partial \text{ρ}}{\partial t}+\nabla \cdot \left(\text{ρ}\overrightarrow{u}\right)=\mathrm{0,}\\ \frac{\partial \left(\text{ρ}\overrightarrow{u}\right)}{\partial t}+\nabla \cdot \left(\text{ρ}\overrightarrow{u}\overrightarrow{u}\right)=-\nabla p+\nabla \cdot \left({\text{τ}}_{\text{μ}}+{\text{τ}}_t\right),\\ \frac{\partial \left(\text{ρ}E\right)}{\partial t}+\nabla \cdot \left(\text{ρ}\overrightarrow{u}h\right)=\nabla \cdot \left[\overrightarrow{u}\left({\text{τ}}_{\text{μ}}+{\text{τ}}_t\right)-\left({\overrightarrow{q}}_{\text{μ}}+{\overrightarrow{q}}_t\right)\right].\end{array}\right.$ (1)

В выражении (1): ρ – плотность среды; $\overrightarrow{u}$ – вектор скорости (u – скорость по оси x, v – скорость по оси y, w – скорость по оси z); p – давление; $E=C_{v}T+0.5\left(u^2+v^2+w^2\right)$ – полная энергия газа на единицу массы; $h=C_{p}T+0.5\left(u^2+v^2+w^2\right)$ – полная энтальпия газа; ${\text{τ}}_{\text{μ}}$ – молекулярная составляющая тензора касательных напряжений; ${\text{τ}}_t$ – турбулентная составляющая тензора касательных напряжений; $q_{\text{μ}}$ – молекулярная составляющая плотности теплового потока; $q_t$ – турбулентная составляющая плотности теплового потока; T – температура; $C_v=\left(C_{p}T-R/m\right)$ – удельная теплоемкость при постоянном объеме; $C_p$ – удельная теплоемкость при постоянном давлении; R – универсальная газовая постоянная; m – молярная масса газа.

Для решения системы уравнений (1) используется конечно-объемный метод, в рамках которого расчетная область дискретизируется на произвольные многогранные элементы, выступающие в роли контрольных объемов (ячеек сетки – рис. 1).

 

Рис. 1. Ячейка расчетной сетки

Fig. 1. Computational grid cell

 

Система (1) в векторной форме записывается так:

$\frac{d}{dt}\underset{\Delta V}{\int }WdV+\underset{\Delta {\sum }_P}{\oint }\left(F-G\right)dS=\underset{\Delta V}{\int }HdV,$ (2)

где W – вектор консервативных переменных, F – вектор конвективных потоков, G – вектор диффузионных потоков, H – слагаемое источников.

$W=\left(\begin{array}{c}\rho \\ \rho u\\ \rho v\\ \rho w\\ \rho E\end{array}\right),F=\left(\begin{array}{c}\rho u_n\\ \rho u{u}_n+p{n}_x\\ \rho v{u}_n+p{n}_y\\ \rho w{u}_n+p{n}_z\\ \rho H{u}_n+p{u}_n\end{array}\right),G=\left(\begin{array}{c}0\\ {\tau }_{nx}\\ {\tau }_{ny}\\ {\tau }_{nz}\\ \tau u+q\end{array}\right),$ (3)

где u_n – скорость по нормали, q – тепловой поток, τ_nj – произведение тензора вязких напряжений и вектора нормали. Для расчета вектора конвективных потоков применяются соответствующие численные схемы, например, AUSM+ [9].

Расчет конвективных и диффузионных потоков подразумевает определение величин на гранях расчетных ячеек, для чего применяется их реконструкция, которая в случае второго порядка точности [9-10] записывается в следующем виде:

$\begin{array}{l}{\varphi }_f^{-}={\varphi }_P^{}+{\alpha }_f^{-}(\Delta {\overrightarrow{R}}_{Pf}\cdot \nabla {\varphi }_P^{}),\\ {\varphi }_f^{+}={\varphi }_E^{}+{\alpha }_f^{+}(\Delta {\overrightarrow{R}}_{Ef}\cdot \nabla {\varphi }_E^{}),\\ \Delta {\overrightarrow{R}}_{Pf}={\overrightarrow{R}}_f-{\overrightarrow{R}}_P=\left(x_f-x_P\right)\overrightarrow{i}+\left(y_f-y_P\right)\overrightarrow{j}+\left(z_f-z_P\right)\overrightarrow{k}=\Delta x_f\overrightarrow{i}+\Delta y_f\overrightarrow{j}+\Delta z_f\overrightarrow{k}\\ \Delta {\overrightarrow{R}}_{Ef}={\overrightarrow{R}}_f-{\overrightarrow{R}}_E=\left(x_f-x_E\right)\overrightarrow{i}+\left(y_f-y_E\right)\overrightarrow{j}+\left(z_f-z_E\right)\overrightarrow{k}=\Delta x_f\overrightarrow{i}+\Delta y_f\overrightarrow{j}+\Delta z_f\overrightarrow{k}\end{array}$ (4)

где ${\varphi }_f^{-}$ и ${\varphi }_f^{+}$ – значение слева и справа от грани, ${\varphi }_P$ и ${\varphi }_E$ – значение в центре ячейки E и P (рис. 2), $\Delta {\overrightarrow{R}}_{Pf}$ и $\Delta {\overrightarrow{R}}_{Ef}$ – расстояние от центра ячейки Е и Р до центра грани f, (x_i, y_i, z_i – декартовые координаты, $\overrightarrow{i},\overrightarrow{j,}\overrightarrow{k}$ – единичные вектора в декартовой системе координат), $\nabla {\varphi }_{Е}$ и $\nabla {\varphi }_{Р}$ – значение градиента в ячейке E и P, ${\alpha }_f^{-}$ и ${\alpha }_f^{+}$ – ограничитель, применяемый для снижения осцилляций в области ударных волн и контактных разрывов [11].

 

Рис. 2. Схема реконструкции величин из центра ячеек на грань

Fig. 2. Scheme of reconstruction of values from the center of cells to the edge

 

Выражение для вычисления градиента по методу Грина-Гаусса [12]:

$\nabla {\phi }_P^{}=\frac{1}{V_P}\sum _{f=face\left(P\right)}{\phi }_f{S}_f,$ (5)

где ${\phi }_f$ – значение на грани, $S_f$ – площадь грани.

Значение на грани f между ячейками P и E может быть вычислено следующим образом.

1. ${\phi }_f=\frac{1}{2}\left({\phi }_P^{}+{\phi }_{Е}^{}\right)$ (6)
2. ${\phi }_f=\lambda {\phi }_P^{}+\left(1-\lambda \right){\phi }_E^{}$ (7)

где $\lambda =\frac{\left|{\overrightarrow{R}}_{EP}\right|}{\left|\Delta {\overrightarrow{R}}_{Pf}\right|+\left|\Delta {\overrightarrow{R}}_{Ef}\right|}$, $\Delta {\overrightarrow{R}}_{EP}$ – вектор, соединяющий центры ячеек P и E, 

$\Delta {\overrightarrow{R}}_{EP}={\overrightarrow{R}}_P-{\overrightarrow{R}}_E=\left(x_P-x_E\right)\overrightarrow{i}+\left(y_P-y_E\right)\overrightarrow{j}+\left(z_P-z_E\right)\overrightarrow{k}=\Delta x_f\overrightarrow{i}+\Delta y_f\overrightarrow{j}+\Delta z_f\overrightarrow{k}$.

Также градиент может быть рассчитан по методу наименьших квадратов (МНК) [12]. В этом случае вектор градиента записывается в виде:

${\left(grad\phi \right)}_E=\left\{{\left(\frac{\partial \phi }{\partial x}\right)}_E,{\left(\frac{\partial \phi }{\partial y}\right)}_E,{\left(\frac{\partial \phi }{\partial z}\right)}_E\right\}={\left(A,B,C\right)}_E$.

(A, B, C) рассчитываются на основе решения приведенной системы:

$\left\{\begin{array}{l}A\sum _{f=1}^F{\left(\Delta x_f\right)}^2+B\sum _{f=1}^F\Delta x_f\Delta y_f+C\sum _{f=1}^F\Delta x_f\Delta z_f=\sum _{f=1}^F\Delta x_f\Delta {\phi }_f,\\ A\sum _{f=1}^F\Delta x_f\Delta y_f+B\sum _{f=1}^F{\left(\Delta y_f\right)}^2+C\sum _{f=1}^F\Delta y_f\Delta z_f=\sum _{f=1}^F\Delta y_f\Delta {\phi }_f,\\ A\sum _{f=1}^F\Delta x_f\Delta z_f+B\sum _{f=1}^F\Delta y_f\Delta z_f+C\sum _{f=1}^F{\left(\Delta z_f\right)}^2=\sum _{f=1}^F\Delta z_f\Delta {\phi }_f.\end{array}\right.$ (10)

Система (10) может быть записана как:

${\left(grad\phi \right)}_E={\left[D\right]}^{-1}\left(B\right)$, здесь матрица $\left[D\right]$ и вектор $\left(B\right)$ представляются следующим образом:

$\left[D_{ij}\right]=\left(\begin{array}{ccc}\sum _{f=1}^F{\left(\Delta x_f\right)}^2& \sum _{f=1}^F\Delta x_f\Delta y_f& \sum _{f=1}^F\Delta x_f\Delta z_f\\ \sum _{f=1}^F\Delta x_f\Delta y_f& \sum _{f=1}^F{\left(\Delta y_f\right)}^2& \sum _{f=1}^F\Delta y_f\Delta z_f\\ \sum _{f=1}^F\Delta x_f\Delta z_f& \sum _{f=1}^F\Delta y_f\Delta z_f& \sum _{f=1}^F{\left(\Delta z_f\right)}^2\end{array}\right),\left\{B_j\right\}=\left(\begin{array}{c}\sum _{f=1}^F\Delta x_f\Delta {\phi }_f\\ \sum _{f=1}^F\Delta y_f\Delta {\phi }_f\\ \sum _{f=1}^F\Delta z_f\Delta {\phi }_f\end{array}\right)$. (11)

Здесь $\Delta {\phi }_f={\phi }_E-{\phi }_P$ является приращением величины ${\phi }_P^{}$ на грани f через ячейки, которые она разделяет. $\Delta x_f,\Delta y_f,\Delta z_f$ определяют расстояние между центрами грани f и ячейки P.

В случае использования МНК матрица $\left[D\right]$ может быть умножена на вес ${\omega }_f$ грани f:

1. ${\omega }_f^2=1$; (12)
2. ${\omega }_f^2=\frac{1}{\left|\Delta {\overrightarrow{R}}_{Ef}\right|}$; (13)
3. ${\omega }_f^2=\frac{1}{{\left|\Delta {\overrightarrow{R}}_{Ef}\right|}^2}$; (14)
4. ${\omega }_f^2=\frac{\lambda \left|S_f\right|}{V_P\left|\Delta {\overrightarrow{R}}_{Ef}\right|}$. (15)

Исследование точности методов вычисления градиента выполняется на блочно-структурированных сетках, применяемых для моделирования течения в квадратном канале, вблизи пластины и при обтекании профиля крыла (рис. 3).

 

Рис. 3. Геометрия области и расчетные сетки для исследования точности

Fig. 3. Geometry of the area and computational grids for accuracy studies

 

 Градиент численно находится для следующих распределений.

1. Распределение линейного типа: $\phi (x,y,z)=a(x-x_0)+b(y-y_0)+c(z-z_0)+d$. (16)

Точное значение градиента: $\frac{\partial \phi }{\partial x}=a,\frac{\partial \phi }{\partial y}=b,\frac{\partial \phi }{\partial z}=c.$ (17)

2. Распределение квадратичного типа:

$\begin{array}{l}\phi (x,y,z)=a{(x-x_0)}^2+b{(y-y_0)}^2+c{(z-z_0)}^2+d(x-x_0)(y-y_0)+e(x-x_0)(z-z_0)+\\ +f(y-y_0)(z-z_0)+g(x-x_0)+h(y-y_0)+i(z-z_0)+j\end{array}$ (18)

Точное значение градиента:

$\begin{array}{l}\frac{\partial \phi }{\partial x}=2a(x-x_0)+d(y-y_0)+e(z-z_0)+g,\\ \frac{\partial \phi }{\partial y}=2b(y-y_0)+d(x-x_0)+f(z-z_0)+h,\\ \frac{\partial \phi }{\partial z}=2c(z-z_0)+e(x-x_0)+f(y-y_0)+i.\end{array}$ (19)

Для определения точности методов вычисления градиента в каждой ячейке i расчетной сетки используются следующие критерии [13]:

1. Отклонение по длине, %: ${\delta }_i=\frac{\left|\left|\overrightarrow{R_{\tau }}\right|-\left|\overrightarrow{R_{}}\right|\right|}{\left|\overrightarrow{R_{\tau }}\right|}\ast 100\%$; (20)
2. Отклонение по углу, %: ${\phi }_i=\arccos \left(\frac{\left(\overrightarrow{R_{\tau }}\cdot \overrightarrow{R_{}}\right)}{\left|\overrightarrow{R_{\tau }}\right|\cdot \left|\overrightarrow{R_{}}\right|}\right)\ast 100\%$; (21)

В выражениях (20)-(21) – $\overrightarrow{R_{}}$ вычисленное, а $\overrightarrow{R_{\tau }}$ – точное значение градиента. Для рассматриваемой сеточной модели точность расчета градиента определяется путем оценки величин:

${\delta }_{\min }=\underset{i}{\min }{\left(\frac{\left|\left|\overrightarrow{R_{\tau }}\right|-\left|\overrightarrow{R_{}}\right|\right|}{\left|\overrightarrow{R_{\tau }}\right|}\ast 100\%\right)}_i$; (22)

${\delta }_{\max }=\underset{i}{\max }{\left(\frac{\left|\left|\overrightarrow{R_{\tau }}\right|-\left|\overrightarrow{R_{}}\right|\right|}{\left|\overrightarrow{R_{\tau }}\right|}\ast 100\%\right)}_i$; (23)

${\phi }_{\min }=\underset{i}{\min }{\left(\arccos \left(\frac{\left(\overrightarrow{R_{\tau }}\cdot \overrightarrow{R_{}}\right)}{\left|\overrightarrow{R_{\tau }}\right|\cdot \left|\overrightarrow{R_{}}\right|}\right)\ast 100\%\right)}_i$; (24)

${\phi }_{\max }=\underset{i}{\max }{\left(\arccos \left(\frac{\left(\overrightarrow{R_{\tau }}\cdot \overrightarrow{R_{}}\right)}{\left|\overrightarrow{R_{\tau }}\right|\cdot \left|\overrightarrow{R_{}}\right|}\right)\ast 100\%\right)}_i$. (25)

Для численного определения градиента используются:

1) метод Грина-Гаусса на основе интерполяции (6) (Г-Г_1);

2) метод Грина-Гаусса на основе интерполяции (7) (Г-Г_2);

3) МНК, весовая функция (12) (МНК_1);

4) МНК, весовая функция (13) (МНК_2);

5) МНК, весовая функция (14) (МНК_3);

6) МНК, весовая функция (15) (МНК_4).

Градиент вычисляется для следующих выражений.

1. Распределение линейного типа:

$\phi (x,y,z)=5x+7y+9z+10$; (26)

2. Распределение квадратичного типа:

$\phi (x,y,z)=5x^2+7y^2+9z^2+10xy-0.2xz+0.7yz+5.5x+7.3y+8.1z$. (27)

Тестовые расчеты

Тест № 1 – Квадратный канал, структурированная сетка

Рассматривается прямоугольное сечение с равномерной сеткой (рис. 4).

 

Рис. 4. Равномерная сетка

Fig. 4. Uniform grid

 

В табл. 1-2 приведены результаты расчета градиентов по методу Грина-Гаусса и МНК (в том числе, с разными весовыми функциями), получено максимальное значение локальной погрешности не более 5 %.

 

Таблица 1. Погрешность при вычислении градиента для распределения (26), тест № 1

Table 1. Error in calculating the gradient for distribution (26), test № 1

Метод	${\mathit{\delta }}_{\mathbf{m}\mathbf{i}\mathbf{n}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{m}\mathbf{i}\mathbf{n}}$, %	${\mathit{\phi }}_{\mathbf{m}\mathbf{a}\mathbf{x}}$, %
Г-Г_1	0.0	1.5e-7	0.0	0.0
Г-Г_2	0.0	1.5e-7	0.0	0.0
МНК_1	0.0	2.7e-9	0.0	0.0
МНК_2	0.0	2.7e-9	0.0	0.0
МНК_3	0.0	1.3e-9	0.0	0.0
МНК_4	0.0	1.3e-9	0.0	0.0

 

Таблица 2. Погрешность при вычислении градиента для распределения (27), тест № 1

Table 2. Error in calculating the gradient for distribution (27), test № 1

Метод	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
Г-Г_1	0.0	1.99	0.0	1.71
Г-Г_2	0.0	1.99	0.0	1.71
МНК_1	0.0	2.55	0.0	1.85
МНК_2	0.0	2.22	0.0	1.73
МНК_3	0.0	1. 91	0.0	1.51
МНК_4	0.0	1.91	0.0	1.51

 

Тест № 2 – Квадратный канал, структурированная сетка с пристеночными слоями

В случае моделирования течения с учетом турбулентного перемешивания необходимо вблизи поверхности обтекаемого тела строить пристеночные слои, например, как в структурированной сетке на рис. 5.

 

Рис. 5. Структурированная сетка с пристеночными слоями

Fig. 5. Structured grid with near-wall layers

 

В табл. 3-4 приведены результаты расчета градиентов по методу Грина-Гаусса и МНК (с разными весовыми функциями), получено максимальное значение локальной погрешности не более 5 %, что является достаточным при решении задач.

 

Таблица 3. Погрешность при вычислении градиента для распределения (26), тест № 2

Table 3. Error in calculating the gradient for distribution (26), test № 2

Метод	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
Г-Г_1	0.0	0.78	0.0	0.82
Г-Г_2	0.0	0.04	0.0	0.03
МНК_1	0.0	2.51e-07	0.0	0.0
МНК_2	0.0	2.48e-07	0.0	0.0
МНК_3	0.0	4.32e-08	0.0	0.0
МНК_4	0.0	3.83e-08	0.0	0.0

 

Таблица 4. Погрешность при вычислении градиента для распределения (27), тест № 2

Table 4. Error in calculating the gradient for distribution (27), test № 2

Метод	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
Г-Г_1	0.0	2.74	0.0	1.74
Г-Г_2	0.0	2.28	0.0	1.11
МНК_1	0.0	3.02	0.0	3.09
МНК_2	0.0	2.84	0.0	3.01
МНК_3	0.0	2.49	0.0	2.24
МНК_4	0.0	2.31	0.0	1.75

 

Наибольшая величина погрешности получена для распределения (27) в случае расчета градиента методом наименьших квадратов для пристеночной ячейки, характеризующейся наибольшей величиной соотношения сторон (рис. 6).

 

Рис. 6. Ячейка, с наибольшей величиной соотношения сторон

Fig. 6. Сell with the largest aspect ratio

 

Тест № 3 – Пластина, блочно-структурированная сетка

Рассмотрим блочно-структурированную сетку, применяемую при исследовании обтекания пластины (рис. 7).

 

Рис. 7. Пластина, блочно-структурированная сетка

Fig. 7. Plate, block-structured grid

 

В табл. 5-6 приведены результаты расчета градиентов по методу Грина-Гаусса и МНК (с разными весовыми функциями) на данной сетке.

 

Таблица 5. Погрешность при вычислении градиента для распределения (26), тест № 3

Table 5. Error in calculating the gradient for distribution (26), test № 3

Метод	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
ГГ1	3.3e-8	0.61	2.9e-8	0.14
ГГ2	1.6e-8	0.61	2.3e-8	0.14
МНК1	7.4e-5	5.82	5.3e-4	3.99
МНК2	8.6e-4	5.16	5.1e-4	3.76
МНК3	7.9e-4	4.75	5.1e-4	1.61
МНК4	7.9e-4	3.88	5.1e-4	1.21

 

Таблица 6. Погрешность при вычислении градиента для распределения (27), тест № 3

Table 6. Error in calculating the gradient for distribution (27), test № 3

Метод	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
ГГ1	5.8e-8	0.31	5.3e-8	0.33
ГГ2	4.1e-8	0.19	4.9e-8	0.27
МНК1	8.8e-5	8.66	7.5e-4	7.86
МНК2	8.5e-4	7.99	7.3e-4	7.05
МНК3	8.5e-4	5.81	7.3e-4	5.92
МНК4	8.1e-4	4.74	6.9e-4	4.88

 

Получено, что значение градиента вычислено с достаточной интегральной точностью. Наибольшая величина погрешности появляется в случае вычисления по МНК (8.66 %) в ячейке пристеночного слоя (рис. 8).

 

Рис. 8. Ячейка в пристеночном слое

Fig. 8. Cell in the near-wall layer

 

Анализ результатов позволяет сделать вывод о том, что при вычислении градиента для вытянутых ячеек метод Грина-Гаусса имеет большую точность.

Тест № 4 – Профиль крыла, блочно-структурированная сетка

Рассмотрим блочно-структурированную сетку, применяемую при исследовании обтекания профиля NACA0012 (рис. 9).

 

Рис. 9. Профиль NACA0012, блочно-структурированная сетка

Fig. 9. NACA0012 profile, block-structured grid

 

В табл. 7-8 приведены результаты расчета градиентов по методу Грина-Гаусса и МНК (с разными весовыми функциями) на данной сетке.

 

Таблица 7. Погрешность при вычислении градиента для распределения (26), тест № 4

Table 7. Error in calculating the gradient for distribution (26), test № 4

Метод	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
ГГ1	4.51e-7	8.99	6.31e-6	9.77
ГГ2	4.44e-7	8.72	6.04e-6	8.01
МНК1	2.81e-8	6.56	4.39e-7	5.99
МНК2	2.55e-8	5.88	4.22e-7	5.05
МНК3	1.77e-8	4.05	4.01e-7	4.76
МНК4	1.77e-8	3.98	4.01e-7	4.76

 

Таблица 8. Погрешность при вычислении градиента для распределения (27), тест № 4

Table 8. Error in calculating the gradient for distribution (27), test № 4

Метод	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
ГГ1	2.73e-6	8.77	3.12e-5	8.66
ГГ2	2.06e -6	8.28	3.11e-5	8.05
МНК1	3.76e-8	9.14	6.58e-6	9.21
МНК2	3.54e-8	9.05	6.42e-6	8.42
МНК3	3.51e-8	8.52	6.35e-6	7.41
МНК4	3.42e-8	7.99	6.35e-6	5.95

 

Наибольшая величина локальной погрешности в случае МНК наблюдается в ячейках пристеночных слоев, аналогично результатам теста № 3 (для обоих распределений). В случае Грина-Гаусса наибольшая величина погрешности характерна для ячеек с не ортогональными гранями (рис. 10).

 

Рис. 10. Ячейка с неортогональными гранями

Fig. 10. Cell with non-orthogonal edges

 

Получено, что метод Грина-Гаусса и МНК позволяют вычислять значение градиента произвольной функции с достаточной точностью, которая при этом может зависеть от формы контрольного объема.

Гибридный метод вычисления градиента

Для повышения точности вычисления градиента произвольных функций при решении промышленных задач на неструктурированных сетках предлагается применять гибридный метод, где базовыми являются метод Грина-Гаусса и МНК. В этом случае значение градиента определяется путем сложения значений градиента, вычисленных методом Грина-Гаусса и МНК. При этом каждое значение градиента берется с учетом весовой функции (28):

$\nabla {\phi }_P^{}=\beta \nabla {\phi }_P^{LSQ}+(1-\beta )\nabla {\phi }_P^{GG}$, (28)

здесь $\nabla {\phi }_P^{LSQ}$ – значение градиента по МНК, $\nabla {\phi }_P^{GG}$ – значение градиента по методу Грина-Гаусса, $\beta$ – весовая функция.

Ключевой особенностью данного подхода является выбор весовой функции β, которая должна иметь зависимость от геометрической формы контрольного объема. Представим β в виде следующего произведения:

$\beta =\beta {}_{AspectCell}*\beta {}_{curv}$, (29)

здесь $\beta {}_{AspectCell}$ – весовая функция, учитывающая соотношение сторон ячейки, $\beta {}_{curv}$ – весовая функция, учитывающая ортогональность граней ячейки. Выражение для определения $\beta {}_{AspectCell}$ имеет следующий вид:

$\beta {}_{AspectCell}=1-(0.0001*AspectCell),0\le \beta {}_{AspectCell}\le 1$. (30)

$AspectCell=\min \left(\mathrm{10000,}Aspect\right)$, (31)

$Aspect=\frac{F_{\max }}{F_{\min }}$, (32)

где $F_{\max }$ и $F_{\min }$ – скалярное произведение (максимальное и минимальное значение) $\overrightarrow{RP}$ и $\overrightarrow{N}$, $\overrightarrow{RP}$ – вектор от центра ячейки P к центру грани f, $\overrightarrow{N}$ – нормаль к грани f.

При записи выражения (30) принимается предположение, что максимальное значение соотношения сторон контрольного объема может достигать значения 10⁴. В этом случае $\beta {}_{AspectCell}$ принимает нулевое значение, и для вычисления градиента используется метод Грина-Гаусса.

Для вычисления $\beta {}_{curv}$ используется выражение:

$\beta {}_{curv}=\frac{th(\alpha *8-\pi )+1}{2}$, (33)

где α – угол между нормалью грани и вектором, соединяющим центры соседних ячеек, π – число Пи.

Значение функции $\beta {}_{curv}$ асимптотически приближается к единице с ростом угла α, и градиент вычисляется по МНК (рис. 11).

 

Рис. 11. График функции β_curv

Fig. 11. Function β_curv

 

Итоговое значение весовой функции β находится в интервале [0;1], где:

1) при β = 1 градиент вычисляется исключительно методом наименьших квадратов;

2) при β = 0 используется только метод Грина-Гаусса;

3) промежуточные значения определяют взвешенную комбинацию обоих методов.

Точность разработанного авторского метода вычисления градиента исследовалась на рассмотренных выше тестах, полученные результаты приведены в табл. 9-10.

 

Таблица 9. Погрешность при вычислении градиента для распределения (26), гибридная схема

Table 9. Error in calculating the gradient for distribution (26), hybrid scheme

Сетка	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{max}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
Тест №1	0.0	4.28e-8	0.0	0.0
Тест №2	0.0	5.3e-2	0.0	6.1e-2
Тест №3	7.19e-5	1.38	6.41e-5	0.73
Тест №4	3.77e-7	4.65	8.21e-5	3.99

 

Таблица 10. Погрешность при вычислении градиента для распределения (27), гибридная схема

Table 10. Error in calculating the gradient for distribution (27), hybrid scheme

Сетка	${\mathit{\delta }}_{\mathbf{min}}$, %	${\mathit{\delta }}_{\mathbf{m}\mathbf{a}\mathbf{x}}$, %	${\mathit{\phi }}_{\mathbf{min}}$, %	${\mathit{\phi }}_{\mathbf{max}}$, %
Тест №1	0.0	2.03	0.0	1.6
Тест №2	0.0	1.99	0.0	1.43
Тест №3	5.82e-5	2.08	1.93e-5	1.84
Тест №4	1.89e-6	6.55	2.88e-5	5.74

 

Предлагаемый метод расчета градиента позволяет уменьшить максимальные значения погрешности для рассматриваемых распределений. Можно заметить некоторое увеличение минимальных значений погрешности, что объясняется «смешиванием» значений по разным методам. Однако эффект снижения максимальных величин имеет более значимое влияние на точность численного решения.

Все показанные здесь численные алгоритмы и методики расчетов реализованы в рамках российского программного обеспечения ЛОГОС, предназначенного для решения комплексных задач в области вычислительной гидродинамики и аэродинамики [14-17].

Заключение

Рассмотрены результаты численного исследования точности методов вычисления градиента на различных сеточных моделях. На основе классических методов (Грина-Гаусса и МНК) разработан комбинированный авторский алгоритм, использующий аддитивную комбинации результатов расчетов по обоим методам с учетом весовой функции, предложенной авторами работы. Описанный подход приводит к уменьшению максимальных значений погрешности вычисления градиента, что позволяет повысит точность численного решения.

Представленный метод может быть рекомендован при разработке численного алгоритма в рамках CAE-моделирования.

Результаты получены при финансовой поддержке национального проекта «Наука и университеты» в рамках программы Минобрнауки РФ по созданию молодежных лабораторий № FSWE-2024-0001 (научная тема: «Разработка численных методов, моделей и алгоритмов для описания течений жидкостей и газов в естественных природных условиях, и условиях функционирования индустриальных объектов в штатных и критических условиях на суперкомпьютерах экса- и зеттапроизводительности»).

About the authors

R. N. Zhuchkov

Russian Federal Nuclear Center ‒ The All-Russian Research Institute of Experimental Physics

Author for correspondence.
Email: rnzhuchkov@vniief.ru
ORCID iD: 0000-0003-2252-6612
Russian Federation, Sarov

A. S. Kozelkov

Russian Federal Nuclear Center ‒ The All-Russian Research Institute of Experimental Physics; Nizhny Novgorod state technical university n.a. R.E. Alekseev

Email: ASKozelkov@vniief.ru
ORCID iD: 0000-0003-3247-0835
Russian Federation, Sarov; Nizhny Novgorod

N. V. Meleshkin

Russian Federal Nuclear Center ‒ The All-Russian Research Institute of Experimental Physics

Email: nvmeleshkin@yandex.ru
ORCID iD: 0009-0004-3143-7632
Russian Federation, Sarov

A. V. Struchkov

Russian Federal Nuclear Center ‒ The All-Russian Research Institute of Experimental Physics; Nizhny Novgorod state technical university n.a. R.E. Alekseev

Email: andrye134@yandex.ru
ORCID iD: 0000-0002-6979-8968
Russian Federation, Sarov; Nizhny Novgorod

References

Supplementary files