Расчет безотрывного околозвукового обтекания стреловидных крыльев с учетом влияния вязкости Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XIX 19 8 8

№ I

УДК 629.735.33.015.3.025.1 : 532.526

РАСЧЕТ БЕЗОТРЫВНОГО ОКОЛОЗВУКОВОГО ОБТЕКАНИЯ СТРЕЛОВИДНЫХ КРЫЛЬЕВ С УЧЕТОМ влияния вязкости

Н. А. Владимирова, В. В. Вышинский, Г. А. Щекин

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

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

При расчетном исследовании пространственного околозвукового обтекания крыльев важным моментом является учет вязкости. Экспериментальные исследования вязко-невязкого взаимодействия при околозвуковых скоростях показывают, что из-за турбулентного характера течения в пограничном слое и малой интенсивности скачков уплотнения в большинстве случаев крейсерского полета пассажирских самолетов не происходит отрыва пограничного слоя при его взаимодействии со скачками уплотнения. При дозвуковых скоростях, а также в случае сверхкритического обтекания со слабыми скачками уплотнения вязконевязкое взаимодействие является слабым и может быть учтено в соответствии с идеей Прандтля посредством изменения поверхности обтекаемого тела на толщину вытеснения пограничного слоя. Подобная схема расчета реализована при исследовании влияния вязкости в случае безотрывного дозвукового и околозвукового обтекания профиля [1, 2| и тела вращения [3, 4].

При создании метода расчета пространственного околозвукового обтекания крыльев с учетом влияния вязкости авторы учли опыт создания подобных методов в двумерном [1—4] и пространственном [5—8] случаях. В частности, результаты работ [5, 6] позволяют сделать вывод, что использование двумерного метода расчета пограничного слоя в со-

четании с пространственным методом расчета внешнего потенциального обтекания обеспечивает хорошее согласование с экспериментальными результатами в распределении давления лишь в случае изолированных крыльев большого удлинения малой и умеренной стреловидности, однако в случае крыльев малого удлинения и/или большой стреловидности хорошего согласования с экспериментом не наблюдается. В работах [7,8] для расчета течения в пограничном слое на крыле используются пространственные интегральные методы, что позволяет надлежащим образом отразить течение у конца крыла, где из-за поперечных перетеканий в пограничном слое нарушается гипотеза плоских сечений. Сравнение с экспериментальными эпюрами давления, которое проводится в работе [8], показывает, что применение пространственного метода расчета пограничного слоя дает наиболее значительное улучшение (по сравнению с методом, использующим «плоскую» вязкую поправку [6]) на нижней поверхности сверхкритического крыла в области задней кромки.

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

Для расчета внешнего потенциального обтекания прямого и стреловидных крыльев в данной работе используется конечно-разностный метод [10], где в качестве исходного уравнения применяется полное уравнение относительно потенциала скоростей. Как показывают, в частности, результаты работы [5], методы, в которых решается полное уравнение относительно потенциала при точном выполнении условий непро-текания на теле, применимы к решению более широкого класса задач и дают результаты, которые обычно лучше согласуются с экспериментом, чем методы малых возмущений.

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

1. Расчет внешнего потенциального обтекания крыла проводится по методу, изложенному в работе [10]. Рассматривается стационарное изоэнтропийное течение невязкого политропического газа. Программа расчета основана на численном решении полного уравнения для потенциала скоростей Ф, которое в декартовой системе координат х, у, г, где ось х направлена вдоль хорды корневого ,профиля, ось г — в направлении размаха крыла, а ось у— вверх, имеет вид:

(аа - и2) Фхх + (а2 — V2) Фуу + (а3 - да2) Фгг — 2иг> Фху — 2ит Фхг —

— 2^1) Фуг = 0. - (1)

2а

Здесь и=Фх, V = Фу, т — Фг— составляющие скорости в направлении осей координат, а — ух/э/р — местная скорость звука, которая выражается посредством уравнения Бернулли:

а» - + тА- - ^ («2 +^2 + ®2), (2)

М00

Мао — число Маха набегающего потока.

Течение предполагается безотрывным, и на поверхности крыла или эффективного (с учетом толщины вытеснения пограничного слоя) тела выполняется условие непротекания:

д Ф дп

= 0,

(3)

где п—вектор внешней нормали к обтекаемой поверхности. Граничные условия на бесконечности имеют вид:

дФ

дх

Т Т ЙФ I I т •

X2 у3 + г2 оо СОБ “> 1*3 + уз + г2 _ со ^ 81П а;

^1 -0

дг I*2 + у2 + г2 ->■ оо ’

(4)

где а — угол атаки крыла.

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

С целью экономии машинного времени в программе использованы алгоритм двукратного удвоения расчетных сеток и алгоритм «замораживания» дальнего поля на сетке после второго удвоения [11], что позволяет сократить время расчетов в три — пять раз.

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

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

котором исходные уравнения в криволинейной неортогональной системе координат г], £, где ось г) направлена по нормали к поверхности в данной точке, а единичные векторы | и ? принадлежат касательной плоскости и составляют между собой угол 0 (0<6<я), записываются в виде:

где и, V, ш—проекции вектора скорости в пограничном слое на координатные оси |, г], £,

В уравнения (5) входят коэффициенты, зависящие от формы поверхности: hu h2, klt k2, ki2, &2i- В частности, hi и h2— метрические коэффициенты, ki и k2■—значения кривизны линий | = const и £=const. Поверхность крыла предполагается гладкой и достаточно плавной, так что радиус кривизны поверхности в каждой точке много больше толщины пограничного слоя.

Система уравнений замыкается уравнением состояния p = pRT и эмпирической формулой Сазерленда:

Связь между скоростью и температурой устанавливается через интеграл энергии срТ0 = ср Г + и2(/2.

При ламинарном течении в пограничном слое и' г>' = г/ Н' =

— р' V' = 0. Для расчетов течения в турбулентном слое по аналогии с коэффициентами молекулярной вязкости ц и теплопроводности К вводятся коэффициенты турбулентной вязкости [Хг и турбулентной теплопроводности /(<. В случае изотропной модели турбулентности коэффициент турбулентной вязкости является скалярной функцией. При этом

^ (р uh2 sin 6) + (р vt hx Л2 sin 0) 4- ^ (p rwhl sin 6) = 0;

V! = v + p’v'lp, u2t = u? -(- w2 + 2 ww cos 6, H = cpT 4- u)j2.

I T\

T \ 3/2 Tm + 110,4 H T+ 110,4 •

проекции вектора турбулентного касательного напряжения трения на оси і и £ определяются посредством:

ди

тге = — ря'г»' = рі^,

—, дт

хп-=

Величина турбулентного теплового потока определяется как:

„ _ 7ПТ_Кі ді _ ді

Чі— р дт] —

где і—энтальпия: і = срТ + сопэ!, а Рг* = ^ р ср1К{ ~ 0,9—местное турбулентное число Прандтля.

Для определения коэффициента турбулентной вязкости использует» ся одна из двух моделей турбулентности: двухслойная модель Себеси [12] или усовершенствованная единая модель Сполдинга [13]. В переходной области коэффициент турбулентной вязкости вычисляется по формуле где продольный коэффициент перемежаемости у записы-

вается в виде:

1_ехр[- 4,605(1-4)’],

| = ^ соответствует началу переходной области, а % — — ее концу.

Граничные условия для уравнений движения при отсутствии вдува» отсоса и турбулентности внешнего невязкого потока имеют вид:

т] = 0, и = V = да = и’ V' = V' да' — V' Н' = р' V' = 0, т] = 8, и = ие, да = тюе, и! V' = •г/да' = г>'/У' =■- р' V' = 0,

Для уравнения энергии на внешней границе г| = 6, Н=Не, на поверхности крыла т] = 0, либо Н = Ту,/Т0—при заданной температуре поверхности 7\е(|, £), либо =0 — при отсутствии теплообмена меж-

О Т£)

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

Обратное влияние пограничного слоя на внешнее невязкое течение-характер изуется толщиной вытеснения б*, величина которой определяется из решения уравнения баланса расхода через произвольный замкнутый объем. При отсутствии вдува и отсоса это уравнение имеет вид;

||(Р,и,А, МП в• («* — *п)| + Д [р, V,Н, гт в- (8* _ ?,,,)] - 0, где 8п_|(1-^)<гч, «„ = ](!-^1-),/,.

о о

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

3. Метод реализован в виде программы расчета на ЭВМ типа БЭСМ-6. Время выполнения одной большой итерации (расчет потен-

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

Для проверки точности получаемых результатов был проведен контрольный расчет обтекания прямоугольного крыла большого удлинения (1=20), составленного из 12%-ного симметричного профиля, координаты которого приведены в работе [10]. Центральное сечение г=0 такого крыла обтекается практически плоским потоком. Поэтому к этому сечению в полной мере применима гипотеза плоских сечений, и характеристики течения, рассчитанные для пространственного и плоского случаев, Должны быть близки. На рис. 1 и 2 приведены результаты расчета толщины вытеснения б* и коэффициента трения с/ в сечении 2 = 0 для режима закритического обтекания Моо = 0,75, угла атаки а = 4°, числа Ке = 8,7-106 при фиксированном положении перехода хп=30% на

Прямое крыло X =20, г=0;Ис-6,7!0е;х„-30 %

Рис. 1

■-----пронстранстбенный пограничный слой (модель Себеси)

° ". " ( - Сполдинга)

-----дйу мерный • - .

Рис. 2

верхней и нижней поверхностях крыла. Сплошными линиями нанесены результаты расчета по предлагаемому методу, штриховыми линиями — по двумерному дифференциальному методу Алберса — Грегга, штрих-пунктирными— по двумерному интегральному методу Грина. Описание последних методов и их сравнительный анализ могут быть найдены в работе [4]. При расчете пространственного пограничного слоя использовались две модели турбулентности: модель Себеси [12] — эти результаты нанесены сплошной линией, и усовершенствованная модель Сполдинга [13] — эти результаты нанесены кружочками.

Как видно, лучшее согласование с двумерными методами в распределениях б*(х) на турбулентном участке дает модель Сполдинга (см. рис. 1), которая используется в дальнейших расчетах. В зависимостях с^(х) смена моделей почти не проявляется. Применение двух двумерных методов позволяет установить приемлемую при сравнении разбежку результатов. На нижней поверхности удается получить особенно хорошее соответствие в значениях б* и Cf между пространственным и двумерными методами на всем протяжении профиля, особенно если учесть, что фиксирование точки перехода производится с точностью до ближайшего расчетного узла сетки.

По разработанной методике проведены расчеты вязкого обтекания крыла с углом стреловидности % = 30°, удлинением к = 8 и сужением г] = 2,5, составленного из симметричного профиля, координаты которого приведены в [10], с максимальной относительной толщиной с= 12%. на режиме докритического обтекания Моо = 0,5 и су = 0,46, при числе Рейнольдса по корневой хорде крыла Ке = 3-106. Линия перехода ламинарного пограничного слоя в турбулентный фиксировалась в расчете на расстоянии хп=3% на верхней поверхности и хп = 5% на нижней поверхности от передней кромки. (Значения абсциссы х отнесены к местной хорде крыла).

Распределение давления на верхней поверхности носит пиковый характер, свойственный профилям с передним положением максимальной толщины профиля. На рис. 3 в трех сечениях по размаху г = 0; 0,3; 0,6 (координата г отнесена к полуразмаху крыла) приведены распределения трех компонент скорости и, V, ъV на верхней и нижней поверхностях крыла. Значения компонент скорости отнесены к скорости набегающего потока. Как видно, максимальные значения поперечной компоненты

Верхняя поверхность

скорости ы), характеризующей степень пространственности обтекания, составляют в среднем 0,1-ь0,15, достигая величины 0,5 на линии растекания и = и = 0, расположенной вблизи передней кромки крыла.

На рис. 4 приведены распределения вдоль хорд местных коэффициентов трения ср Для сравнения в тех же сечениях приведены результаты расчета для двумерного пограничного слоя с использованием гипотезы плоских сечений по методу Алберса—Грегга. Как видно, расчет плоского пограничного слоя (штриховые линии) на верхней поверхности дает заниженные значения величины Cf, причем это отличие увеличивается по мере продвижения от корня (г = 0) к концу крыла. На нижней поверхности расчет с/ по плоскому методу дает близкие, но слегка завышенные (в сравнении с расчетом пространственного пограничного-слоя) результаты. Следует отметить, что в окрестности точки перехода ламинарного течения в турбулентное наблюдается скачок в значении местного коэффициента трения, обусловленный изменением закона трения при переходе ламинарного пограничного слоя в турбулентный.

На рис. 5 приведены зависимости от координаты г вдоль размаха коэффициентов сопротивления трения сечений крыла (интеграл по-хорде) для верхней Ср в и нижней Ср н поверхностей, а также суммарный коэффициент трения Ср в + Ср н. Последняя величина слабо зависит от 2, слегка увеличиваясь к концу крыла. Здесь же штрихпунктирной линией нанесена величина полного коэффициента сопротивления трения крыла Ске =0,0071, что составляет 83% от величины индуктивного-сопротивления Сжг = 0,0086.

Влияние вязкости посредством толщины вытеснения на распределение нагрузки по рассмативаемому крылу невелико. Коэффициент полной подъемной силы су, рассчитанный с учетом вязкости (три вязких итерации) отличается лишь на 4% от значения су, полученного при расчете потенциального обтекания.

ЛИТЕРАТУРА

I. Брутян М. А., Серебрийский Я. М. Приближенный метод расчета подъемной силы и продольного момента профиля с учетом вязкости при малых скоростях. — Ученые записки ЦАГИ, 1976, т. 7, № 1.

2. Melnik R. Е., Chow R., Mead Н. R. Theory of viscous transonic flow over airfoils at high Reynolds numbers.—AIAA Paper 77-680, 1977.

3. M a r i n g D. E. A theoretical study of body drag in subcritical axisym-metric flow.— Aeron. Quart, 1976, vol. 27, N 3.

4. Вышинский В. В. Метод расчета околозвукового безотрывного обтекания тел вращения с учетом вязкости. — Труды ЦАГИ, 1981, вып. 2109.

5. Hinson В. L. and Burges К. P. An evaluation of three-dimensional transonic codes using new correlation-tailored test data.—AIAA Paper 80-0003, 1980.

6. Mason W. H., et al. An automated procedure for computing the three-dimensional transonic flow over wing-body combinations including viscous effects. —AFFDL-TR-77-122, 1977, vol. I.

7. Street C. L. Viscous-inviscid interaction for transoniq wing-body configurations including wake effects. — AIAA J., 1982, vol. 20, N 7.

8. Waggoner E.G. Transonic three-dimensional viscous-inviscid interaction for wing-body configuration analysis. — AIAA Paper 82-0163,

1982.

9. Щ e к и н Г. А. Численный расчет трехмерного пограничного слоя в ламинарной и турбулентной областях течения на крыле при сверхзвуковых скоростях полета. — В сб.: Экспериментальное и теоретическое исследования аэродинамических характеристик JIA и его частей.— М.: МАИ,

1983.

10. В л а д и м и р о в а Н. А. Исследование обтекания прямых и стреловидных крыльев большого удлинения при околозвуковых скоростях. —■ Ученые записки ЦАГИ, 1983, т. 14, № 4.

II. Владимирова И. А., Вышинский В. В., Ляпунов С. В., Серебрийский Я. М. Об ускорении сходимости методов расчета плоского и пространственного трансзвукового обтекания тел в неограниченном потоке. — Ученые записки ЦАГИ, 1985, т. 16, № 4.

12. Cebeci Т. Calculation of three-dimensional boundary layers.

I. Swept infinite cylinders and small cross flow. — AIIAA J., 1974, vol. 12, N 6.

13. Колина H. П., Солодкин E. E. Программа на языке ФОРТРАН для численного интегрирования уравнений пространственного пограничного слоя на линии растекания и на бесконечном скользящем цилиндре.— Труды ЦАГИ, 1980, вып. 2046.

Рукопись поступила 26/V 1986 г.