
\documentclass[10pt]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}

\def\udcs{519.6} %Здесь автор определяет УДК своей работы
\def\mscs{??X??} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{1}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}

\def\logo{{\bf\huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}

\def\semrtop
     {
  \vbox{
     \noindent\logo
%     \hspace{80mm}\raisebox{1ex}{ISSN 1813-3304 }

     \vspace{5mm}

     \begin{center}
     {\huge СИБИРСКИЕ \ ЭЛЕКТРОННЫЕ} \\[2mm]
     {\huge МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ} \\[2mm]
     {\large Siberian Electronic Mathematical Reports} \\[1mm]
     {\LARGE\tt{http://semr.math.nsc.ru}}\\[0.5mm]
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
     \vspace{-3mm}
     \noindent
     \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa} \\
     \hline\hline
     \end{tabular}

     \vspace{1mm}
     {
      %\flushleft\it Том 16, стр. 144--144 (2019) 
      \hspace{65mm}{\rm\small УДК \udcs}}
%     \newline
%     	{\rm\small DOI~10.33048/semi.2019.16.xxx}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}{\rm\small MSC\ \ \mscs }
  }%\vbox
}

\usepackage{graphics}
\graphicspath{{Pictures/}}

\begin{document}
\renewcommand{\refname}{References}
\thispagestyle{empty}

\title[]{Модификация IDW метода для численного решения задач аэродинамики на сетках большой размерности}
\author{А.В. САРАЗОВ, А.С. КОЗЕЛКОВ}%

\address{Alexey Vladimirovich Sarazov  % Контактные данные всех авторов и место работы указываются только на английском языке
	\newline\hphantom{ask} Russian Federal Nuclear Center—All-Russian Research Institute of Experimental Physics,
	\newline\hphantom{ask} Sarov, Russia,
	\newline\hphantom{ask} Department of Applied Mathematics, Nizhny Novgorod State Technical University n.a. R.E. Alekseev,
	\newline\hphantom{ask} Nizhny Novgorod, Russia}
\email{alex{\_}sar@bk.ru}%

\address{Andrey Sergeevich Kozelkov  % Контактные данные всех авторов и место работы указываются только на английском языке
	\newline\hphantom{ask} Russian Federal Nuclear Center—All-Russian Research Institute of Experimental Physics,
	\newline\hphantom{ask} Sarov, Russia,
	\newline\hphantom{ask} Department of Applied Mathematics, Nizhny Novgorod State Technical University n.a. R.E. Alekseev,
	\newline\hphantom{ask} Nizhny Novgorod, Russia,
	\newline\hphantom{ask} Moscow Aviation Institute, National Research University,
	\newline\hphantom{ask} Moscow, Russia}
\email{askozelkov@mail.ru}%


\thanks{\sc Sarazov A., Kozelkov A.,
	Modification of the IDW method for numerical simulation of aerodynamics problems on large grids}
\thanks{\copyright \ 2024 Саразов А.В., Козелков А.С}
%\thanks{\rm Результаты получены в рамках программы Минобрнауки РФ по созданию молодёжных лабораторий (тема № FSWE-2021-0008 «Гидродинамические и тепломассообменные процессы в элементах ядерных энергетических установок») и при поддержке гранта Президента РФ по государственной поддержке ведущих научных школ РФ НШ-70.2022.1.5.}
\thanks{\rm Результаты получены при поддержке национального проекта «Наука и университеты» в рамках программы Минобрнауки РФ по созданию молодежных лабораторий № FSWE-2024-0001 (научная тема: «Разработка численных методов, моделей и алгоритмов для описания течений жидкостей и газов в естественных природных условиях, и условиях функционирования индустриальных объектов в штатных и критических условиях на суперкомпьютерах экса- и зеттапроизводительности»), а также при поддержке научной программы Национального центра физики и математики, направление № 2 «Математическое моделирование на супер-ЭВМ экса- и зеттапроизводительности. Этап 2023-2025».\\
	The results were obtained with the assistance of the national project «Science and universities» within the framework of the Ministry of Science and Higher Education of the Russian Federation and the program to create young-adult laboratories. FSWE-2024-0001 (research project: «Developing numerical methods, models and algorithms to describe fluid and gas hydrodynamic characteristics in natural conditions, and industrial object operation in standard and emergency conditions on petascale mainframes»), sponsored by the scientific program of the National Center for Physics and Mathematics, section $\sharp$2 «Mathematical Modeling on Zetta-scale and Exa-scale Supercomputers. Stage 2023-2025, as well as within the framework of implementing the Program to create and develop a world class scientific center «Supersonic» for 2020–2025, funded by the Ministry of Science and Higher Education of the Russian Federation (Agreement . 075-15-2022-309 dated 20 April 2022).}

%\thanks{\it Поступила XX января 2023 г., опубликована XX декабря 2023 г.}%

\semrtop \vspace{1cm}
\maketitle {\small
\begin{quote}
\noindent{\sc Abstract. }Работа посвящена вопросам численного моделирования задач обтекания подвижных тел. В качестве базового подхода для расчёта выбран численный метод, основанный на методе конечных объёмов. Приводятся уравнения Навье-Стокса, описывающие течение вязкого сжимаемого газа, и схемы дискретизации. Движение границ, как правило, влечёт за собой изменение контрольных объёмов, поэтому возникает необходимость в использовании методов изменения расчётной сетки. В работе детально рассматривается IDW метод деформирования расчётных сеток. Вычислительная сложность базового метода эквивалента оценке $O\left( n_i\,n_b\right)$ ($n_i$  и $n_b$ – количество внутренних вершин и вершин на поверхности), что вызывает трудности использования в промышленно ориентированных задачах. Ключевым элементом ускорения процедуры изменения расчётной сетки является применение метода быстрых мультиполей. Иерархические структуры данных совместно с методами теории графов позволили снизить сложность алгоритма до $O\left(n_i\,log(n_b)\right)$ и обеспечить качественное деформирование расчётной сетки. Дополнительно предложен метод распределения вычислительной нагрузки, что вкупе повысило скорость работы алгоритма до 80{\%} на характерных задачах авиационной промышленности. \newline

The work is dedicated to the issues of numerical modeling of the flow problems around moving bodies. A numerical method based on the finite volume method was chosen as the basic approach for calculation. The basic equations of the Navier-Stokes system of equations describing the flow of a viscous compressible gas and discretization schemes are given. The movement of boundaries, as a rule, leads to a change in control volumes, so there is a need to use methods for changing the computational mesh. The article discusses in detail the IDW method for deformation of computational mesh. The computational complexity of the basic method is estimated as $O\left( n_i\,n_b\right)$ (where $n_i$  and $n_b$ are the number of internal vertices and boundary vertices), which causes difficulties for industrially oriented problems. A key element in accelerating the procedure for changing the computational mesh is the use of the fast multipole method. Hierarchical data structures with graph theory methods applied to the chosen approach made it possible to reduce the complexity of the algorithm to $O\left(n_i\,log(n_b)\right)$ and ensure high-quality deformation of the computational mesh. Additionally, a method of distributing the computational load is proposed, which increased the algorithm working speed up to 80{\%} on the characteristic problems of the aviation industry.
\medskip

\noindent{\bf Keywords:} numerical simulation, Navier-Stokes equations, IDW, LOGOS software package
 \end{quote}
}


\section{Введение}

Важной составляющей методов расчёта задач аэродинамики с изменяющимися границами помимо численных схем является тот или иной алгоритм деформирования расчётной сетки. К таким задачам можно отнести задачу определения границы флаттера крыла \cite{1} и задачу определения аэродинамического равновесия крыла \cite{2}. В данном типе задач сетка деформируется не настолько, чтобы ее перестраивать на каждом расчётном шаге и достаточно использования методов деформирования сетки, что позволяет существенно сократить время моделирования. Известные методы, используемые для деформации расчётных сеток, можно разделить на две группы. Первая группа представляет набор методов, использующих физические аналогии и информацию о связи вершин для их перемещения (mesh connectivity approaches). Например, в работе \cite{3} используются уравнения равновесия упругого тела. Отдельным случаем такого подхода является метод пружиной аналогии, в котором рёбра сетки моделируются пружинками, жёсткость которых обратно пропорциональна их длине \cite{4}. Также может применяться аналогия с закрученной пружиной, где жёсткость пружинки описывается функцией угла между гранями \cite{5}. Ряд методов основывается на решении уравнений с оператором Лапласа и бигармоническим оператором \cite{6}. Все методы данной группы используют информацию о структуре сетки для того, чтобы построить систему уравнений порядка $n_i$, равного количеству свободных вершин. Затем эта система решается для нахождения деформированного состояния сетки. На практике при расчётах многих задач такой метод требует большой вычислительной мощности.

Ко второй группе относятся методы, использующие известные функции перемещения подвижных вершин, а также изменяемой границы. Один из таких методов – метод радиальных базисных функций, RBF метод (radial basis function) \cite{7}. Он хорошо «справляется» с большими деформациями в том числе и для неструктурированных сеток. Однако, RBF метод, как и выше описанные, требует решения системы уравнений размерность, которой сопоставима с количеством точек на подвижной поверхности. Хотя в данном случае приходится решать систему меньшего размера, чем в методах, описанных выше, матрица системы оказывается плотно заполненной и плохо обусловленной. Решение этой системы методом Гаусса требует $O\left(n_i^3\right)$ операций, что по-прежнему отнимает большое количество вычислительного времени. В результате, решение системы, как правило, требует использования хорошего предобуславливателя, что заметно усложняет применение RBF метода. После решения системы уравнений, требуется вычислить значение интерполяционной функции во всех точках, что эквивалентно $O\left(n_i\,n_b\right)$ операциям ($n_b$ – количество точек на изменяемой поверхности).

Альтернативным методом, относящимся ко второй группе алгоритмов, и, нашедшим множество применений, называется метод обратных взвешенных расстояний \cite{8}, или IDW метод (inverse distance weighting). Он не требует решения системы уравнений для изменения пространственной сетки и является устойчивым к большим деформациям. В отличие от RBF метода, для перемещения вершин сетки здесь используются выражения в виде функции смещения границы. Такой явный подход уменьшает вычислительные затраты и облегчает параллельную реализацию алгоритма.

IDW метод хорошо зарекомендовал себя ввиду простоты реализации безотносительно к типу расчётной сетки. Однако, практика его использования сводится к модельным задачам, расчётные сетки которых содержат малое количество геометрических примитивов, поскольку вычислительная сложность алгоритма сравнима с RBF методом, а именно составляет $O\left(n_i\,n_b\right)$. По этой причине применение подхода в исходном виде для промышленно ориентированных задач нецелесообразно, так как в инженерной практике применяются, как правило, дискретные модели большой размерности.

В данной работе детально рассмотрен IDW метод \cite{8} и его модификации \cite{9} для сокращения общего времени моделирования с сохранением точности расчёта. В статье представлен модифицированный алгоритм обратных взвешенных расстояний для задачи обновления координат вершин расчётной сетки с сохранением топологии связей. Алгоритм характеризуется оптимальным использованием вычислительных ресурсов, а также позволяет увеличить скорость выполнения процедуры изменения расчётной сетки с сохранением качества получаемых результатов \cite{10}. Модификация алгоритма проводилась путём применения рёберно-узловой топологии, иерархических структур данных и подхода мультиполей \cite{11}, \cite{12}, что позволяет снизить вычислительную сложность алгоритма.
На примерах решения задач определения границы флаттера крыла и определения аэродинамического равновесия крыла показана применимость и преимущество предложенной модификации алгоритма по сравнению с исходным и модифицированным подходом \cite{9} в части быстродействия. Показано, что применение распределённых иерархических структур данных совместно с методами теории графов позволяет дополнительно сократить время расчёта задачи до 80{\%}.

\section{Основные уравнения, способ дискретизации и учёт граничных условий}
Нестационарные трёхмерные турбулентные течения вязкого теплопроводного газа описываются системой уравнений Навье-Стокса (\ref{eq_1}), осредненных по Рейнольдсу \cite{13}, которая в консервативной форме, в декартовых координатах, имеет вид (знаки осреднения опущены):

\begin{equation}
\label{eq_1}
\left\{ \begin{array}{l}
\frac{\partial{\rho}}{\partial{t}}+\nabla{\rho\vec{u}}=0\\\\
\frac{\partial{\left(\rho\vec{u}\right)}}{\partial{t}}+\nabla\left(\rho\vec{u}\vec{u}\right)=-\nabla{p}+\nabla\left( \tau_\mu+\tau_t\right)\\\\
\frac{\partial{\rho E}}{\partial{t}}+\nabla\left[\vec{u}\tau-\left(\vec{q_\mu}+\vec{q_t}\right)\right]
\end{array} \right.,
\end{equation}


В системе уравнений (\ref{eq_1}) используются общепринятые обозначения: $\rho$ – плотность; $\vec{u}$ – вектор скорости течения с компонентами $u$, $v$, $w$; $p$ – давление; $E=C_vT+0.5\left( u^2+v^2+w^2\right)$ – полная энергия газа; $h=C_pT+0.5\left( u^2+v^2+w^2\right)$ – полная энтальпия газа; $\tau_{\mu}$ и $\tau_t$ – молекулярная и турбулентная составляющие тензора касательных напряжений соответственно; $q_{\mu}$ и $q_t$ – молекулярная и турбулентная составляющие плотности теплового потока соответственно; $T$ – температура; $C_v=\left( C_pT-R/m\right)$ – удельная теплоёмкость при постоянном объёме; $C_p$ – удельная теплоёмкость при постоянном давлении; $R$ – универсальная газовая постоянная; $m$ – молярная масса газа.

Величины молекулярной составляющей тензора касательных напряжений ньютоновской среды удовлетворяют реологическому закону Ньютона, связывающий тензор вязких напряжений и тензор скоростей деформаций. Компоненты вектора плотности теплового потока связаны с локальным градиентом температуры законом Фурье \cite{13}.

\begin{equation}
\label{eq_2}
\tau_{\mu}=2\mu\left( T\right),  S=\frac{1}{2}\left(\nabla\vec{u}+\left[\nabla\vec{u}\right]^t\right), q_{\mu}=\chi\left( T\right) \nabla T
\end{equation}

Коэффициент динамической вязкости $\mu\left(T\right) $ и коэффициент теплопроводности $\chi\left(T\right) $ в зависимости от температуры потока определяются по формуле Сазерленда \cite{14}:

\begin{equation}
\label{eq_3}
\mu=\mu_0\left(\frac{T}{T_0}\right)\,\frac{T_0+T_s}{T+T_s},  \chi=\chi_0\left(\frac{T}{T_0}\right)\,\frac{T_0+T_s}{T+T_s}
\end{equation}

где $\mu_0$ и $\chi_0$ соответственно динамическая вязкость и коэффициент теплопроводности при температуре $T_0$, $T_s$  – константа Сазерленда.

Система уравнений (\ref{eq_1}) является незамкнутой из-за неизвестной связи одних из основных переменных этой системы $\tau_t$ и $q_t$ с осредненными параметрами течения. Эта связь может быть установлена с помощью дополнительных соотношений, которые в общем случае называются моделями турбулентности. В инженерной практике используют линейные дифференциальные модели турбулентности \cite{15}, \cite{16}, где применяются эмпирические соотношения для коэффициента турбулентной вязкости $\mu_t$, а также вихреразрешающие модели турбулентности крупных вихрей LES (Large Eddy Simulation), DES (Detached Eddy Simulation) \cite{17, 18}.
Для решения данная система уравнений Навье-Стокса должна быть дополнена граничными условиями различных типов, описание которых можно найти в \cite{13, 14}.
Численное решение системы вида (\ref{eq_1}) основано на методе контрольных объёмов. Контрольные объёмы (ячейки сетки) являются произвольными многогранниками, заполняющими расчётную область без зазоров и наложений. Каждый многогранник ограничен произвольным числом граней. Общий вид ячейки показан на рисунке~\ref{fig:Picture1}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture1}
	\caption{Общий вид ячейки сетки}	
	\label{fig:Picture1}
\end{figure}

Сделаем пояснения к рисунку~\ref{fig:Picture1}. Здесь $P$ – центр ячейки, $F$ – одна из граней ячейки, выделенная цветом, $\vec{n}_f$ – нормаль к грани $F$.

Для численного решения методом конечных объёмов выпишем систему уравнений Навье-Стокса для трёхмерного течения теплопроводного вязкого газа в векторной форме:

\begin{equation}
\label{eq_4}
\frac{d}{dt}\int\limits_{\Delta V}WdV+\oint\limits_{\Delta S} \left(F-G \right) dS=\int\limits_{\Delta V}HdV
\end{equation}

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

%\begin{equation}
$$ 
W=\left(\begin{array}{l}
\rho\\
\rho u\\
\rho v\\
\rho w\\
\rho E
\end{array}
\right),  F=\left( \begin{array}{l}
\rho u_n\\
\rho uu_n+pn_x\\
\rho vu_n+pn_y\\
\rho wu_n+pn_z\\
\rho Hu_n+pu_n
\end{array}
\right),  G=\left( \begin{array}{l}
0\\
\tau_{nx}\\
\tau_{ny}\\
\tau_{nz}\\
\tau u+q
\end{array} \right)
$$
%\end{equation}

Здесь $u_n$ – нормальная составляющая скорости, $q$ – тепловой поток, $\tau_{ij}$ – компоненты тензора вязких напряжений.

В системе уравнений (\ref{eq_4}) учёт движения границ контрольного объёма осуществляется путём введения скорости движения грани:

\begin{equation}
\label{eq_5}
\frac{\partial}{\partial{t}}\int\limits_{\Delta V}WdV+\oint\limits_{\Delta S} \left(F-\vec{\dot{x}}W-G \right) dS=\int\limits_{\Delta V}HdV
\end{equation}

где $\vec{\dot{x}}$ – скорость движения границы контрольного объёма \cite{19}.

По сравнению с уравнением для неподвижных сеток, уравнение вида (\ref{eq_5}) включает скорость изменения контрольного объёма и скорость движения его границы. Отметим, что сохранение порядка аппроксимации численной схемы гарантируется соблюдением условия геометрической консервативности (Geometric Conservation Law, GCL), как показано в \cite{20}.

Движение границ в свою очередь требует построения нестационарных схем расчёта для уравнений типа (\ref{eq_5}). Однако непосредственное использование выражения (5) для создания численной схемы не рекомендуется по причине того, что шаг по времени будет определяться условием устойчивости для каждой ячейки сетки. Это приведёт к существенному снижению глобального шага по времени ввиду наличия в расчётной сетке ячеек, обладающих малым объёмом. По этой причине выражение (5) требует модификации.

Один из широко применяемых подходов моделирования нестационарных течений исходит из предположения, что течение потока в каждый момент времени рассматривается как результат установления некоторого стационарного процесса \cite{21}. Поэтому вместо решения уравнений вида (\ref{eq_5}) рассматриваются модифицированные уравнения:

\begin{equation}
\label{eq_6}
\frac{\partial}{\partial{t}}\int\limits_{\Delta V}WdV+\frac{\partial}{\partial{\tau}}\int\limits_{\Delta V}WdV+\oint\limits_{\Delta S} \left(F-\vec{\dot{x}}W-G \right) dS=\int\limits_{\Delta V}HdV
\end{equation}

где $\tau$ - псевдовремя.

При интегрировании будем использовать теорему о среднем. В качестве среднего значения по объёму принимается значение в центре ячейки, а в качестве среднего значения на грани – значение в центре грани. Тогда уравнение (\ref{eq_6}) запишется следующим образом:

\begin{equation}
\label{eq_7}
\frac{\partial}{\partial{t}}\left(WV\right)_P+\frac{\partial}{\partial{\tau}}\left(WV\right)_P+\sum_{f=1}^{N}\left[F\left( W\right)-\vec{\dot{x}}W-G\left(W, \nabla W\right) \right]_fS_f-\left(H\left(W, \nabla W\right)\right)_PV_P^{n+1}=0
\end{equation}

где $n$ - шаг по времени, $N$ - количество граней в ячейке $P$.

Выражение (\ref{eq_7}) можно представить в операторной форме:

\begin{equation}
\label{eq_8}
\frac{\partial}{\partial{t}}\left(WP\right)_P+\frac{\partial}{\partial{\tau}}\left(WP\right)_P+L\left(W\right)=0
\end{equation}

где

\begin{equation}
\label{eq_9}
L\left(W\right)=\sum_{f=1}^{N}\left[F\left(W\right)-\vec{\dot{x}}W-G\left(W, \nabla W\right) \right]_fS_f-\left(H\left(W, \nabla W\right)\right)_PV_P^{n+1}
\end{equation}

Разностный оператор $L\left(W\right)$ включает аппроксимацию конвективных и диффузионных потоков и источника. Построение разностного оператора $L\left(W\right) $ основано на неявных аппроксимациях конвективных и диффузионных потоков. Для аппроксимации конвективных потоков, как правило, используются односторонние разности. Поэтому аппроксимация дифференциального оператора разностным оператором обладает первым порядком по пространственной переменной. Поскольку уравнения являются нелинейными, то для определения решения на новом временном слое используются итерации по нелинейности. Обозначим индексом $\gamma$  параметры, найденные на предыдущей итерации. Тогда, представим значение оператора (\ref{eq_9}) в виде:

\begin{equation}
\label{eq_10}
L\left(W^{\gamma+1}\right)=L\left(W^{\gamma+1}\right)-L\left(W^{\gamma}\right)+L\left(W^{\gamma}\right)=\Delta{L_1\left(W^{\gamma+1}, W^\gamma\right)}+L_2\left(W^{\gamma}\right)
\end{equation}

где

\begin{equation}
\label{eq_11}
L_2\left(W^{\gamma}\right)=\sum_{f=1}^{N}\left[F\left(W\right)-\vec{\dot{x}}W-G\left(W, \nabla W\right) \right]_fS_f-\left(H\left(W, \nabla W\right)\right)_PV_P
\end{equation}

\begin{equation}
\label{eq_12}
\Delta L_1\left(W^{\gamma+1}, W^{\gamma}\right) =\sum_{f=1}^{N}\left[\frac{\partial{F}}{\partial{W}}-\vec{\dot{x}}-\frac{\partial{G}}{\partial{W}}\right]_f^{\gamma}\left(W^{\gamma+1}-W^{\gamma}\right)-\frac{\partial{H^{\gamma}}}{\partial{W_P}}\left(W^{\gamma+1}-W^{\gamma}\right) 
\end{equation}

В результате получим дельта – форму разностных уравнений

\begin{equation}
\label{eq_13}
\left(\frac{2V_P^{n+1}-V_P^n}{\Delta t}+\frac{V_P^{n+1}}{\Delta \tau}\right)\Delta W_P^{\gamma+1}+\sum_{f=1}^{N}\left[\frac{\partial F}{\partial W} - \vec{\dot{x}}-\frac{\partial G}{\partial W}\right]
\Delta W_f^{\gamma+1}S_f-\frac{\partial{H}}{\partial{W_P}}\Delta{W_P^{\gamma+1}}V_P^{n+1}=R\left(W\right) 
\end{equation}

Правая часть уравнений является невязкой уравнений баланса

\begin{equation}
\label{eq_14}
R\left(W\right)=\left(2V_P^{n+1}-V_P^{\gamma}\right)\frac{W_P^n-W_P^{\gamma}}{\delta t}-\sum_{f=1}^{N}\left[F-\vec{\dot{x}}W-G\right]_f^{\gamma}\Delta s_f+H_P^{\gamma}V_P^{n+1}
\end{equation}

Для вычисления конвективных потоков используются две наиболее известные схемы: схема Роу \cite{22} и схемы AUSMP \cite{23}. Данные схемы относят к схемам типа Годунова. В этих схемах сначала проводится реконструкция (восстановление) решения, в результате чего вычисляются примитивные параметры потока слева и справа от грани ячейки. Для вычисления диффузионных потоков по явной схеме используются параметры и производные, определенные на внутренних гранях интерполяцией по их значениям в ячейках.

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

\section{Деформирование расчётной сетки}

Метод обратных взвешенных расстояний позволяет вычислить перемещение вершин расчётной сетки на основе выражения \cite{8}:

\begin{equation}
\label{eq_15}
s\left(x\right)=\frac{\sum w_j\left(x\right)s_j }{\sum w_j\left(x\right)}
\end{equation}

\begin{equation}
\label{eq_16}
w_j\left(x\right)=\left||x-x_{b_j}|\right|^p,\,p<0,
\end{equation}

где $w_j\left(x\right)$ – весовая функция, $s_j$ – перемещение вершины на границе.

Ключевым параметром выражения (\ref{eq_16}) является константа $p$. Уменьшение данного параметра приводит к затуханию деформирования расчётной сетки вдали от подвижных границ. Данная реализация хорошо подходит для небольших деформаций сдвига и, как правило, для таких задач $p\in \left[ -4; -1\right] $ . Однако, для произвольного изменения формы поверхностной сетки алгоритм изменения координат элементов сетки, основанный на выражении (\ref{eq_15}), не рекомендуется ввиду отсутствия поправок на вращение.

Улучшение исходного алгоритма смещения вершин сетки достигается путём корректировки весовой функции и введения функции перемещения \cite{9}. Идея введения подобной функции лежит в предположении, что перемещение $s_j\left(x\right) $ определяется не только перемещением $j$-ой вершины на поверхности, но и положением произвольной точки $x$, для которой производится расчёт смещения. Рассмотрим полное перемещение произвольной точки $r$ твёрдого тела, как показано на рисунке~\ref{fig:Picture2}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture2}
	\caption{Поворот и перемещение твёрдого тела}	
	\label{fig:Picture2}
\end{figure}

Очевидно, что:

\begin{equation}
\label{eq_17}
r’=r’_b+\left(r’-r’_b\right)=r’_b+M\left(r’-r’_b\right),
\end{equation}

где $M$ – матрица поворота на угол между нормалями $n’=Mn$ (рисунок~\ref{fig:Picture2}). Тогда для смещения точки $r$ получим:

\begin{equation}
\label{eq_18}
s\left(r\right)=r’-r=r’+M\left(r’-r’_b\right)=Mr-r+r’_b=Mr-r+b,
\end{equation}

где $b=r’_b-Mr_b=r_b+s\left(r_b\right)-Mr_b$.

Аналогичным образом функция $s_j\left(x\right) $ можно представить, как смещение рассматриваемой точки в результате переноса и поворота расчётной области относительно $j$-ой вершины:

\begin{equation}
\label{eq_19}
s_j\left(x\right)=M_jx-x+b_j
\end{equation}

При этом вектор $b_j$ будет находиться из условия $s_j=M_jx_{b_j}-x_{b_j}+b_j$. Таким образом, применение данной функции к точкам на поверхности точно описывает их перемещения. Это соответствует корректной постановке задачи интерполяции, то есть $s_j\left(x_{b_j}\right)=s_j$.С учётом выражения для функции смещения формула (\ref{eq_15}) будет иметь вид:

\begin{equation}
\label{eq_20}
s\left(x\right)=\frac{\sum w_j\left(x\right)s_j\left(r\right)}{\sum w_j\left(x\right)}
\end{equation}

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

В свою очередь матрицы поворота определяются на основе изменения направления нормалей поверхностной сетки в точках $x_{b_j}$. Нормаль в вершине на поверхности может быть вычислена, как среднее от нормалей смежных граней. Но данный подход даёт хорошие результаты в случае малого отличия площадей смежных граней. Для нерегулярных дискретных моделей соразмерность соседних граней, как правило, не гарантируется, поэтому в зонах локального сгущения поверхностной сетки правильней нормаль вычислять как средне взвешенную сумму нормалей граней, где вес зависит от площади самой грани. Из найденных нормалей до и после сдвига границы на данном шаге, вычисляется соответствующая матрица поворота $M_j$.

В качестве весовой функции берётся более общее по сравнению с (\ref{eq_15}), (\ref{eq_16}) выражение в виде суммы двух степенных членов, что позволяет сохранить качество призматического слоя расчётной сетки в результате деформирования и сглаживать интерполяционную функцию вдали от границы расчётной области. Согласно работе \cite{9} весовая функция может иметь следующий вид:

\begin{equation}
\label{eq_21}
w_j\left(x\right)=A_j\left[\left(\frac{L}{\left||x-x{b_j}|\right| }\right)^a+\left(\frac{\alpha L}{\left||x-x{b_j}|\right| }\right)^b\right] 
\end{equation}

где $A_j$ – вес вершины пропорциональный площади окружающих его граней, $L$ – характерный размер области, $\alpha L$ – масштаб приграничной области, $a$ и $b$ – некоторые показатели степени.

Численный эксперимент показывает \cite{9}, что значения $a=3$, $b=5$ дают хорошие результаты в широком диапазоне деформаций. Значение параметра $L$ также может варьироваться для разных задач, но в общем случае определяется автоматически как максимальное расстояние между элементами сетки. Параметр $\alpha$ отвечает за отношение величины смещения вершин в приграничной области к смещению вдали от неё. Для лучшей работы метода в задачах с большими смещениями границ рекомендуется увеличивать значения этого параметра.

Анализ выражений (\ref{eq_15}), (\ref{eq_16}) и (\ref{eq_20}), (\ref{eq_21}) показывает, что алгоритм деформирования расчётной сетки, основанный на данной модели, весьма требователен к вычислительным мощностям. Так для двумерных задач число вершин на подвижной границе сравнительно невелико, и базовый алгоритм IDW достаточно быстро справляется с задачей. Тем не менее, для трёхмерных задач число точек на поверхности растёт как $O\left(n^{\frac{2}{3}}\right)$, где $n$ – общее число вершин, поэтому вычислительная сложность алгоритма составляет $O\left(n\cdot n^{\frac{2}{3}}\right) $ , и как следствие – эффективность алгоритма стремительно падает. Используя идеи алгоритмов Барнса-Хата \cite{24} или метода быстрых мультиполей \cite{11, 12}, вычислительная сложность алгоритма IDW можно свести к $O\left(n\log n\right)$.

Суть метода быстрых мультиполей (применительно к методу IDW) состоит в предположении, что для вычисления вектора перемещения необязательно использовать полный набор численных характеристик поверхностных точек ввиду свойств степенной функции вида (\ref{eq_21}). Вместо этого предлагается для вычисления перемещения вершин сетки, лежащих на расстоянии нескольких характерных длин обтекаемого тела, применять некий набор точек с весами близкими по характеристикам полному набору. Это предложение может быть обосновано тем, что, как правило, ячейки в ядре расчётной сетки являются более крупными и погрешности подобной аппроксимации перемещения не окажут влияния на качество расчётной сетки. В свою очередь для призматического слоя расчёт вектора перемещения весьма важен. И наибольшее влияние на перемещение оказывают ближайшие точки, лежащие на подвижной границе, когда влияние прочих поверхностных вершин можно учесть аппроксимацией, как показано на рисунке~\ref{fig:Picture3}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture3}
	\caption{Фрагмент расчётной сетки}
	\label{fig:Picture3}
\end{figure}

Сделаем пояснения. На рисунке~\ref{fig:Picture3} представлен фрагмент поверхностной сетки, на которой точки, не принадлежащие поверхности, выделены красным цветом. Будем считать, что $L_2\gg1>L_1\approx 0$  и левая вершина лежит достаточно близко к поверхности. Ближайшая грань к выделенной точке слева отмечена цветом вместе с её вершинами. Согласно (\ref{eq_21}) смещение левой точки будет определяться решающим вкладом вершин вышеупомянутой грани, в то время как влияние выделенных цветом вершин справа существенно ниже. Справедливо так же полагать, что смещения вершин, составляющих ближайшую грань на поверхности для выделенной точки справа, оказывают влияние соразмерное вкладу всех остальных выделенных вершин. По этой причине ниже приведён алгоритм, позволяющий вычислить перемещение произвольной точки, используя подходы аппроксимации расчётных величин.

Ввиду того, что расположение вершин расчётной сетки в общем случае произвольно, поиск набора точек подвижной поверхностной сетки, который вносит наибольший вклад в выражение (\ref{eq_20}), представляет собой довольно сложную задачу перебора. По этой причине предлагается метод, основанный на разбиении множества поверхностных точек на подмножества \cite{25}.

Наиболее подходящей иерархической структурой данных для представления вложенных множеств является бинарное дерево. Ключевой задачей построения бинарного дерева является разделение поверхностных вершин. Ввиду того, что поверхностная сетка зачастую имеет неравномерное распределение вершин, требуется подход, учитывающий данный факт. Создание нового узла дерева предполагает разделение набора поверхностных вершин, соответствующих узлу-родителю, на два равномощных подмножества точек. Дополнительным требованием, соблюдение которого влечёт повышение скорости работы алгоритма, является формирование каждого узла дерева таким образом, чтобы поверхностные вершины представляли собой компактные вложенные подмножества. Для построения наиболее сбалансированного бинарного дерева предлагается учесть рёберно-узловую топологию поверхностной сетки, которая применяется для разбиения графа сетки \cite{25}.

Построение бинарного дерева с учётом связанности вершин поверхностной сетки, ассоциированных с узлами иерархической структуры, позволяет уменьшить радиус описывающей сферы, как показано на рисунке~\ref{fig:Picture4}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture4}
	\caption{Сравнение разбиения графа сетки а) с применение топологии связи, б) без применения топологии связи}
	\label{fig:Picture4}
\end{figure}

Очевидно, что связанность вершин в общем случае приводит к уменьшению радиуса описывающих сфер $\left(R_i<\tilde{R_i}\right)$.
Разбиение множества поверхностных вершин осуществляется до тех пор, пока количество точек на нижнем уровне не станет меньше определённого значения (как правило, 30 – 60). Узлы бинарного дерева соответствующие данным вершинам принято называть листами.
Ключевой этап уменьшения вычислительной сложности алгоритма деформирования расчётной сетки состоит в построении особого набора точек, качественно представляющего ассоциированные вершины узла бинарного дерева. Поэтому по завершении создания структуры данных для каждого узла выполняются следующие шаги:
\begin{itemize}
	\item[\sf ---] 	создаётся сфера, содержащая всё множество ассоциированных вершин, вычисляется её радиус и центр;
	\item[\sf ---] 	формируется тестовая сфера с центром, совпадающим с центром первой, но с увеличенным радиусом;
	\item[\sf ---] 	на тестовой сфере формируется набор точек и для каждой вычисляется значение, определяемое выражением:
	\begin{equation}
	\label{eq_22}
	f\left(\vec{r}\right)=\sum_{i=1}^{N}\frac{a_i}{\left||\vec{r}-\vec{r_i}|\right|^3},
	\end{equation}
\end{itemize}

где $a_i$ – нормированный вес вершины, $\vec{r_i}$ – координаты ассоциированной вершины подвижной поверхности, $\vec{r}$ – координаты вершины тестовой сферы.

Общепринятым требованием для генерации выборки точек на тестовой сфере является их равномерное распределение. Например, это могут быть вершины додекаэдра или икосаэдра, вписанного в тестовую сферу. На основании вычисленных значений функции вида (\ref{eq_22}) для каждого узла дерева на основе регрессионного анализа строятся 4 $q$-точки, аппроксимирующие функцию $f\left(\vec{r}\right)$:

\begin{equation}
\label{eq_23}
f\left(\vec{r}\right)\approx\frac{1}{4}\left(\frac{1}{\left||\vec{r}-\vec{q_1}|\right|^3}+\frac{1}{\left||\vec{r}-\vec{q_2}|\right|^3}+\frac{1}{\left||\vec{r}-\vec{q_3}|\right|^3}+\frac{1}{\left||\vec{r}-\vec{q_4}|\right|^3}\right)
\end{equation}

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

Таким образом, ускоренный алгоритм вычисления перемещения произвольной вершины состоит из следующих шагов. Вычисляется расстояние от вершины до центра корня бинарного дерева. В том случае, если расстояние в 3 раза больше, чем радиус тестовой сферы, перемещение вычисляется на основе (\ref{eq_22}) с использованием только q-точек корневого узла. В противном случае выполняется спуск по бинарному дереву, и процедура повторяется для каждого из потомков. В том случае, если узел дерева является листом, вклад в перемещение производится всем набором ассоциированных вершин. По завершении спуска по бинарному дереву вычисляется итоговое перемещение на основании вкладов каждого из узлов. Такой подход позволяет достаточно быстро произвести вычисление перемещения любой внутренней вершины благодаря сформированным наборам $q$-точек.

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

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture5}
	\caption{Сопоставление индексов поверхностных вершин узлам дерева}
	\label{fig:Picture5}
\end{figure}

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

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

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

Следует отметить, что в общем случае характер изменения формы границы заранее не известен, поэтому обновление набора $q$-точек и соответствующих ему массивов данных следует делать перед каждой процедурой деформирования расчётной сетки. Данная процедура позволяет существенно повысить качество деформированной расчётной сетки. Однако, задачи, имеющие практический интерес, характеризуются достаточно подробными расчётными сетками. Поэтому соответствующее бинарное дерево так же будет содержать большое количество узлов. Так, например, для поверхности, содержащей 300 тысяч вершин, бинарное дерево будет иметь 13 слоёв (будем считать, что количество вершин в листах 32), и для проведения регрессионного анализа для узлов каждого слоя требуется полный набор поверхностных точек. Тогда вычислительную сложность этапа обновления данных бинарного дерева можно оценить как $O\left(n_b\cdot n_t\cdot n_l\right)$, где $n_t$  – количество точек на тестовой сфере, $n_l$ – количество слоёв структуры данных. Следовательно, процедура обновления полного набора данных для $q$-точек становится затратной с вычислительной точки зрения для подробных дискретных моделей, что является узким местом.

Устранение этого недостатка лежит в распределении нагрузки этапа регрессионного анализа по вычислительным ядрам. Подход разделения вычислительных мощностей основан на свойствах иерархической структуры данных. Предложенное выше бинарное дерево по построению симметрично и обладает детерминированным количеством уровней. Это позволяет реализовать быстрые алгоритмы доступа к памяти, которые применяются для равномерного распределения нагрузки на вычислительные ядра. Одним из таких подходов является алгоритм формирования числового кода для каждого узла таким образом, чтобы исходная иерархическая структура данных представляла собой двоичное дерево поиска (числовые коды узлов на рисунке~\ref{fig:Picture5}). Процедура быстрого поиска применяется для формирования списка ключей узлов и доступа к информации соответствующих вершин (рисунок~\ref{fig:Picture6}).

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture6}
	\caption{Распределение вычислительной нагрузки}
	\label{fig:Picture6}
\end{figure}

На рисунке~\ref{fig:Picture6} приведён пример разбиения узлов бинарного дерева на 3 вычислительных ядра. Разделение осуществляется с учётом объёма вычислений, связанных с ассоциированными вершинами, для каждого узла. Так, вычислительное ядро с индексом 0 обрабатывает только корневой узел, ядро с индексом 1 обрабатывает прямых потомков корневого узла, а ядро с индексом 2 производит вычислительные операции для листов бинарного дерева. В общем случае подобный подход позволяет существенно сократить время, требуемое на регрессионный анализ, что будет показано на примере решения характерных задач авиационной промышленности.

\section{Численные эксперименты}

Представленные выше методы и алгоритмы реализованы в пакете ЛОГОС \cite{26, 27, 28}. Пакет программ ЛОГОС – российский пакет программ инженерного анализа (CAE), предназначенный для решения широкого круга индустриальных задач \cite{26, 27, 28}. Все расчёты, представленные в данной статье, выполнены с использованием многопроцессорной ЭВМ с распределённой памятью в пакете программ ЛОГОС.

\subsection{Задача аэроупругого равновесия крыла HIRENASD}

Проводится расчёт аэроупругих характеристик модели крыла коммерческого самолёта в условиях трансзвукового потока при числе Рейнольдса 23,5e+6 и аэродинамических нагрузках, соответствующих реальным условиям полёта большого транспортного самолёта \cite{2}. Данная задача является многодисциплинарной и требует применения двух расчётных методик, а именно: методики расчёта задач аэродинамики и методики расчёта напряжённо-деформированного состояния. Геометрия модели задачи аэродинамики представлена на рисунке~\ref{fig:Picture7}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture7}
	\caption{Схема расчётной модели}
	\label{fig:Picture7}
\end{figure}

Дискретная модель задачи аэродинамики содержала 10,4 млн. ячеек, $\sim$300 тыс. точек на поверхности крыла. Общее количество вершин равно $\sim$12,5 млн. На поверхности крыла и фюзеляжа количество призматических слоёв равно 18. Среднее значение $y+$ у поверхности крыла составляет 1.

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

Для описания материала в расчётах использовалась модель упругого материала с параметрами: плотность $\rho=7920$ кг/м$^3$, модуль упругости $E=1,813e{+}11$ Па, коэффициент Пуассона $\nu=0,33$. Структурная модель содержит незначительные модификации, заменяющие фюзеляж. Узлы модели, принадлежащие заменителю фюзеляжа, зафиксированы на протяжении всего расчёта.

Условия набегающего потока при моделировании трансзвукового обтекания соответствовали условиям аэродинамического эксперимента \cite{2}:

\begin{itemize}
	\item[\sf ---] 	$M_{\infty}=0,8$ --- число Маха;
	\item[\sf ---] 	$\alpha=2^{\circ}$ --- угол атаки;
	\item[\sf ---] 	$T_{\infty}=204$ K --- температура набегающего потока;
	\item[\sf ---] 	$P_{\infty}=201000$ Па --- статическое давление набегающего потока;
	\item[\sf ---] 	$V_{\infty}=229,06$ м/с --- скорость набегающего потока;
	\item[\sf ---] 	$Re=23.5e{+}6$ --- число Рейнольдса.
\end{itemize}

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

На рисунке~\ref{fig:Picture8} приведено распределение давления на поверхности деформированного крыла и фюзеляжа.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture8}
	\caption{Распределение поля давления}
	\label{fig:Picture8}
\end{figure}

На рисунке~\ref{fig:Picture9} приведено поле напряжений по Мизесу в прочностной части расчета для крыла в состоянии аэроупругого равновесия.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture9}
	\caption{Поле напряжений по Мизесу}
	\label{fig:Picture9}
\end{figure}

Приведённые рисунки демонстрируют корректное распределение расчётных полей, характерные для состояния аэроупругого равновесия крыла при условиях, соответствующих постановке задачи \cite{2}.

В численном эксперименте оценивались безразмерные коэффициенты подъёмной силы ($C_l$), лобового сопротивления ($C_d$) и распределение коэффициента давления ($C_p$) в сечении крыла, а также деформация вдоль размаха крыла ($dS$) в состоянии аэроупругого равновесия.

В таблице~\ref{table:table_1} представлено сравнение с экспериментальными данными интегральных характеристик, а также максимального смещения крыла, полученных в численном эксперименте.

\begin{table}[h!]
	\caption{Сравнение интегральных характеристик крыла HIRENASD}
	\begin{tabular}{|c|c|c|}
		\hline
		& Эксперимент & ЛОГОС \\
		\hline	
		$C_l$ & 0,358 & 0,361 \\
		\hline	
		$C_d$ & 0,0115 & 0,0113 \\
		\hline	
		$dS$, м & 0,0306 & 0,0319 \\
		\hline	
	\end{tabular}
	\label{table:table_1}
\end{table}

На рисунке~\ref{fig:Picture10} представлено сравнение полученного коэффициента давления с экспериментальными данными.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture10}
	\caption{График коэффициента давления по сечению}
	\label{fig:Picture10}
\end{figure}

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

На рисунке~\ref{fig:Picture11} представлено перемещение контрольной точки на поверхности законцовки крыла относительно первоначального положения в зависимости от шага обмена расчётных модулей.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture11}
	\caption{Перемещение крайней точки крыла (аэродинамика)}
	\label{fig:Picture11}
\end{figure}

Аэроупругое равновесие крыла достигнуто за 10 итераций обмена. Отметим, что на каждом шаге обмена аэродинамический модуль решает задачу стационарного обтекания деформированного крыла. Форма поверхности крыла определяется путём решения задачи напряжённо-деформированного состояния крыла при воздействии аэродинамических нагрузок. При данной постановке сходимость аэродинамической части задачи на каждом шаге обмена достигается за 500 итераций. При использовании 96 вычислительных ядер среднее время выполнения итерации составляет 1,1 с. Следовательно, можно считать, что на решение только задачи обтекания без накладных расходов на алгоритм деформирования расчётной сетки и расчёт напряжённо-деформированного состояния потребовалось 5000 итераций или $\sim$1,5 часа.

Дополнительно в численном эксперименте оценивалась производительность нескольких модификаций метода IDW для различного количества вычислительных ядер. Среди рассмотренных модификаций для решения задачи применялись:

\begin{itemize}
	\item[\sf ---] 	модификация № 1 предполагает применение одного бинарного дерева для всех подвижных поверхностей без распределения вычислительной нагрузки для регрессионного анализа;
	\item[\sf ---] 	модификация № 2 предполагает применение одного бинарного дерева для всех подвижных поверхностей с использованием распределения вычислительной нагрузки для регрессионного анализа;
	\item[\sf ---] 	модификация № 3 предполагает построение бинарных деревьев для каждой подвижной поверхности с использованием распределения вычислительной нагрузки для регрессионного анализа.
\end{itemize}

Полученные результаты представлены в таблицах~\ref{table:table_2},~\ref{table:table_3}.

\begin{table}[h!]
	\caption{Среднее время выполнения деформирования расчётной сетки для задачи аэроупругого равновесия крыла HIRENASD}
	\begin{tabular}{|c|c|c|c|}
		\hline
		Число ядер & Модификация №1, с & Модификация №2, с & Модификация №3, с \\
		\hline	
		24 & 53,4 & 50,1 & 43,6 \\
		\hline	
		48 & 28,5 & 24,8 & 22,9 \\
		\hline	
		96 & 17,2 & 14,3 & 11,6 \\
		\hline	
		192 & 10,1 & 7,5 & 6,4 \\
		\hline	
		384 & 7,6 & 4,8 & 4,3 \\
		\hline	
	\end{tabular}
	\label{table:table_2}
\end{table}

\begin{table}[h!]
	\caption{Среднее время выполнения расчёта q-точек для задачи аэроупругого равновесия крыла HIRENASD}
	\begin{tabular}{|c|c|c|}
		\hline
		Число ядер & Параллельный режим, с & Скалярный режим, с \\
		\hline	
		24 & 0,083 & 1,97 \\
		\hline	
		48 & 0,041 & 1,97 \\
		\hline	
		96 & 0,021 & 1,97 \\
		\hline	
		192 & 0,012 & 1,97 \\
		\hline	
		384 & 0,0057 & 1,97 \\
		\hline	
	\end{tabular}
	\label{table:table_3}
\end{table}

В свою очередь время выполнения неоптимизированного метода IDW \cite{8} при использовании 96 вычислительных ядер потребовало более 25 минут. Представленные результаты показывают, что модификация №3 алгоритма деформирования расчётной сетки позволяет получать решение задачи, хорошо согласующееся с экспериментальными данными, а накладные расходы несущественны.

\subsection{Расчёт флаттера крыла AGARD 445.6}

Рассматривается задача определения границы флаттера крыла AGARD 445.6 в потоке вязкого газа \cite{1}. Локализация границы флаттера крыла проводилась для чисел Маха 0,499 – 1,141 при нулевом угле атаки. Задача определения границ флаттера так же является многодисциплинарной задачей, поэтому применяется технология связанного расчёта. Геометрия модели представлена на рисунке~\ref{fig:Picture12}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture12}
	\caption{Схема расчётной модели}
	\label{fig:Picture12}
\end{figure}

В сечении, параллельном направлению потока, модель имеет симметричный профиль NACA65A004. Основные размеры модели крыла показаны на рисунке~\ref{fig:Picture13}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture13}
	\caption{Профиль и проекция крыла AGARD 445.6}
	\label{fig:Picture13}
\end{figure}

Дискретная модель задачи аэродинамики содержала 5,7 млн. неструктурированных ячеек, ~230 тыс. точек на поверхности крыла. Общее количество вершин составило $\sim$6 млн. Количество призматических слоёв равно 18. Среднее значение $y{+}$ у поверхности крыла равно 1.

Параметры набегающего потока рассчитывались на основе чисел Маха и Рейнольдса с варьированием опорного давления. Метод варьирования опорного давления подразумевает изменение давления набегающего потока при постоянной скорости звука для получения различных значений плотности и динамического давления в окрестности границы флаттера крыла. Воздух рассматривается как идеальный газ, для описания турбулентных характеристик использовалась модель Спаларта – Аллмараса.

Для описания поведения материала в расчётах использовалась модель ортотропного материала. Оси ортотропии заданы вдоль размаха крыла и вдоль хорды. Плотность материала $\rho=381,98$ кг/м${^3}$, модуль упругости в продольном направлении $E_1=3,15e{+}9$ Па, модуль упругости в поперечном направлении $E_2=4,2e{+}8$ Па, коэффициент Пуассона в продольном направлении и в поперечном направлении с плоскости крыла, в поперечном направлении перпендикулярно плоскости крыла $\nu=0,31$, модуль сдвига в каждой плоскости $G=4,12e{+}8$ Па.

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

Для количественной характеристики флаттера используются такие величины, как скоростной индекс флаттера (Flutter Speed Index, FSI) и отношение частот (Frequency Ration, FR). На рисунках~\ref{fig:Picture14},~\ref{fig:Picture15} представлены графики коэффициентов в зависимости от числа Маха в сравнении с экспериментальными и расчётными данными \cite{29, 30, 31}.

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture14}
	\caption{Сравнение зависимостей коэффициента FSI от числа Маха}
	\label{fig:Picture14}
\end{figure}

\begin{figure}[h!] 
	\centering
	\includegraphics[width=1.0\textwidth]{Picture15}
	\caption{Сравнение зависимостей коэффициента FR от числа Маха}
	\label{fig:Picture15}
\end{figure}

Как видно из представленных графиков, результаты моделирования хорошо согласуются с экспериментальными данными и данными, полученными в других комплексах программ \cite{29, 30, 31}.

Для более качественного описания процесса изменения формы крыла под действием аэродинамических нагрузок рекомендуемый шаг по времени для каждого расчётного случая подбирался равным $0,005T$, где $T$ – период колебания в эксперименте. Каждый шаг по времени вычисляется на основе 20 итераций по нелинейности, где среднее время одной итерации составило 1,3 с при использовании 96 вычислительных ядер. Также для более точного вычисления количественных характеристик в численном эксперименте желательно получить 5, 6 полных колебаний крыла. Так, например, для режима, соответствующего числу Маха равного 0,678 частота колебания крыла на границе флаттера составляет 17,9 Гц. Следовательно, для определения параметров границы флаттера необходимо моделирование до 0,28 – 0,34 с. Таким образом, ожидаемое время моделирования аэродинамической части связанной задачи без учёта накладных расходов алгоритма деформирования крыла составляет $\sim$8 часов при использовании 96 вычислительных ядер и 20 итераций по нелинейности.
Дополнительно в численном эксперименте оценивалась производительность нескольких модификаций метода IDW для различного количества вычислительных ядер. Модификации рассмотрены ранее.

Полученные результаты представлены в таблицах~\ref{table:table_4},~\ref{table:table_5}.

\begin{table}[h!]
	\caption{Среднее время выполнения деформирования расчётной сетки для задачи локализации флаттера крыла AGARD 445.6}
	\begin{tabular}{|c|c|c|c|}
		\hline
		Число ядер & Модификация №1, с & Модификация №2, с & Модификация №3, с \\
		\hline	
		24 & 20,1 & 18,9 & 16,2 \\
		\hline	
		48 & 10,8 & 9,5 & 9,3 \\
		\hline	
		96 & 7,2 & 5,9 & 5,1 \\
		\hline	
		192 & 4,6 & 3,3 & 2,9 \\
		\hline	
		384 & 3,5 & 2,4 & 2,2 \\
		\hline	
	\end{tabular}
	\label{table:table_4}
\end{table}

\begin{table}[h!]
	\caption{Среднее время выполнения расчёта $q$-точек для задачи локализации флаттера крыла AGARD 445.6}
	\begin{tabular}{|c|c|c|}
		\hline
		Число ядер & Параллельный режим, с & Скалярный режим, с \\
		\hline	
		24 & 0,035 & 1,05 \\
		\hline	
		48 & 0,017 & 1,05 \\
		\hline	
		96 & 0,0092 & 1,05 \\
		\hline	
		192 & 0,0046 & 1,05 \\
		\hline	
		384 & 0,0026 & 1,05 \\
		\hline	
	\end{tabular}
	\label{table:table_5}
\end{table}

Для определения границы флаттера крыла AGARD 445.6 тестирование неоптимизированного метода IDW не проводилось ввиду бесперспективности исходного алгоритма для подобного класса задач. Представленные результаты показывают, что модификация №3 алгоритма деформирования расчётной сетки позволяет получать  решение исходной задачи, согласующееся с экспериментальными и референсными данными \cite{29, 30, 31}.

\section{Анализ полученных результатов}

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

Анализ полученных результатов показал, что построение отдельных иерархических структур данных для каждой границы в отдельности позволяет сократить время работы алгоритма до 40{\%}. Также весомый вклад в быстродействие алгоритма вносит подход распределения вычислительной нагрузки на этапе регрессионного анализа для $q$-точек. В работе представлены результаты ускорения только расчёта $q$-точек (таблицы~\ref{table:table_3},~\ref{table:table_5}). Подобный подход к распараллеливанию также применён и для оценки смещений, матриц поворота, весов аппроксимирующего набора, что вкупе с подходом построения требуемого количества структур данных позволяет достигнуть дополнительного ускорения алгоритма от 20{\%} до 80{\%}.

Следует отметить, что распределение по вычислительным ядрам проводится для ячеек, а не по вершинам дискретной модели. Также обе представленные в работе задачи решались с использованием неструктурированных расчётных сеток со сгущением вблизи подвижных элементов, поэтому разделение вершин по процессорному полю не является равномерным, что приводит к нелинейной зависимости времени выполнения алгоритма деформации от количества процессоров (таблицы~\ref{table:table_2},~\ref{table:table_4}). Тем не менее, время работы метода IDW согласуется с линейным законом.

Важность предложенных подходов к ускорению алгоритма интерполяции IDW методом продемонстрирована на двух характерных задачах. Выбранные задачи разные, поскольку первая задача является квазистационарной, в то время как вторая – нестационарной. По этой причине подходы к моделированию для решения этих задач были разные. Так для задачи определения аэроупругого равновесия крыла выполнен один расчёт, в процессе которого алгоритм деформирования расчётной сетки использовался 10 раз. Накладные расходы на обновление координат вершин для данной задачи составили примерно 2 минуты, что составляет 2{\%} от общего времени моделирования аэродинамической части связанной задачи, что приемлемо. В свою очередь для локализации границы флаттера применялся подход варьирования опорного давления при постоянной скорости звука. Это подразумевает выполнение ряда расчётов для каждого числа Маха с целью получения, как незатухающих, так и затухающих колебаний крыла. Анализ полученных результатов показал, что один расчёт для определения флаттера крыла требует приблизительно 1100 изменений расчётной сетки, что составляет $\sim$1,5 часа или 15{\%} от общего времени моделирования в части аэродинамики. По этой причине выбор подходов к модификации исходного алгоритма очень важен. Так применение отдельных бинарных деревьев для каждой подвижной границы и распределение вычислительной нагрузки на этапе регрессионного анализа позволило снизить общее время расчёта от 3{\%} до 13{\%}, что подчёркивает эффективность модификаций.

Предложенный в настоящей статье адаптированный алгоритм изменения расчётной сетки внедрён в пакет программ ЛОГОС \cite{26, 27, 28} и применяется в качестве штатного подхода для задач с подвижными граничными условиями. Пакет программ ЛОГОС предназначен для решения задач вычислительной гидро-, аэродинамики и теплопереноса и в настоящее время уже широко используется для решения различных классов задач индустриального назначения \cite{26, 27, 28}.

\section{Заключение}

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

Основой ускорения алгоритма является применение иерархических структур данных, таких как бинарные деревья. По причине детерминированности (по построению) структуры данных дальнейшая адаптация алгоритма производилась с применением быстрых алгоритмов доступа к памяти и распределённых вычислений. Также, помимо технических средств ускорения рассматривался подход, предлагающий построение бинарного дерева для каждой подвижной границы в отдельности.

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


\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).


\bibitem{1}  Yates Carson E., Jr.,
{\it AGARD standart aeroelastic configurations for dynamic response I-Wing 445.6},
// Interdisciplinary Research Office NASA Langley Research Center Hampton (USA), 1988, AGARD Report № 765/ 23665-5225, pp. 1-83.

\bibitem{2} Reimer L., Boucke A., Ballmann J., Behr M.,
{\it Computational Analysis of High Reynolds Number Aero-Structural Dynamics (HIRENASD) Experiments},
// IFASD-2009-130, International Forum on Aeroelasticity and Structural Dynamics.- Seattle, WA, June 21-25, 2009.

\bibitem{3} Cavallo P.A., Hosangadi A., Lee R.A., Dash S.M.,
{\it Dynamic unstructured grid methodology with application to aero/propulsive flowfields},
// 15th AIAA Applied aerodynamics Conference, Atlanta, Georgia, 1997, AIAA Paper, $ \sharp $ 1997-2310.

\bibitem{4} Batina J.T.,
{\it Unsteady Euler algorithm with unstructured dynamic mesh for complex-aircraft aerodynamic analysis}
// AIAA J., 1991, v.29, pp. 327-333.

\bibitem{5} Farhat C., Degand C., Koobus B., Lesoinne M.,
{\it Torsional springs for two-dimensional dynamic unstructed fluid meshes},
// Comput. Method appl. M., 1998, v. 163, pp. 231-245.

\bibitem{6} Helenbrook B.,
{\it Mesh deformation using the biharmonic operator},
// Int. J. Numer. Meth. Eng., 2003, v. 56, pp. 1007-1021.

\bibitem{7} Buhmann M.D.,
{\it Radial basis functions: Theory and implementations},
 Cambridge University Press, Cambridge, United Kingdom, 2003.

\bibitem{8} Witteveen J.A., Bijl H.,
{\it Explicit mesh deformation using inverse distance weighting interpolation},
// 19th AIAA Computational Fluid Dynamics Conference, AIAA, San Antonio, Texas, 2009, AIAA Paper, $ \sharp $ 2009-3996.

\bibitem{9} Luke E., Collins E., Blades E.,
{\it A fast mesh deformation method using explicit interpolation},
// Journal of Computational Physics, 2012, v. 231, pp. 586–601.

\bibitem{10} Борисенко О.Н., Черенкова М.В., Гиниятуллина А.Г., Чухманов Н.В, Смолкина Д.Н., Тимаева Т.Е., Блажнова К.А.,
{\it Улучшение качества граней объемных ячеек при генерации неструктурированных сеток в пакете программ «Логос»},
// Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. 2022. Саров. Вып.3. с. 73-85.

\bibitem{11} Greengard L., Rokhlin V.,
{\it A fast algorithm for particle simulations},
// Journal of Computational Physics, 1987, v. 73, pp. 325–348.

\bibitem{12} Sivak S.A., Royak M.E., Stupakov I.M.,
{\it Application of the Fast Multipole Method to Optimization of the Boundary Element Method of Solving the Helmholtz Equation},
// J. Appl. Ind. Math., 2021, v. 15, pp. 490–503, https://doi.org/10.1134/S199047892103011X. 

\bibitem{13} Fletcher C.A.J.,
{\it Computational Techniques for Fluid Dynamics. Second Edn.},
Springer, 1991, Vol I, 401 pp.

\bibitem{14} Лойцянский Л.Г.,
{\it Механика жидкости и газа.},
М.: Наука, 1973. [In Russian].

\bibitem{15} Spalart P.R., Allmaras S.R.,
{\it A one-equation turbulence model for aerodynamic flows}
// AIAA Paper, 1992, № 0439.

\bibitem{16} Menter F.R., Kuntz M., Langtry R.,
{\it Ten years of industrial experience with SST turbulence model},
// Turbulence, Heat and Mass, 2003, v. 4, pp. 625-632.

\bibitem{17} Travin A., Shur M., Strelets M., Spalart P.R.,
{\it Physical and numerical upgrades in the detached-eddy simulation of complex turbulent flows},
// Fluid Mechanics its Application, v.65. Advances in LES of Complex Flows, R. Friedrich, W. Rodi, (editors)/ Proceedings of Euromech Coll. 412. Kluwer-Boston-London, 2002, pp. 239-245.

\bibitem{18} Kozelkov A.S., Krutyakova O. L., Kurulin V.V., Lashkin S.V., Tyatyushkina E. S.,
{\it Application of numerical schemes with singling out the boundary layer for the computation of turbulent flows using eddy-resolving approaches on unstructured grids},
// Computational Mathematics and Mathematical Physics, 2017, v. 57, №6, pp. 1036–1047.

\bibitem{19} Саразов А.В.,
{\it Особенности вычисления конвективных потоков для задач аэродинамики при использовании сеток с подвижными узлами},
// Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов, 2023, Саров, Вып.3, С. 55-66.

\bibitem{20} Thomas P.D., Lombard C.K.,
{\it Geometric conservation law and its application to flow computations on moving grids},
// AIAA Jurnal, 1979, v. 17, № 10, pp. 1030-1037.

\bibitem{21} Blazek J.,
{\it Computational fluid dynamics: principles and applications},
Elsevier, 2001, 470 p., ISBN 008-0-4300-90.

\bibitem{22} Roe P.L.,
{\it Characteristic-based schemes for the Euler equations},
// Annu. Rev. FluidMech., 1986. Vol. 18, pp. 337–365.

\bibitem{23} Weiss J M, Simth W A.,
{\it Preconditioning applied to variable and const density flows},
// AIAA, 1995, v. 33, № 11, pp. 2050?2057.

\bibitem{24} Barnes J., Hut P.,
{\it A hierarchical O(NlogN) force-calculation algorithm},
// Nature, 1986, v. 324, № 4, pp. 446–449.

\bibitem{25} Hendrickson B., Leland R.,
{\it The Chaco user’s guide: Version 2.0},
// Tech. Rep. Sandia National Laboratories, Albuquerque, NM, July 1995, SAND95-2344.

\bibitem{26} Korotkov A., Kozelkov A.,
{\it Three-dimensional numerical simulations of fluid dynamics problems on grids with nonconforming interfaces},
// Siberian Electronic Mathematical Reports, http://semr.math.nsc.ru.

\bibitem{27} Kozelkov A.S., Krutyakova O.L., Kurulin V.V., Strelets D.Yu., Shishlenin M.A.,
{\it The accuracy of numerical simulation of the acoustic wave propagations in a liquid medium based on Navier-Stokes equations},
// Siberian Electronic Mathematical Reports, 2021b, http://semr.math.nsc.ru.

\bibitem{28} Struchkov A.V., Kozelkov A.S., Volkov K., Kurkin A.A., Zhuchkov R.N., Sarazov A.V.,
{\it Numerical simulation of aerodynamic problems based on adaptive mesh refinement method},
// Acta Astronautica, v. 172, July 2020, pp. 7-15.

\bibitem{29} Stanford B., Jacobson K., Chwalowski P.,
{\it Ongoing Aeroelastic Prediction and Validation Activities},
// NASA Langley Research Center. 10.2514/6.2022-1557, 2022.

\bibitem{30} Gupta,
{\it Development of a Finite Element Aeroelastic Analysis Capability},
// Journal of Aircraft, v. 33, № 5, 1996.

\bibitem{31} Pahlavanloo P,
{\it Dynamic Aeroelastic Simulation of the AGARD 445.6 Wing using Edge},
// Tech. Rep. FOI, FOI-R-2259-SE, Stockholm, April 2007.


\end{thebibliography}
\end{document}
