МЕТОД УЛУЧШЕНИЯ АНАЛИТИЧЕСКИХ РЕЗУЛЬТАТОВ ДЛЯ ВРАЩАЮЩЕЙСЯ ЛОПАСТИ В СРАВНЕНИИ С ПРЯМЫМ ЧИСЛЕННЫМ РАСЧЕТОМ Текст научной статьи по специальности «Физика»

ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION NATURAL SCIENCE. 2021. No.4

Научная статья УДК 533.6

doi: 10.18522/1026-2237-2021-4-19-25

МЕТОД УЛУЧШЕНИЯ АНАЛИТИЧЕСКИХ РЕЗУЛЬТАТОВ ДЛЯ ВРАЩАЮЩЕЙСЯ ЛОПАСТИ В СРАВНЕНИИ С ПРЯМЫМ ЧИСЛЕННЫМ РАСЧЕТОМ

Евгений Алексеевич Казаков1, Межлум Альбертович Сумбатянш

1,2Южный федеральный университет, Ростов-на-Дону, Россия

1Eugene.A.Kazako V@yandex. ru

2masumbatyan@sfedu.ruB

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

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

Благодарности. Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (РФФИ), проект № 19-29-06013.

Для цитирования: Казаков Е.А., Сумбатян М.А. Метод улучшения аналитических результатов для вращающейся лопасти в сравнении с прямым численным расчетом // Известия вузов. Северо-Кавказский регион. Естественные науки. 2021. № 4. С. 19-25.

Статья опубликована на условиях лицензии Creative Commons Attribution 4.0 International (CC-BY 4.0).

Original article

METHOD FOR IMPROVING ANALYTICAL RESULTS FOR A ROTATING BLADE VERSUS DIRECT NUMERICAL CALCULATION

Eugene A. Kazakov1, Mezhlum A. Sumbatyan2b

1,2Southern Federal University, Rostov-on-Don, Russia 1Eugene.A.Kazako V@yandex. ru 2 masumbatyan@sfedu.ru M

Annotation. The paper proposes an analytical method for calculating the value of the driving force of a thin rotating plate. The classical approach represents the plate as a set of plane problems corresponding to the section of the blade. The theory of a thin wing describes the translational motion of the plate. Unlike them, the proposed method allows us to consider a rotating thin plate as a whole. A two-dimensional integral equation permitting to describe a plate of any shape and twist is obtained. The singular terms of the obtained equation qualitatively correspond to the singular terms of the classical equation for the translational motion of the thin plate, and additional regular terms are responsible for the rotational character of the motion. A numerical experiment of rotating thin plate was performed in the AN-

© Казаков Е.А., Сумбатян М.А., 2021

SYS CFX simulation software. The results of the numerical experiment are compared with the data obtained by the classical method of the plane problem and the proposed analytical method. The proposed method permits reducing the deviation of the draft value obtained analytically from the value obtained numerically by three times - from 200 to 75 %.

Keywords: aerodynamics, thin plate, rotation, theory of thin lifting surface, perturbation theory, analytical method

Acknowledgments. This work was financially supported by the Russian Foundation for Basic Research (RFBR), project No. 19-29-06013.

For citation: Kazakov E.A., Sumbatyan M.A. Method for Improving Analytical Results for a Rotating BladeVersus Direct Numerical Calculation. Bulletin of Higher Educational Institutions. North Caucasus Region. Natural Science. 2021;(4):19-25. (In Russ.).

This is an open access article distributed under the terms of Creative Commons Attribution 4.0 International License (CC-BY 4.0).

Введение

Стандартный способ аналитических расчетов для вращающейся лопасти - применение двумерной теории к набору плоских задач, соответствующих сечениям лопасти. Поскольку двумерная теория обтекания тонкой дуги дает явные аналитические результаты [1], такой метод традиционно используется в теоретических работах и при проектировании, несмотря на то что он дает существенную ошибку. Как правило, грубые теории дают сильно завышенное значение силы тяги по сравнению с истинным значением. Более точные теории описаны классиками гидро- и аэродинамики [2], однако не дают достаточную в сравнении с натурными экспериментами точность, тем более что эти теории непросто реализовать на практике. В связи с этим для этой задачи развивались различные приближенные подходы [3, 4]. Между тем вопрос о возможности построить точную с теоретической точки зрения теорию, аналогичную теории малых возмущений в теории тонкого крыла, и надежную с практической точки зрения в смысле близости с экспериментом, по-прежнему остается открытым. Некоторый вариант подхода к этой задаче именно в данной трактовке реализуется в настоящей работе.

Основной целью работы является построение расчетной модели в рамках классической теории тонкого крыла и теории малых возмущений, которая позволяет с большей точностью, чем двумерная, рассчитать значение силы тяги для вращающейся лопасти. Предлагается новый полуаналитический метод, который сокращает погрешность получаемых результатов до 75-80 % вместо 200 %, которые наблюдаются при применении двумерной задачи.

Постановка задачи

Рассматривается тонкая слабоизогнутая лопасть, вращающаяся вокруг оси x1 c угловой скоростью ш, так что в текущий момент времени ось лопасти проходит через ось Х2. Рабочая поверхность лопасти описывается уравнением х1 = f(x2,Хз) (рис. 1). В соответствии с классической трактовкой принимается линеаризованная аэродинамическая теория по малым возмущениям. Предполагается, что малыми являются закрученность лопасти и угол атаки а, т.е. величины |3//3х2| и |3//3хз|.

Рис. 1. Привязка координатных осей к тонкой лопасти / Fig. 1. Arrangement of the coordinate axes of a thin blade

Таким образом, поверхность лопасти располагается близко к плоскости вращения Х2Х3, и угол между вектором нормали к поверхности лопасти и осью x1 также мал. Это позволяет перенести граничное условие непроницаемости поверхности лопасти на соответствующую ей проекцию S в плоскости вращения Х2Х3. Предполагается, что течение жидкости безвихревое везде, за исключением вихревой пелены, срывающейся с задней кромки лопасти. Задача является стационарной во вращающейся системе координат, связанной с лопастью.

Метод решения

Для применения стандартных методов теории потенциала запишем кинематическое условие рассматриваемой задачи:

р(хх, г, Ы + e,t) = (p(xi, г, а>д£ + д£ = 0*д£ = —шдЁ-

Интеграл Коши - Лагранжа в линеаризованном виде

р -t р р -t р дв' ( )

Граничное условие непроницаемости на лопасти имеет вид (v ■ п) = (щ ■ п), где верхний штрих означает возмущение соответствующей величины, нижний индекс b - лопасть (blade). С учетом того,

что в линейном приближении v = { vi, v2, v3}, п = {l, — j^, —a}, vb = {0,0, wr] « {0,0, wx2], получаем

(v ■ п) = (щ ■ п) * vi = —wx2a, (2)

где отброшены малые второго порядка малости.

Основная формула теории потенциала [5]:

р'(х) = JJS

any any

dsy (3)

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

ТТ А Г Г 1 1 (х = (х1,х2,хз) = (х1,г,в)

на. При этом функция Грина имеет вид и = — = —;-- , { , . , ...

4пг 4п\х-У\ (у = (У1,У2,Уз) = (У1,Р,Ч>)

С учетом соотношения (1)

&пу У) И Хдфдуц' дфду2' дфду3) \ ' ду2 ) И дфду1

Из (3)следует

дср'(х)

= Л • ^^ - -д^— и (у, х) (Бу = Л \у(у) -д^^С(у,х)] (1зу, (4)

дв дтр дпу дтрду-1 у ду± дтр у' 4 '

где у(у)- функция, связанная со скачком аэродинамического давления на лопасти; х = (х^^,х2,хз) -произвольная фиксированная точка в жидкости; у = (у1, у2, Уз) - точка, по которой происходит интегрирование по поверхности. В явном виде у (у) выглядит следующим образом:

д(р'

дф

д(р'

дф

рш

Очевидно, что в цилиндрической системе координат

(У2= ^cosxp ,(X2=rCOSd,

ly3 = ^sin-ф, {x3 = rsin0. ( '

Из (2) и (5) следует, что v'(y) = —шу2а, поэтому

„ dv' 1 , s ду2 Maud(cos-) , .

G—t =--(-Ы ■ а)-fr=----= -Cúa(-sin-é).

дф 4ЛЦ v y д-ф 4лц д-ф v

Эта функция нечётна по поэтому по крайней мере для лопастей, симметричных по X3, а следовательно, и по ^, второй интеграл в (4) равен нулю. Тогда (4) принимает вид

дср(х)

дв 4п JS ' (y ) ду1

Атг

ttSY(.y)

Л

_^(У1-Х1)2 + (У2-Х2)2+(Уз-Хз)2_

dy2 dy3.

У 1=0

Это эквивалентно равенству $ ^ = — JJ '(y)

Xidy2 dy3

дв 4TTJJs'K^ [Х12+(у2-Х2)2+(Уз-Хз)2]3/2

1

ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION NATURAL SCIENCE. 2021. No.4

Отсюда

dV

дх1дв

= ±ff

Y(y)dyz ¿Уз

X!

=0 4nJJS[(y2-X2)2+(y3-X3)2]3l2 ■

(6)

Интегрируем (6) по в. Поскольку d9 =

Х2

х2+х:

■dx3, то интегрирование по в можно заменить на инте-

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

д2<р'

дх1 дх3

*2

х!

=0 4n(x2+x2)JJS

ИЧУ(У2'Уз)

dy2 йуз

[(Х2-У2)2 + (Хз-Уз)2]3'2

„'I —д<Р' ^l\xi=0 =

дх1

х1=0

Р œ

= — —

4п I

J%3

dz

И

Г(У2,Уз)^У2 dУз

x\+Z2"S [(Х2-У2)2+(г-уз)2]3/2'

Последнее выражение с учетом граничного условия (2) принимает вид

r œ

4п

J*3

dz

x\+z2

И

У(У2,Уз)^У2 ¿Уз Ч(х2-У2)2 + (г-Уз)2]3/2

= ша.

(7)

Вычислим интеграл (7) по частям:

I

œ

dz

dv =

(г2+х2)[(г-уз)2+(х2-у2)2] <2 dz

U =

Z2+X-

du = —■

2z dz (z2+x2)

Тогда

I =

[(г-уз)2+(Х2-У2)2]3'2'

Хз-Уз

V =

z-Уз

(*2 -У2)2 V(г-Уз)2+(*2-У2)2

(х2+Х2)(Х2-У2)2^(Хз-Уз)2 + (Х2-У2)2 + (*2-У2)2 Jx (z2 +X2f ^г—ЬУ+ХЪ-^

œ

г(г-уз) dz

(8)

Обозначим последний интеграл через Поскольку вдоль хорды лопасти величины х3,у3 малы, если лопасть узкая, то положим в (8) х3 ~ 0, у3 ~ 0. Получим

Г™ 9

I г2 dz I dz 9 | dz

I г . -><2

и

(z2+X2) Vz2+(X2-V2)2

= f'

0

(г2+х1Уг2 + (х2-У2)2

œ

-*2 J

0

(Z2+X2) Vz2+(X2-V2)2'

Обозначим первый интеграл в последней строке ¡2 и воспользуемся табличным значением для второго интеграла [6]. В результате получим

где

h = ¡2

h =

xl

zVz2+(X2-yz)2

2[(x2-y2)2-x2]x2(z2+x2)

+ ■

(*2-У2)2-2Х2

z=0

ln

Jx2-(X2-y2)2+X2

\Х2-У2\

h =

X2l(X2-y2)2-xi

■- arctg

X2

2[(Х2-У2)2-Х!]Х2 2) 2[(Х2-У2)2-Х2]

при \Х2-У2\<Х2 , при \Х2-У2\ > Х2.

[(Х2-У2)212-1]

(9)

Окончательно из (7) получаем следующее двумерное интегральное уравнение:

ИУ(У2,Уз)К(х2, х3,у2,y-i)dy2dy3 = ша, (Х2, х3) Е S,

где ядро имеет вид

К = I =

Уз-Хз

+ ■

(Х2-У2)2Ь-1

(Х2 +Х2)(Х2 -У2^(Хз-Уз)2+(Х2-У2)2 (х2 -У2)2 [(х2 -У2)2 -Х^] '

(10)

При этом ¡2 выражается формулами (9). Поскольку

= 1[_

(Х2-У2)2[(Х2-У2)2-Х^] X2 [(Х2-У2)2-Х2 (Х2-у2)2

то из (10) выделяется главная сингулярная часть

К =

Уз-Хз

+ ■

(х2 +Х2 ) (Х2 -y2)2V (Хз -Уз)2+ (Х2 -У 2 ) 2 Х2(Х2-У2)2 (Х2-У2)2-Х| X2 [(Х2-У2)2-Х2]

l

2

2

СО

l

l

Х2\Х2-(Х2-У2)2

(Х2-У2)2-Х2

l

l

l

l

l

I

l

2

где первые два слагаемых - сингулярные члены, третье и четвертое - регулярные члены ядра. В итоге первый член является аналогом первого члена для классической теории тонкого крыла, второй - аналогом второго [7], третий и четвёртый члены, не имеющие особенностей, отвечают за вращательный характер движения пластинки как целого вместо поступательного в классической теории крыла.

Интегральное уравнение (7) решалось методом дискретных вихрей [8] в разработанной компьютерной программе в среде программирования С++.

Решение прямым численным методом

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

Для случая вращательного движения тонкой пластинки классическую формулу двумерной теории [1] необходимо модифицировать с учётом изменения скорости набегающего потока вдоль размаха лопасти Р = 2npv2ba, где b - полухорда. С учётом того, что v(%2) = имеем

Р = j^ 2пр(шх2)2Ьа dx2, Р = 2прЬаш2(Ь3 — 13).

В качестве геометрических параметров исследуемого объекта были выбраны линейные размеры лопасти квадрокоптера DJI Mavic Air: длина хорды 2b = 1 см; длина по размаху L — I = 6 см; расстояние от оси вращения до лопасти I = 0,5 см. Вычислительный эксперимент рассматривал цилиндрический объём воздуха, вращающийся вокруг собственной оси с угловой скоростью, равной угловой скорости вращающейся лопасти. Цилиндр был разделён на два центрально-симметричных субдомена плоскостью, содержащей в себе жёсткую тонкую пластинку (рис. 2), соответствующую размерам лопасти. Такой подход позволяет избежать технических ошибок при описании идеально тонкой пластинки, а также даёт возможность рассматривать силу на каждой её стороне независимо. Угол атаки пластинки задавался наклоном разделительной плоскости относительно оси, проходящей через хорду.

Рис. 2. Геометрия модели в среде ANSYS CFX / Fig. 2. Geometry of the model in the ANSYS CFX software

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

ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION .NATURAL SCIENCE. 2021. No.4

В данной работе используется представленная в ANSYS CFX модель турбулентности k-epsilon, либо ламинарный расчёт без учёта турбулентности, система уравнений Навье - Стокса, включающая в себя уравнения неразрывности, моментов, полной либо тепловой энергии и уравнение состояния идеального газа. С учётом рассматриваемых параметров исследуемой лопасти сжимаемость среды не оказывает существенного (более одного процента) влияния на итоговое значение силы тяги. Характерный размер элемента в численном эксперименте - порядка сантиметра, вблизи лопасти, разделительной плоскости и боковой поверхности цилиндрического домена - порядка миллиметра. Общее число элементов в сетке - 1,4 млн.

Результаты расчётов

В первой серии экспериментов (табл. 1) варьировалось значение угла атаки а при неизменной угловой скорости ш = 1257 рад/с, плотность невозмущённого идеального газа р= 1,1765 кг/м3.

Таблица 1 / Table 1

Значения силы тяги для случая ш = 1257 рад/с / The value of the draft force, ш = 1257 rad/s

Р, H C++/ANSYS Классика/ANSYS

а,0 ANSYS CFX C++ Классика

1 0,02872166 0,05428 0,093513 1,89 3,255835491

2 0,0599567 0,1086 0,187026 1,81 3,119351132

3 0,0718167 0,1628 0,280539 2,26 3,906319839

4 0,1125846 0,2171 0,374052 1,93 3,322408216

5 0,1602236 0,2714 0,467564 1,69 2,91819682

8 0,2317718 0,4342 0,748103 1,87 3,227756785

10 0,265236 0,5428 0,935129 2,05 3,525648856

Во второй серии расчётов (табл. 2) при фиксированном угле атаки а=5° варьировалось значение угловой скорости, прочие параметры оставались неизменными.

Таблица 2 / Table 2

Значения силы тяги для случая «=5° / The value of the draft force for «=5°

P, H C++/ANSYS Классика/ANSYS

œ, рад/с ANSYS CFX C++ Классика

628,5 0,0354135 0,0678 0,116891 1,92 3,30074689

1257 0,1602236 0,2714 0,467564 1,69 2,91819682

1885,5 0,362547 0,611 1,05202 1,69 2,901747911

2514 0,649678 1,0862 1,870258 1,67 2,878746087

3142 0,980484 1,6967 2,921348 1,73 2,979495841

Заключение

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

ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION NATURAL SCIENCE. 2021. No.4

Список источников

1. Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1973. 848 с.

2. Жуковский Н.Е. Вихревая теория гребного винта. М.; Л.: Гостехиздат, 1950. 239 c.

3. Mescheryakov K., Sumbatyan M.A., Bondarchuk A.A. A boundary integral equation over the thin rotating blade of a wind turbine // Engineering Analysis with Boundary Elements. 2016. Vol. 71. P. 20-26.

4. Sumbatyan M.A., Mescheryakov K.I. Two-dimensional integral equation for a thin wind turbine blade rotating in the round tunnel // Mechanics Research Communications. 2017. Vol. 85. P. 1-4.

5. Сумбатян М.А., Скалия А. Основы теории дифракции с приложениями в механике и акустике. М.: Физ-матлит, 2013. 328 c.

6. Прудников А.П., Брычков Ю.А., Маричев О.И. Интегралы и ряды. Элементарные функции. М.: Наука, 1981. 800 c.

7. Sumbatyan M.A., Bondarchuk A.A., Mescheryakov K.I. An efficient numerical algorithm in the classical 3D theory of thin lifting surface in a flow of non-viscous incompressible fluid // Mechanics Research Communications. 2018. Vol. 89. P. 18-22.

8. Белоцерковский С.М., Лифанов И.К. Численные методы в сингулярных интегральных уравнениях. М.: Наука, 1985. 254 c.

References

1. Loitsyanskiy L.G. Mechanics of liquids and gases. Moscow: Nauka Publ.; 1973. 848 p. (In Russ.).

2. Zhukovsky N.E. Vortex theory of screw propeller. Moscow, Leningrad: Gostekhizdat Publ.; 1950. 239 p. (In Russ.).

3. Mescheryakov K., Sumbatyan M.A., Bondarchuk A.A. A boundary integral equation over the thin rotating blade of a wind turbine. Engineering Analysis with Boundary Elements. 2016;71:20-26.

4. Sumbatyan M.A., Mescheryakov K.I. Two-dimensional integral equation for a thin wind turbine blade rotating in the round tunnel. Mechanics Research Communications. 2017;85:1-4.

5. Sumbatyan M.A., Scalia A. Fundamentals of diffraction theory with applications in mechanics and acoustics. Moscow: Fizmatlit Publ.; 2013. 328 p. (In Russ.).

6. Prudnikov A.P., Brychkov Yu.A., Marichev O.I. Integrals and series. Elementary functions. Moscow: Nauka Publ.; 1981. 800 p. (In Russ.).

7. Sumbatyan M.A., Bondarchuk A.A., Mescheryakov K.I. An efficient numerical algorithm in the classical 3D theory of thin lifting surface in a flow of non-viscous incompressible fluid. Mechanics Research Communications. 2018;89:18-22.

8. Belotserkovsky M., Lifanov I.K. Numerical methods in singular integral equations. Moscow: Nauka Publ.; 1985. 254 p. (In Russ.).

Информация об авторах

Казаков Е.А. - аспирант, кафедра теоретической и компьютерной гидроаэродинамики, Институт математики, механики и компьютерных наук им. И.И. Воровича.

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

Information about the authors

Kazakov E.A. - Postgraduate, Department of Theoretical and Computational Hydroaerodynamics, Vorovich Institute of Mathematics, Mechanics and Computer Science.

Sumbatyan M.A. - Doctor of Science (Physics and Mathematics), Professor, Department of Theoretical and Computational Hydroaerodynamics, Vorovich Institute of Mathematics, Mechanics and Computer Science.

Статья поступила в редакцию 02.09.2021; одобрена после рецензирования 17.09.2021; принята к публикации 26.11.2021. The article was submitted 02.09.2021; approved after reviewing 17.09.2021; acceptedfor publication 26.11.2021.