Обтекание телесного аэродинамического профиля потоком вязкой жидкости Текст научной статьи по специальности «Физика»

ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.

NATURAL SCIENCE. 2019. No. 4

УДК 532.5: 532.13 DOI 10.23683/0321-3005-2019-4-12-18

ОБТЕКАНИЕ ТЕЛЕСНОГО АЭРОДИНАМИЧЕСКОГО ПРОФИЛЯ ПОТОКОМ ВЯЗКОЙ ЖИДКОСТИ*

© 2019 г. Я.А. Бердник1, М.А. Сумбатян1

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

FLOW OF VISCOUS FLUID AROUND AN AIRFOIL PROFILE

Ya.A. Berdnik1, M.A. Sumbatyan1

Southern Federal University, Rostov-on-Don, Russia

Бердник Янина Александровна - аспирант, кафедра теоретической и компьютерной гидроаэродинамики, Институт математики, механики и компьютерных наук им. И.И. Воровича, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов н/Д, 344090, Россия, e-mail: yberdnik@sfedu.ru

СумбатянМежлум Альбертович - доктор физико-математических наук, профессор, кафедра теоретической и компьютерной гидроаэродинамики, Институт математики, механики и компьютерных наук им. И.И. Воровича, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов н/Д, 344090, Россия, e-mail: sumbat@math.sfedu.ru

Yanina A. Berdnik - Postgraduate, Department of Theoretical and Computational Hydroaerodynamics, Vorovich Institute of Mathematics, Mechanics and Computer Sciences, Southern Federal University, Milchakova St., 8a, Rostov-on-Don, 344090, Russia, e-mail: yberdnik@sfedu.ru

Mezhlum A. Sumbatyan - Doctor of Physics and Mathematics, Professor, Department of Theoretical and Computational Hydroaerodynamics, Vorovich Institute of Mathematics, Mechanics and Computer Sciences, Southern Federal University, Milchakova St., 8a, Rostov-on-Don, 344090, Russia, e-mail: sumbat@math.sfedu.ru

Рассматривается задача обтекания симметричного телесного аэродинамического пр офиля потоком вязкой жидкости с применением точных уравнений Навье - Стокса. Используется итерационный метод по малым возмущениям скоростей относительно скоростей основного потока. Путём применения преобразования Фурье к уравнениям системы дифференциальных уравнений она сводится к системе линейных алгебраических уравнений относительно образов Фурье. Для интегральных уравнений относительно функции вязкого трения на каждом шаге используется метод коллокации, сводящий интегральное уравнение, полученное из системы дифференциальных уравнений, к системе линейных алгебраических уравнений. Параметр сопротивления, в зависимости от различных профилей NACA с симметричным граничным контуром, сравнивается с решением П. Кунца. Итерационный метод, впервые предложенный авторами для обтекания тонкого телесного профиля потоком вязкой несжимаемой жидкости, сходится за несколько шагов. Функция профиля учитывается только на первом шаге итерации, а все последующие итерации работают в том же алгоритме, что и в задаче для тонкой пластинки, решённой авторами ранее. Точность метода проверена сравнением с имеющимися данными для Re = 2000. Максимальная относительная погрешность не превышает 1 % по сравнению с результатами П. Кунца.

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

This work presents a study on the problem of a flow of viscous fluid around a symmetric airfoil using the exact Navier -Stokes equations. An iterative method is used for small perturbations of velocities relatively to the velocities of the basic flow. By applying the Fourier transform to equations of a system of differential equations, it is reduced to a system of linear algebraic equations concerning Fourier images. A collocation method is used at each step for integral equations relatively to the viscous friction function, reducing the integral equation obtainedfrom the system ofdifferential equations to a system of linear algebraic equations. The drag coefficient, depending on different symmetrical NACA airfoils, is compared with Peter Kunz's solution. Thus, iterative method, which was first proposed by the authors for the flow of viscous incompressible fluid around a thin airfoil,

* Работа выполнена в рамках базовой части государственного задания Министерства образования и науки РФ № 9.5794.2017/8.9.

ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.

NATURAL SCIENCE.

2019. No. 4

converges in several steps. Whereas the airfoil function is taken into account only at the first iteration step, all following iterations use the same algorithm as the one for the problem of a thin plate solved by the authors earlier. The accuracy of the method is verified by comparison of the results obtained and the available data for Re = 2000. As a result, the maximum relative error does not exceed 1 % compared to the results of P. Kunz.

Keywords: flow around an airfoil, viscous fluid, drag coefficient, NACA airfoils, iterative method, integral equations, singular integrals, numerical methods.

Введение

Главные проблемы проектирования и изготовления микровоздушных транспортных средств (или микролетательных аппаратов) -сложность их конструирования, обеспечение летных качеств и возможность их полного контролирования. Многие годы преобладал в большей степени теоретический, нежели практический, интерес к изучению аэродинамики очень маленьких физических масштабов и небольших чисел Рейнольдса. Однако технологический прогресс с конца 90-х гг. XX в. сделал полёт при малых числах Рейнольдса реальным. Таким образом, растущий интерес к разработке микровоздушных транспортных средств создал потребность в улучшении понимания эффективной аэродинамики. Режим полёта микросамолёта или микровертолёта сопряжен с решением многочисленных проблем аэродинамического анализа и последующего конструирования. По этой причине существует относительно небольшое количество экспериментальных и теоретических работ для аэродинамических поверхностей, работающих при Re<10 000. Экспериментальные данные аэродинамического профиля были опубликованы в [1, 2] при числах Рейнольдса от 20 000 до 30 000, но с точки зрения природы вязких эффектов они значительно отличаются от диапазона ниже Re=10 000, представляющего значительный интерес.

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

Цель данного исследования - вычисление коэффициента сопротивления c¿ для различных симметричных профилей NACA и сравнение полученных данных с работой [3], где автор, в свою очередь, сравнивает свои результаты, полученные с использованием прямого численного метода, с экспериментальными [4, 5] и получает хорошее совпадение. Данная работа логически продолжает исследование [6], в котором предлагается новый итерационный полуаналитический метод изучения обтекания тонкой пластинки конечной длины однородным потоком вязкой несжимаемой жидкости, и использует результаты, полученные в этом исследовании.

Постановка задачи и вывод основных формул

Исследуется обтекание симметричного телесного аэродинамического профиля безграничным однородным потоком вязкой несжимаемой жидкости (рис. 1). Задача исследуется в двумерном приближении. Искомой величиной является коэффициент сопротивления c .

Рис. 1. Обтекание симметричного телесного аэродинамического профиля с хордой l = 2a набегающим потоком вязкой жидкости со скоростью U0 / Fig. 1. Flow around a symmetrical two-dimensional airfoil with a chord l = 2a of the incoming flow of viscous fluid at a velocity U0

ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKII REGION.

Течение жидкости рассматривается в декартовой системе координат Oxy, где горизонтальная ось Ox направлена вдоль набегающего потока, ось Oy - вертикально вверх. В силу симметрии течения рассматривается только верхняя полуплоскость. Течение жидкости полагается стационарным и ламинарным и описывается уравнениями Навье — Стокса и уравнением неразрывности:

NATURAL SCIENCE.

2019. No. 4

до до2 „ —1 + —2 = 0 , о

дх ду

до до

—1 + о —1

дх

— vAо +

1Ф

1 р дх

= 0,

до

до

2 + о —2 - vAo + — — = 0 .

ду

1—Р = 0 . (1)

дх ду р ду

В (1) р - массовая плотность жидкости; V -коэффициент кинематической вязкости; р - гидродинамическое давление; ц = ц и и2 = ц - продольная и поперечная компоненты вектора скорости соответственно.

Граничные условия прилипания на поверхности профиля имеют следующий вид:

и • Т = 0, (2)

и • п = 0, (3)

где т и п — соответственно единичные векторы касательной и нормали к граничной кривой профиля у =f( x).

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

сумм скоростей основного набегающего потока и малых на его фоне искомых возмущений (ц, ц):

их (х у) = и (х у) + и'х (х у), ^

иу(х у) = V(X у) + иу(х у).

При этом количество штрихов над функцией соответствует номеру итерации.

На нулевом шаге предполагается, что вектор скорости равен {и0,0}. Таким образом, на первом шаге: и(х,у) = и0, V(х, у) = 0. Также вводится функция тока соотношениями их = ,

ду

т =

1 f '(х)

Jl + f '2(х) Vi + f '2(х)f

f '(х) _1

Vi + f '2 (х)

U = \Uо +

ду' ду'

ду дх

(5)

(6)

(7)

Штрих над функцией ^ x) означает производную. Предполагаем, что телесный профиль очень тонкий, следовательно, величина /'(х) мала.

В первом приближении, подставляя разложения (4), записанные в терминах функции тока, в исходные уравнения Навье — Стокса, придем к линейному уравнению в частных производных с постоянными коэффициентами:

U

дАу'

дх

-vA2y' = 0 .

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

У

' (а, у) = С,е

х\у

+ С е

va —U^ia

, (8)

где а - параметр преобразования Фурье.

В линейном приближении сносим граничное условие (3) на ось x, т.е. берем его при у = 0. При

этом оно принимает вид

U0 +

ду

ду

f + У = 0 дх

(9)

= . Выражение для функции тока преду дх

ставляется аналогично разложению скоростей на первой итерации: у(х, у ) = ^(у ~) + у' (х, у) ,

т(у) = и0у .

Векторы касательной и нормали к поверхности крыла т и п , а также вектор скорости потока и имеют вид

Условие (9) выполняется не только в области профиля, но и вне его на оси x, если считать, что в этой области / = 0. В самом деле, последнее слагаемое в (9) соответствует вертикальной компоненте вектора скорости, которая вне профиля в силу симметрии равна нулю.

Используя принцип замораживания функции / (х), применяя преобразования Фурье к соотношению (9) и учитывая (8), получаем

_ /а(С1 + С2) + /'(х)_ а\С\ _ _ =

V У

= _Ц)/' (а). (10)

V

ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIIREGION.

NATURAL SCIENCE.

2019. No. 4

Введём в рассмотрение вектор напряжений Коши р^ = РцП] . Так как т • п = 0, то после некоторых математических преобразований получаем выражение для рт:

Рт= (11)

= М

дщ

( (

2 ■

V 1^1

Чп1 +-

дщ

дхо

1 (дщ ди2 V 2 ) I дх2 дх1 /

Отбрасывая функцию /' в силу малости и принимая во внимание (10), получаем из выражения (11) уравнение

4ia(ia(C1 + C2 J - U0f '(a jJ+

+ f2a2C1 +(2a2 -"¿H)^ 1 = 1

) ) М

где Т' = р'т - касательная сила трения на поверхности профиля.

Запишем соотношения (10), (12) в виде системы линейных алгебраических уравнений (СЛАУ) относительно неизвестных коэффициентов С1 и С2:

(13)

Решаем эту СЛАУ методом Крамера. Главный определитель системы (13)

А = -2a'

- ia- f 'Ja -

2 U0ia

(ia + f '|a| J 2a2 +

2 , U0ia

v

..2 rt

= 2а2 / >р - 2а2 / '|а| + ^ _ /' |«|.

V V V V

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

2

лости по сравнению с а :

А*Uo ,

v

0 (a2 - ia | a | f 'J.

Аналогично, оставляя только члены главного порядка малости по а и отбрасывая квадраты величин типа /', получаем

А1 =

f ff ^ — + 4iaU 0f' М

V

2 U 0ia - ia - f Ja--

- U of'( 2a2 + ^

A T А1 * — М

-aa- f ' a2 -^

f ff

U o2 f'

ia

(12) А2 = 2a2Uof' + - + 4iaUof' |(ia + f '|a\J

М

IL

М

(ia + f '|a\J.

Соотношение (2) с учетом выражений (5)-(7) принимает вид

U • Т = o

- 2a2Cj 2a2 + u°1h^|c2 =— + 4HUo f '(aJ,

( I-Ü~iä) ~ _____"

(- ia-f '(x|HJc1 + - ia-f '(xJJa2 - C2 =-U of '(aJ- f'(XJ+®

U + f' = o.

ду

дх

Отсюда следует

X

2л

J (C1 (aJ + C2 (afJe~iaxiada -

(14)

of

\

e~1Hda = -U,

0 :

71 НС^Ь а2 -^С,(а) 2л-х[ V V )

|х| < а.

Первый интеграл в (14) можно отбросить, так как он содержит / '2( х) - малое 2-го порядка малости по /' . Выражение, входящее во второй интеграл, связано с величиной

На + Л -

2 Uoia

А 2 =

(15)

2 Urda

Н - Ja2--0—

= — (- га) М

Используя (15), получаем

2

Uof'

ia\a\.

T_

2л ^

= -Un

(- ia)

Unia

2 U(V H -Ja--0

v

U

0 rt

f ia\a\

U-°(a2 - ia\a\f J

-da =

v

v

v

V

V

V

V

2

V

ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIIREGION. i——^

NATURAL SCIENCE.

2019. No. 4

Г

1 и

или--J —

V

1—. 11+Uo-

via

+ -

U,

j

v

f '

2n.

U,

■ = —U,

о

if'+isigrUa))

о

В результате, применяя формулу обращения

для Т' , приходим к интегральному уравнению относительно силы трения на граничном контуре профиля:

v

jfTSei 4 j

1 —

1 +

U

о

и —00 f' ix) + isigria)

2U0n—a

TT +a Jai^—x)

= —U0 + — J f' i§)d4—A-

2n—a f' ix) + isigria)

Va eiai4—x)da =

da .

При этом интеграл в правой части уравнения вычисляется в явном виде

+o e'ai4—x)

i +

f '(*)

da + J

+o e~iai4—x)

— i + f' ix)

da =

=7 тО^. 2 Г <&*+2*í sOyf,

O 1 + f' ix) O 1 + f' 2 ix)

2f' ix)Si? — x)

+

1

2

1 + /2 (х) х 1 + /2 (х)'

Отбрасывая первый член в последнем выражении как содержащий нелинейный множитель / 2 , перепишем правую часть интегрального уравнения в виде _ Щ + и° • 2 | /.

2* _а £,_х

В результате основное интегральное уравнение принимает следующий вид:

v

J ^ d#j-

1 —

1 + Uo-

via ia e

nUо —a и —о isigria) + f' ix)

= —uo + J

Uo a f4

d4.

4x)da =

(16)

* _а х

Любопытно, что при /(х) = 0 уравнение (16) вырождается в интегральное уравнение обтекания тонкой пластинки вязким потоком несжимаемой жидкости [6], как и должно быть с физической точки зрения.

Уравнение (16) решается численно методом коллокации [7, 8] - для нахождения первого при-

ближения к решению для телесного симметричного профиля.

Полученное первое приближение подставляется в виде «входного» для итерационного процесса, использованного ранее для тонкой пластинки [6]. Таким образом, толщинная функция профиля учитывается только на первом шаге итерации, а последующие итерации работают в том же алгоритме, что и для тонкой пластинки.

Этот метод тестируется на профилях серии NACA00t [9], граничный контур которых симметричен и описывается уравнением

у = /(х) = ± 5г (0,2969 4х _ 0,1260 х _

_ 0,3516 х2 + 0,2843 х3 _ 0,1015 х4),

где 0 < x < 1 - горизонтальная координата; t максимальная толщина профиля (совпадает с номером из двух последних цифр номера NACA, делённым на 100. Например, толщина профиля NACA0002 равна 02/100=0,02).

Отметим, что сравнение идёт по коэффициенту сопротивления Cd [10] (сравни также с [6]):

О =

F,

1pUо2/'

W = SWaz = SU¡ 0,664

0/2W,

Fd = 2W,

Fd = 2SU0/2 W,

где W - полная сила граничного трения с одной

ч

стороны профиля, WBlaz = 0,664 pU0.

2a v

- ана-

логичная величина, полученная Г. Блазиусом приближенно в аналитическом виде [10].

В итоге для коэффициента лобового сопротив-

4S

ления получаем Cd = — P '

2aUn

(таблица).

Сравнение полученных результатов с результатами Kunz [3] / Comparison of the results obtained with the results of Kunz [3]

Профиль NACA Cd [3] Cd

NACA0000 (пластинка) 0,0656 0,0667

NACA0002 0,0682 0,0690

NACA0004 0,0708 0,0709

NACA0006 0,0734 0,0727

NACA0008 0,0760 0,0755

v

0

0

V

ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIIREGION.

Таким образом, максимальное отклонение результатов, представленных в таблице, приблизительно составляет 1 %. Графическое сравнение представлено на рис. 2.

006

0 NACA NACA NACA NACA

0002 0004 0006 0008

Рис. 2. Изменение параметра сопротивления Cd в зависимости от номера профиля NACA.

Сплошная линия соответствует расчётам Kunz [3], пунктирная - методу, предложенному авторами / Fig. 2. Change of drag parameter Cd versus NACA profiles. The solid line corresponds to Kunz calculations [3], the dashed line corresponds to the method proposed by the authors

Заключение

В работе предлагается новый итерационный метод расчета ламинарного обтекания тонкого симметричного телесного профиля потоком вязкой жидкости при небольших числах Рейнольдса. Метод состоит в последовательной линеаризации задачи на фоне приближения предыдущего шага. На каждом шаге решается линейная краевая задача для уравнения в частных производных 4-го порядка с переменными коэффициентами. Применяя принцип замораживания и преобразование Фурье вдоль потока, приходим к некоторым интегральным уравнениям относительно функции вязкого трения на контуре профиля. Метод тестируется на семействе профилей класса NACA00t [9].

ЕСТЕСТВЕННЫЕНАУКИ. 2019. № 4 NATURAL SCIENCE. 2019. No. 4

Максимальное отклонение от экспериментальных данных для Re = 2000 по силе сопротивления составляет не более 1 %.

Литература

1. Schmitz F. W. Aerodynamics of the Model Airplane. Huntsville, AL: NASA Marshall Space Flight Center, 1967. 201 p.

2. Althaus D. Profilpolarenfur den Modellflug. Band 2. GmbH: Neckar-Verlag, 1988. 176 p.

3. Kunz P.J. Aerodynamics and Design for Ultra-Low Reynolds Number Flight: A Dissertation for the Degree of Doctor of Philosophy. Stanford, 2003. 180 p.

4. Sunada S., Kawachi K. Comparison of Wing Characteristics at an Ultra Low Reynolds Number // J. of Aircraft. 2002. Vol. 39 (2). P. 331-338.

5. Azuma A., Okamato M., Yasuda K. Aerodynamic Characteristics of Wings at Low Reynolds Number. Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications // Progress in Aeronautics and Astronautics, AIAA. Reston, VA, 2001. Vol. 195, ch. 17.

6. Berdnik Y., Beskopylny A. The approximation method in the problem on a flow of viscous fluid around a thin plate // Aircraft Engineering and Aerospace Technology. 2019. Vol. 91 (6). P. 807-813.

7. Калиткин Н.Н. Численные методы. М.: Наука, 1978. 512 с.

8. Sumbatyan M.A., Scalia A. Equations of mathematical diffraction Theory. Boca Raton, FL: CRC Press /Taylor&FrancisGroup, 2005. 291 p.

9. NACA airfoil. URL: https://en.wikipedia.org / wiki/NACA_airfoil (дата обращения: 27.08.2019).

10. Schlichting H. Boundary-layer theory. N.Y.: McGraw-Hill, 1955. 535 p.

References

1. Schmitz F.W. Aerodynamics of the Model Airplane. Huntsville, AL: NASA Marshall Space Flight Center, 1967, 201 p.

2. Althaus D. Profilpolarenfur den Modellflug. Band 2. GmbH: Neckar-Verlag, 1, 1988, 176 p.

3. Kunz P.J. Aerodynamics and Design for Ultra-Low Reynolds Number Flight: A Dissertation for the Degree of Doctor of Philosophy. Stanford, USA, 2003, 180 p.

ISSN 0321-3005 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIIREGION.

4. Sunada S., Kawachi K. Comparison of Wing Characteristics at an Ultralow Reynolds Number. Journal of Aircraft. 2002, vol. 39 (2), pp. 331-338.

5. Azuma A., Okamato M., Yasuda K. Aerodynamics Characteristics of Wings at Low Reynolds Number. Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications. Progress in Aeronautics and Astronautics, AIAA, Reston, VA, 2001, vol. 195, ch. 17.

6. Berdnik Y., Beskopylny A. The approximation method in the problem on a flow of viscous fluid around a

ECTECTBEHHblEHAyKM. 2019. № 4 NATURAL SCIENCE. 2019. No. 4

thin plate. Aircraft Engineering and Aerospace Technology. 2019, vol. 91 (6), pp. 807-813.

7. Kalitkin N.N. Chislennye metody [Numerical methods]. Moscow: Nauka, 1978, 512 p.

8. Sumbatyan M.A., Scalia A. Equations of Mathematical Diffraction Theory. Boca Raton, FL: CRC Press /Taylor&FrancisGroup, 2005, 291 p.

9. NACA airfoil. 2006. Available at: https://en.wi-kipedia.org/wiki/NACA_airfoil (accessed 27.08.2019).

10. Schlichting H. Boundary-Layer Theory. New York: McGraw-Hill, 1955, 535 p.

Поступила в редакцию /Received_12 сентября 2019 г. / September 12, 2019