Известия РАН. Механика жидкости и газа, 2020, № 4, стр. 105-116

Численное моделирование обтекания жесткого винта в косом потоке

И. В. Абалакин a, В. Г. Бобков a*, Т. К. Козубская a, В. А. Вершков ab, Б. С. Крицкий b, Р. М. Миргазов b

a Институт прикладной математики им. М.В. Келдыша РАН
Москва, Россия

b Центральный аэрогидродинамический институт им. проф. Н.Е. Жуковского (ЦАГИ)
Жуковский, Россия

* E-mail: veld13@gmail.com

Поступила в редакцию 29.11.2019
После доработки 18.12.2019
Принята к публикации 19.12.2019

Полный текст (PDF)

Аннотация

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

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

Ключевые слова: жесткий несущий винт, вертолет, течение вязкого газа, аэродинамические характеристики винта, численное моделирование, вращающаяся система координат, EBR схема, неструктурированная гибридная сетка

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

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

В данной работе для подтверждения корректности разработанной ранее авторами численной методики моделирования обтекания вращающегося винта вертолета [1] проводится сравнение результатов моделирования жесткого модельного несущего винта на режиме косой обдувки. При этом результаты, полученные с применением методики, реализованной авторами в программном комплексе NOISEtte [2], сравниваются как с результатами натурного эксперимента, так и с данными расчетов, проведенных с использованием коммерческого пакета программ ANSYS CFX (лицензия ЦАГИ № 501024).

Работ, посвященных численному исследованию задачи в подобной постановке, в отечественной литературе не так много. К ним можно отнести статьи [3] и [4], в которых описаны сеточные методы для расчетов аэродинамических характеристик несущего винта вертолета на основе газодинамической системы уравнений Эйлера, а также приведены результаты численного эксперимента. В статьях [58] представлены результаты решения данной задачи с использованием численного моделирования на основе сеточных методов, где в качестве математической модели для описания движения вязкого турбулентного потока сжимаемого газа выбиралась система осредненных по Рейнольдсу уравнений Навье–Стокса, замкнутая различными моделями турбулентности.

В отличие от [3, 4, 8], в данной работе при проведении расчетов с помощью программного комплекса NOISEtte математическая модель, описывающая аэродинамику винта, строится на основе системы осредненных по Рейнольдсу уравнений Навье–Стокса для сжимаемого газа, записанных в неинерциальной, связанной с вращением винта, системе координат. Особенностью реализованной в программном комплексе NOISEtte вычислительной методики является использование “экономного” численного метода повышенной точности с определением переменных в узлах неструктурированных гибридных сеток. Применение гибридных неструктурированных сеток существенно упрощает процедуру построения сеточных моделей, гибко адаптирующихся к изменению геометрии исследуемой конфигурации и, в частности, лопасти винта. “Экономный” численный метод для неструктурированных сеток строится на основе EBR схемы [912], повышенная точность которой обеспечивается за счет квазиодномерной реберно-ориентированной аппроксимации невязких потоков.

1. ФИЗИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ О ВЯЗКОМ ОБТЕКАНИИ ЧЕТЫРЕХЛОПАСТНОГО НЕСУЩЕГО ВИНТА НА РЕЖИМЕ КОСОГО ОБТЕКАНИЯ

Моделируемая конфигурация (рис. 1) соответствует конфигурации, испытанной в серии экспериментов [13], проведенных в ЦАГИ им. проф. Н.Е. Жуковского в 1970-х годах, и представляет собой четырехлопастный жесткий несущий винт вертолета с комплектом плоских лопастей, прямоугольной формы в плане, жестких на изгиб и кручение.

Рис. 1.

Общая конфигурация модельного несущего винта: геометрия, размеры и система координат.

Лопасти построены на базе аэродинамического профиля NACA-0012, хорда лопасти составляет b = 0.15 м, крутка отсутствует, радиус комля лопасти составляет r0 = 0.2 м. Общий шаг лопастей φ составляет 8°. Радиус винта R = 1.2 м, центральное тело представляет собой эллипсоид вращения, имеющий горизонтальный и вертикальный радиусы r1 = 0.04 и r2 = 0.2 м соответственно.

Как в физическом, так и в вычислительном эксперименте, скорость вращения винта равнялась 360 об./мин, что соответствует концевой скорости лопасти Vtip = 45.24 м/с. В работе исследовалось три режима обтекания с различными скоростями набегающего потока Vflow = 6.79, 11.31 и 20.36 м/с с нулевым углом атаки винта. Физические параметры невозмущенного потока выбирались следующими: плотность воздуха ρ0 = 1.225 кг/м3, давление 103 025 Па, температура воздуха 23°C. Число Рейнольдса в расчете определяется как Re = ρ0Vtipb0. При значении динамической вязкости воздуха μ0 = 1.827 × 10–5 Пас ⋅ с, соответствующей температуре 23°C, число Рейнольдса равняется Re = 4.55 × 105.

Расчетная область представляет собой цилиндр радиусом, равным 10 радиусам винта. Ось цилиндра совпадает с осью вращения винта. Нижняя граница расчетной области располагается на расстоянии 15 радиусов винта от плоскости вращения, а верхняя – на расстоянии 10 радиусов.

2. ИССЛЕДУЕМЫЕ АЭРОДИНАМИЧЕСКИЕ ХАРАКТЕРИСТИКИ

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

Определим аэродинамические силы лопасти, следуя монографии [14]. Рассмотрим две системы координат: неподвижную систему каоординат (x, y, z) и систему координат (x', y', z'), связанную с вращающимся винтом. Так как расчет проводится во вращающейся вокруг оси Oz системе координат (см. рис. 1), можно отождествить оси неподвижной и вращающейся систем координат. Таким образом, неинерциальная, вращающаяся система координат характеризуется поворотом неподвижной системы координат на угол ψ(t) + π/2 вокруг оси Oz, где азимутальный угол ψ(t) определяется как ψ(t) = –|ω|t = –ωt.

Так как расчет проводится во вращающейся вокруг оси Oz системе координат (см. рис. 1), можно отождествить оси неподвижной и вращающейся систем координат.

Пусть p(x, y, z, t) есть распределение давления по поверхности винта S, а вектор n – внешняя к поверхности S единичная нормаль с компонентами nx, ny, nz. Тогда сила тяги винта T вычисляется как интеграл давления по поверхности винта

$T(t) = \int\limits_S^{} {p{{n}_{z}}(t)ds} $

Момент силы винта определяется проекциями сил на оси неподвижной системы координат и радиус-вектором (рычагом) координат лопасти во вращающейся системе координат r = (0, 0, z) = = (0, 0, z')

${\mathbf{M}}(t) = \int\limits_S^{} {p(t)({\mathbf{r}} \times {\mathbf{n}})ds} $
откуда следует выражение для аэродинамического крутящего момент винта

${{M}_{k}}(t) = {{M}_{y}}(t) = \int\limits_S^{} {zp{{n}_{x}}(t)ds} $

Коэффициент силы тяги cT и коэффициент аэродинамического крутящего момента mk вычисляются обезразмериванием силы тяги и крутящего момента

${{c}_{T}}(t) = \frac{{2T(t)}}{{{{\rho }_{0}}A{{{(\omega R)}}^{2}}}},\quad {{m}_{k}}(t) = \frac{{2{{M}_{k}}(t)}}{{{{\rho }_{0}}A{{{(\omega R)}}^{2}}R}}$
где ρ0 – плотность невозмущенного воздуха, A = πR2 – площадь диска винта, R – радиус винта, а ω – модуль угловой скорости лопасти.

Также введем коэффициент нормальной силы cn

${{c}_{n}} = \frac{1}{{{{\rho }_{0}}A{{{(\omega R)}}^{2}}}}\int\limits_L^{} {pn_{z}^{'}dl} $
определяемый в нормальном сечении лопасти L, где $n_{z}^{'}$ = $n_{x}^{'}$cosφ – $n_{z}^{'}$sinφ – компонента нормали в системе координат, связанной с лопастью.

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

3. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

3.1. Уравнения Навье–Стокса в неинерциальной системе координат

Для расчета внешнего обтекания потоком газа винта, вращающегося с угловой скоростью ω, используется осредненная по Рейнольдсу система уравнений Навье–Стокса для сжимаемого газа с замыканием по модели турбулентности Спаларта–Аллмараса. Система уравнений рассматривается в неинерциальной вращающейся системе координат. Вращение осей координат происходит вокруг неподвижной оси винта с постоянной угловой скоростью, равной скорости вращения лопасти. Это означает, что обтекаемые потоком газа лопасти винта неподвижны, а направление внешнего потока меняется в зависимости от азимутального угла ψ.

Вид системы уравнений Навье–Стокса без модели турбулентности, записанной для сжимаемого газа в неинерциальной вращающейся системе координат относительно вектора абсолютной скорости, как показано в работе [1], имеет вид

(3.1)
$\begin{gathered} \frac{{\partial \rho }}{{\partial t}} + {\text{div}}\rho ({\mathbf{u}} - {\mathbf{V}}) = 0 \hfill \\ \frac{{\partial \rho {{{\mathbf{u}}}^{{^{{^{{}}}}}}}}}{{\partial t}} + {\text{div}}\rho ({\mathbf{u}} - {\mathbf{V}}) \otimes {\mathbf{u}} + \nabla p = {\text{Div}}{\mathbf{S}} - \rho ({\mathbf{\omega }} \times {\mathbf{u}}) \hfill \\ \frac{{\partial {{E}^{{^{{^{{}}}}}}}}}{{\partial t}} + {\text{div}}({\mathbf{u}} - {\mathbf{V}})E + {\text{div}}{\mathbf{u}}p = {\text{div}}{\mathbf{q}} + {\text{div}}{\mathbf{Su}} \hfill \\ \end{gathered} $
где u – вектор абсолютной скорости в неподвижной системе координат, E – полная энергия, q – поток внутренней энергии, а V = ω × r – вектор линейной скорости вращения, определяемый вектором угловой скорости ω. С точки зрения наблюдателя, находящегося в неподвижной системе координат, система уравнений (3.1) описывает изменение консервативных переменных за счет их переноса во вращающейся со скоростью V среде, градиента давления и повтора вектора скорости на азимутальный угол ψ(t) = –|ω|t.

В настоящей работе методика моделирования турбулентных течений на основе RANS требует введения понятия коэффициента турбулентной вязкости μT, для определения которой система RANS дополняется эволюционным дифференциальным уравнением Спаларта–Аллмараса для величины $\tilde {\nu }$

$\frac{{\partial \rho \tilde {\nu }}}{{\partial t}} + \frac{{\partial \rho \tilde {\nu }{{u}_{i}}}}{{\partial {{x}_{i}}}} = {{D}_{\nu }} + {{G}_{\nu }} - {{Y}_{\nu }}$

Вид членов Dν, Gν, Yν, описывающих, соответственно, диффузию, генерацию и диссипацию турбулентности, можно найти, например, в статье [15]. Коэффициент турбулентной вязкости μT определяется согласно модели Спаларта–Аллмараса следующим образом

${{\mu }_{T}} = \rho \tilde {\nu }\frac{{{{\chi }^{3}}}}{{{{\chi }^{3}} + 357.911}},\quad \chi = \frac{{\rho \tilde {\nu }}}{\mu }$

Описанная модель используется при проведении расчетов в программном комплексе NOISEtte. Моделирование турбулентного обтекания несущего винта на режиме горизонтального полета в коммерческом пакете программ ANSYS CFX проводится также на основе уравнений RANS, однако в качестве замыкания используется модель турбулентности SST [16].

3.2. Граничные условия на основе пристеночных функций

Для корректной постановки граничных условий на твердой поверхности при наличии турбулентного пограничного слоя в каждом сеточном узле, лежащем на поверхности, должна быть известна величина пристеночного числа Рейнольдса

${{y}^{ + }} = \frac{{\rho {{u}_{*}}{{h}_{{nw}}}}}{\mu }$
где hnw – характерный размер по нормали к твердой поверхности пристеночных ячеек, а ${{u}_{*}}$ – динамическая скорость, определяемая через напряжение трения на стенке τw как

(3.2)
${{u}_{*}} = \sqrt {\frac{1}{{{{\rho }_{w}}}}{{\tau }_{w}}} $

Для всех используемых в расчетах сеток величина y+ была больше единицы. Это означает, что постановка граничных условий прилипания (скорость жидкости на поверхности совпадает со скоростью движения самой поверхности) не корректна, так как невозможно правильное определение градиентов компонент вектора скорости на твердой стенке. Поэтому в данной ситуации используется механизм пристеночных функций [17, 18], определяющий поле скорости турбулентного потока, которое моделирует течение в неразрешенной сеткой части пограничного слоя. Реализация механизма пристеночных функций зависит от используемых методов дискретизации уравнений Навье–Стокса. Далее приведена конкретная реализация, используемая в представленных расчетах.

Определим вязкую составляющую потока импульса на граничной поверхности через напряжение трения на стенке или с учетом соотношения (3.2) через динамическую скорость следующим образом

(3.3)
${\mathbf{F}}_{w}^{{wl}} = \rho u_{f}^{2}{\mathbf{u}}_{t}^{e}$
где ${\mathbf{u}}_{t}^{e}$ – единичный вектор тангенциальной скорости на граничной поверхности, равный
${\mathbf{u}}_{t}^{e} = {{{\mathbf{u}}}_{{\mathbf{t}}}}/{\text{|}}{{{\mathbf{u}}}_{{\mathbf{t}}}}{\text{|}},\quad {{{\mathbf{u}}}_{{\mathbf{t}}}} = {\mathbf{u}} - {{u}_{n}}{\mathbf{n}}$
а un – нормальная скорость на граничной поверхности, которая полагается равной нормальной компоненте скорости движения стенки Vn = (ω × r) ⋅ n.

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

(3.4)
${\text{|}}{{{\mathbf{u}}}_{t}}{\text{|}} - {{u}_{*}}f\left( {\frac{{{{u}_{*}}\rho {{h}_{{nw}}}}}{\mu }} \right) = 0$
где пристеночная функция f(y+), определяющая профиль скорости в вязком, буферном и логарифмическом подслое, задается законом Райхардта [19]

$f({{y}^{{\text{ + }}}}) = \frac{1}{{0.41}}\ln (1 + 0.41{{y}^{ + }}) + 7.8\left( {1 - {{e}^{{ - \frac{{{{y}^{ + }}}}{{11}}}}} - \frac{{{{y}^{ + }}}}{{11}}{{e}^{{ - \frac{{{{y}^{ + }}}}{3}}}}} \right)$

Нелинейное уравнение (3.4) решается методом Ньютона. На каждой итерации метода Ньютона значения ut, μ, ρ берутся в узлах граничной поверхности. Нулевое приближение динамической скорости выбирается равным $u_{f}^{0}$ = $\tau _{w}^{0}$/|ut|, где напряжение трения на твердой стенке $\tau _{w}^{0}$ вычисляется через нормальную производную от модуля тангенциальной скорости с использованием первой разности τw = μ|ut|/δ. Величина δ = hw/2 определяет расстояние до твердой стенки, на которой выполняются условия прилипания, а тангенциальная скорость задается на граничной поверхности.

Окончательное выражение для вектора потока импульса на граничной поверхности представляет собой сумму потока (3.3) с учетом приграничных функций и конвективной составляющей потока, учитывающей условие un = Vn

${{{\mathbf{F}}}_{w}} = p{\mathbf{n}} + {\mathbf{F}}_{w}^{{wl}}$

Описанная техника пристеночных функций используется при проведении расчетов в программном комплексе NOISEtte. В случае пакета ANSYS CFX в качестве граничного условия на поверхностях винта и центрального тела ставится условие прилипания, корректность которого обеспечивается в расчетах за счет использования достаточно подробных вблизи границы сеток.

4. ЧИСЛЕННАЯ МЕТОДИКА И РАСЧЕТНЫЕ СЕТКИ

4.1. Численная методика в программном комплексе NOISEtte

Пространственная дискретизация в программном комплексе NOISEtte базируется на вершинно-центрированной формулировке, означающей, что все искомые переменные определены в узлах гибридной сетки11, вокруг которых построены расчетные ячейки (дуальная сетка). Для аппроксимации конвективных потоков системы уравнений Навье–Стокса используется метод конечных объемов. Под конечным объемом понимается ячейка дуальной сетки. Поток через грань контрольного объема рассчитывается по схеме Роу. При использовании EBR схемы повышение точности пространственной аппроксимации достигается применением квазиодномерной реконструкции переменных на расширенном шаблоне, ориентированном вдоль ребра сетки. Подробно класс EBR-схем описан в [912].

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

Для пространственной аппроксимации вязких членов в системе уравнений Навье–Стокса используется конечно-элементный метод Галёркина на основе линейных базисных функций.

Интегрирование по времени проводится по неявной трехслойной схеме 2 -го порядка аппроксимации с последующей линеаризацией по Ньютону разностной по пространству системы уравнений. На каждой ньютоновской итерации применяется стабилизированный метод бисопряженных градиентов (BiCGSTAB) [20] для решения системы линейных уравнений.

Как было упомянуто в п. 1, расчетная область задачи представляет собой цилиндр с осью, совпадающей с осью вращения винта. В расчетной области строится неструктурированная гибридная сетка, измельченная в областях, требующих подробного моделирования аэродинамических эффектов. Так, наиболее подробная сетка строится на поверхности лопастей винта вблизи передней и задней кромок, а также на торцевой части лопасти и в области примыкания лопастей к центральному телу (рис. 2б, 2в). От поверхности винта внутри расчетной области строятся слои призматической сетки с увеличивающейся в геометрической прогрессии высотой. С удалением от твердых поверхностей по достижении призм изотропных размеров объемная сетка становится тетраэдральной, причем высота тетраэдров увеличивается с удалением от винта (рис. 2а).

Рис. 2.

Общий вид объемной сетки в расчете NOISEtte: продольное сечение (а), сетка вблизи центрального тела (б) и законцовки лопасти (в).

Таким образом, построенная средствами ICEM CFD из программного пакета ANSYS [21] сетка содержит 9.3 млн узлов и 24.3 млн объемных элементов.

4.2. Вычислительная постановка в пакете ANSYS CFX

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

Для проведения расчетов с помощью программы ICEM CFD строится неструктурированная сетка с размером расчетной области 30 × 30 × 23 м. Так как расчеты пакетом ANSYS CFX проводятся в неподвижной системе координат, для моделирования вращения винта в центре расчетной области строится цилиндрическая подобласть радиусом 2 м и высотой 2 м, содержащая геометрическую модель винта. Передача данных между внутренней и внешней подобластями реализуется на основе скользящих интерфейсов. Общее число ячеек расчетной сетки составляет 66 млн, причем из них во внутренней области содержится 51 млн ячеек. Для разрешения пристеночной области вблизи поверхности строится 21-слойная призматическая сетка, обеспечивающая величину числа y+ не более 1. На передней (по потоку) границе расчетной области устанавливается входное граничное условие, на задней – выходное. Боковые, верхняя и нижняя границы расчетной области моделируются свободной границей с заданным фиксированным значением статического давления.

5. РЕЗУЛЬТАТЫ РАСЧЕТОВ

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

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

Общая картина течения во многом определяется концевыми вихрями, сходящими с лопастей, сносящиеся вниз по потоку и взаимодействующие как с набегающими лопастями, так и вихревыми структурами, присутствующими в потоке. При этом положения ядер концевых вихрей повторяют траектории движения концов лопастей с учетом вращения винта и набегающего потока. Это можно видеть на рис. 3, где приведена визуализация вихрей многосвязной изоповерхностью Q-критерия [22] и траектории концов лопастей в абсолютной системе координат в плоскости вращения винта Oxy в фиксированные моменты времени, соответствующие азимутальным углам ψ = 0 (рис. 3а) и ψ = 45 (рис. 3б). Траектории концов лопастей определяются следующим законом

(5.1)
$\begin{gathered} \left( \begin{gathered} x \\ y \\ \end{gathered} \right) = \left( \begin{gathered} R\cos (\omega t + {{\psi }_{k}}) - t{{V}_{{flow}}} \hfill \\ R\sin (\omega t + {{\psi }_{k}}) \hfill \\ \end{gathered} \right), \\ k = 0,1,2,3 \\ \end{gathered} $
где ψk = πk/2 – начальные азимутальные углы соответствующих лопастей винта. Также на рис. 3 изображены контуры избыточного давления pg на поверхности винта.

Рис. 3.

Траектории концевых вихрей – траектории концов лопастей и визуализация изоповерхности Q-критерия при разных азимутальных положениях лопастей: а, б – ψ = ωt = 0, 45°; 1–4 – номер лопасти.

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

Рис. 4.

Распределения коэффициента давления вдоль хорды лопасти, полученные с использованием NOISEtte (1) и ANSYS CFX (2) при различных значениях относительного радиуса r/R и азимутального угла ψ: a–г – (r/R, ψ) = = (0.5, 0), (0.7, 90°), (0.9, 180°), (0.99, 270°).

Значения коэффициента нормальной силы по размаху лопасти для различных азимутальных положений представлены на рис. 5. Результаты, полученные в вычислительных экспериментах, достаточно близки друг к другу и к экспериментальным значениям, однако на обратном ходе лопасти для азимутальных положений ψ = 270°– 360°(0) наблюдается расхождение численных данных, полученных с помощью различных численных методик в NOISEtte и ANSYS CFX, с данными физического эксперимента. Возможными причинами этого являются несовершенство методики измерения давления на поверхности лопасти в присутствии поля центробежных сил, а также неточность измерения азимутального положения лопасти в натурном эксперименте [13]. Справедливость второго утверждения подтверждается сравнением разницы давлений, измеренной в двух характерных точках лопасти на носке сечения с относительной координатой по хорде 0.03, расположенных на верхней и нижней поверхностях эффективного сечения r/R = 0.8. Эволюция данной характеристики за один оборот винта, полученная в эксперименте и в расчете, приведена на рис. 6. Видно, что для обеих скоростей обтекания винта характер поведения и величина экстремумов в расчете получены достаточно точно, но наблюдается систематический сдвиг в виде запаздывания лопасти приблизительно на 10°.

Рис. 5.

Сравнение распределения коэффициентов нормальной силы вдоль размаха лопасти, полученных в эксперименте (1) (см. [13]), с использованием NOISEtte (2) и ANSYS CFX (3) при различных значениях азимутального угла: а–г – ψ = 0, 90°, 180°, 270°.

Рис. 6.

Разницы давления в контрольных точках на передней кромке лопасти, полученные в натурном эксперименте (1) (см. [13]) и в расчете NOISEtte (2) при скорости набегающего потока 6.79 м/с (а) и 11.31 м/с (б).

На рис. 7 приведено сравнение коэффициентов силы тяги (рис. 7а) и коэффициентов компонент момента винта (рис. 7б–7г), полученные в расчетах с использованием кодов NOISEtte и ANSYS CFX. Достаточно хорошее совпадение этих аэродинамических характеристик подтверждается малыми значениями относительно разницы величин максимумов и минимумов, а также разницы в их азимутальном положении, приведенными в табл. 1.

Рис. 7.

Сравнение зависимостей коэффициентов силы тяги (а) и продольной (б), поперечной (в) и нормальной (г) компонент момента винта от азимутального положения лопасти, полученных с использованием NOISEtte (1) и ANSYS CFX (2).

Таблица 1.

Относительная разница экстремумов и их азимутального положения для аэродинамических характеристик, полученных с использованием NOISEtte и ANSYS CFX

Величина Тип экстремума Относительная разница, % Разница в азимутальном положении Δψ (°)
ct min 7.27 2.64
max 0.14 1.68
mx min 3.26 2.68
max 2.85 0.68
my min 2.08 0.17
max 8.92 0.81
mz min 9.40 2.45
max 6.87 0.85

Расчеты, результаты которых приведены в статье, были получены с использованием вычислительных ресурсов ОВК НИЦ “Курчатовский институт” [23]. Вычисления производились с использованием комплекса программ NOISEtte в параллельном режиме с использованием от 200 до 480 вычислительных ядер, при этом использовалось двухуровневое распараллеливание MPI+OpenMP [2].

ЗАКЛЮЧЕНИЕ

Проведено численное моделирование течения вязкого газа около вращающегося жесткого несущего винта вертолета на режиме косого обтекания с использованием “экономной” вычислительной методики, разработанной и реализованной в исследовательском программном комплексе NOISEtte. Проведенное сравнение полученных в численном эксперименте аэродинамических характеристик показало удовлетворительное согласование как с экспериментальными данными, так и с результатами аналогичных расчетов с использованием коммерческого пакета ANSYS CFX. Использование оригинальных численных схем повышенной точности позволило получить подобные результаты, используя расчетную сетку существенно меньшего размера, чем в расчетах ANSYS CFX: 9.2M узлов в расчетах NOISEtte и 66M ячеек в расчетах ANSYS CFX. Валидация, проведенная путем сравнения с экспериментальными данными, подтвердила корректность численной методики программного комплекса NOISEtte применительно к моделированию обтекания вращающегося винта турбулентным потоком и определения его аэродинамических характеристик.

Список литературы

  1. Абалакин И.В., Аникин В.А., Бахвалов П.А., Бобков В.Г., Козубская Т.К. Численное исследование аэродинамических и акустических свойств винта в кольце // Изв. РАН. МЖГ. 2016. № 3. С. 130–145. https://doi.org/10.7868/S0568528116030026

  2. Gorobets A. Parallel Algorithm of the NOISEtte Code for CFD and CAA Simulations // Lobachevskii Journal of Mathematics. 2018, Vol. 39, No. 4, pp. 524–532. https://doi.org/10.1134/S1995080218040078

  3. Карабасов С.А. Использование гибридного метода для моделирования шума от высокоскоростных лопастей вертолета // Матем. моделирование. 2006. Т. 18. № 2. С. 2–23.

  4. Копьев В.Ф., Титарев В.А., Беляев И.В. Разработка методологии расчета шума винтов с использованием суперкомпьютеров // Уч. зап. ЦАГИ. 2014. Т. XLV. № 2. С. 78–106.

  5. Игнаткин Ю.М., Константинов С.Г. Исследование аэродинамических характеристик несущего винта вертолета методом CFD // Труды МАИ. Вып. 57. 2012.

  6. Batrakov A., Garipova L., Kusyumov A., Mikhailov S., Barakos G. Computational fluid dynamics modeling of helicopter fuselage drag // J. Aircraft. 2015. V. 52. № 5. P. 1634–1643. https://doi.org/10.2514/1.C033019

  7. Garipova L., Batrakov A., Kusyumov A., Mikhaylov S., Barakos G. Aerodynamic and acoustic analysis of helicopter main rotor blade tips in hover // Interactional Journal of Numerical Methods for Heat and Fluid Flow. 2016. V. 26. № 7. P. 2101–2118. https://doi.org/10.1108/HFF-08-2015-0348

  8. Копьев В.Ф., Зайцев М.Ю., Воронцов В.И., Карабасов С.А., Аникин В.А. Расчет шума несущего винта вертолета и его экспериментальная проверка на режиме висения // Акустический журнал. 2017. Т. 63. № 6. С. 651–664. https://doi.org/10.7868/S0320791917060077

  9. Абалакин И.В., Козубская Т.К. Схема на основе реберно-ориентированной квазиодномерной реконструкции переменных для решения задач аэродинамики и аэроакустики на неструктурированных сетках // Матем. моделирование. 2013. Т. 25. № 8. С. 109–136.

  10. Бахвалов П.А. Схема с квазиодномерной реконструкцией переменных на сетках из выпуклых многоугольников для решения задач аэроакустики // Матем. моделирование. 2013. Т. 25. № 9. С. 95–108.

  11. Abalakin I., Bakhvalov P., Kozubskaya T. Edge-based reconstruction schemes for prediction of near field flow region in complex aeroacoustics problems // Int. J. Aeroacoust. 2014. V. 13. № 3–4. P. 207–234. https://doi.org/10.1260/1475-472X.13.3-4.207

  12. Abalakin I., Bakhvalov P., Kozubskaya T. Edge-based reconstruction schemes for unstructured tetrahedral meshes // Int. J. Numer. Meth. Fluids. 2016. V. 81. № 6. P. 331–356. https://doi.org/10.1002/fld.4187

  13. Павлов Л.С. Распределение давления в сечениях прямоугольного крыла (лопасти) при криволинейном движении в несжимаемой среде // Уч. зап. ЦАГИ. 1979. Т. X. № 2. С. 104–108.

  14. Джонсон У. Теория вертолета. Т. 1. М.: Мир, 1983. 502 с.

  15. Spalart P.R., Allmaras S.R. A one-equation turbulence model for aerodynamic flows // 30th Aerospace Sciences Meeting and Exhibit, Aerospace Sciences Meetings. AIAA Paper 1992-0439. https://doi.org/10.2514/6.1992-439.

  16. Menter F.R. Two-equation eddy-viscosity turbulence models for engineering applications // AIAA Journal. 1994. V. 32. № 8. P. 1598–1605. https://doi.org/10.2514/3.12149

  17. Белов И., Исаев С. Моделирование турбулентных течений: учеб. пособ. Санкт-Петербург: Балт. гос. техн. ун-т, 2001. 108 с.

  18. Wilcox D.C. Turbulence Modeling for CFD, Third ed. La Canada, California: DCW Industries, 2006. 522 p.

  19. Reichard H. Vollständige darstellung der turbulenten geschwindigkeitsverteilung in glatten leitungen // ZAMM – Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 1951. V. 31. P. 208–219.

  20. Y. Saad. Iterative methods for sparse linear systems // Philadelphia: Society for Industrial and Applied Mathematics, 2003. 528 p.

  21. ANSYS, Inc. ANSYS, 1970–2018. URL www.ansys.com.

  22. Cucitore R., Quadrio M., Baron A. On the effectiveness and limitations of local criteria for the identification of a vortex // European Journal of Mechanics – B/Fluids. 1999. V. 18. № 2. P. 261–282. https://doi.org/10.1016/s0997-7546(99)80026-0

  23. ОВК НИЦ “Курчатовский институт”. URL http://computing.nrcki.ru/.

Дополнительные материалы отсутствуют.