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

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

_______ 1 ^ ^ 1

№ 3

УДК 629.735.33.015.3.025.1 : 532.526

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНОГО ОБТЕКАНИЯ ВЯЗКИМ ПОТОКОМ ТЕЛЕСНОГО КРЫЛА КОНЕЧНОГО РАЗМАХА

И. А. Гурьянов, В. Н. Котовский, М. И. Ништ

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

Общее решение задачи вязкого обтекания в широком диапазоне изменения параметров движения для сплошной среды обеспечивает полная система уравнений Навье—Стокса. И хотя этот подход в последнее время разрабатывается довольно успешно, однако для решения пространственной турбулентной задачи требуется преодоление ряда математических проблем и необходимо существенное повышение быстродействия ЭВМ [1]. Поэтому в настоящее время при изучении обтекания крыльев/ расчеты ограничены только некоторыми частными случаями [2].

В данной работе моделирование обтекания крыла вязким потоком строится согласно гипотезе Прандтля путем объединения нелинейного метода расчета невязкого течения с интегральным методом расчета пограничного слоя. Это позволяет проводить исследования обтекания крыла в широком диапазоне углов атаки, вплоть до закритических.

1. Постановка задачи. При моделировании обтекания крыла вязким потоком, как и в работе [3], течение делится на две области: область невязкого течения вне крыла и пограничного слоя и область вязкого течения в пограничном слое. Считается, что крыло из состояния покоя внезапно приводится в движение со скоростью ио, которая в дальнейшем не изменяется.

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

ности, а также начальные условия. Вихревая пелена, образующаяся за крылом, движется вместе со средой. На бесконечном удалении от крыла и его следа возмущения, вызванные ими, затухают. На кромках крыла, с которых сходят вихревые пелены, выполняется гипотеза Чаплыгина—Жуковского о конечности скоростей. Кроме того, поскольку движение производится из состояния покоя, в каждый расчетный, момент времени выполняется теорема Томпсона о равенстве нулю циркуляции по любому замкнутому контуру, проведенному через одни и те же частицы жидкости и охватывающему крыло и его след. Рассматриваются только положительные углы атаки, поэтому отрыв пограничного слоя моделируется только с верхней поверхности крыла.

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

Рассматривается обтекание крыла без скольжения, при этом моделируется только правая половина крыла (рис. 1), а влияние левой учитывается на основе условий симметрии. Расчетная схема характеризуется числом сечений на полукрыле N в которых выполяется условие непротекания поверхности, и числом поперечных вихревых отрезков п в сечении на каждой верхней и нижней поверхности крыла.

Поперечные вихревые отрезки характеризуются номером ^. Отсчет производится от задней кромки крыла вдоль нижней, а затем и верхней поверхности по задней кромки (1<^<2п) в сечении z = const. Нумерация продольных линий (сечений z = const на крыле) ведется слева направо. Если k —номер линии продольных вихревых отрезков, то k=1 соответствует корневой хорде, а k = N + 1 —торцу крыла. Сечения контрольных точек обозначаются номером р.

Суммарный вихревой слой на поверхности крыла, включающий присоединенные и свободные вихри, моделируется вихревыми отрезками с безразмерной циркуляцией г^:+1 г и продольными Г^*1 кг. Вне поверхности крыла существуют три системы свободных вихрей: кормовая I, сходящая с верхней поверхности II и боковая III. Кормовая система состоит из поперечных отрезков 8*'* + 1 и продольных Д'1) *, а система II — из поперечных 3^^+1 и продольных Д*2) '5. Боковая пелена моделируется одним поперечным вихревым отрезком а(3) * = 8$)я + 1 (s=1, 2, ... ,г).

Расчет потенциального обтекания крыла сводится к решению системы линейных алгебраических уравнений относительно величин циркуляции суммарных поперечных вихрей Г^:+17 на крыле и свободных поперечных вихрей ЗУ){+1 в системе I. Эта система получается в результате выполнения условия непротекания поверхности крыла в контрольных точках в расчетные моменты времени:

N 2п N )

£ г^1 г л^+1 р+£ 5**+1 а:: 1 Р=н;; :

Аа1 "'= 1 * = 1 I (2,1)

'1 = 1,2, ... ,2п+ 1; р = \, 2, ... , N.

Коэффициенты А левых частей этих уравнений выражаются через 'безразмерные функции т (х, у, г) [4]:

АЙ+‘р -К!+ ,-ИюХ‘.‘+1 -ЧГ1 *)? «?; I (2.2)

к = 1,2...... №, г = 1, 2______2п + 1. | ’

■Здесь й? —вектор нормали к поверхности крыла. Правые части имеют вид

г-1 ( , N N+1

НО V = - 4*/0* - £ (£ 8Я1 Н ®?£+1р + £ дГ +

^ = 1\б=1 А=2

+ £^* + 1 + £ д5?> 4 «#?_!*, + N+2 ®?,да+ь) , (2.3)

*=1 к=2 )

где /о, — закон изменения нормальной составляющей переносной скорости в контрольных точках.

Система уравнений (2.1) замыкается условием о постоянстве

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

2п 7-1

X ГЙ+1' + = - £ (44+1 + йЧ+0 ; (2.4)

",=1 *=1

к -:-=1, 2, ... , N.

Из анализа системы уравнений (2.1) и (2.4) видно, что она оказывается переопределенной, так как общее число контрольных точек для всего крыла на N больше, чем общее число неизвестных циркуляций. Устранение переопределенности этой системы осуществляется введением регуляризирующих переменных [5] с номером ^ = 2л+ 1.

Решение уравнений (2.1) и (2.4) осуществляется в каждый момент времени и проводится независимо. Последовательно вычисляются

коэффициенты левых частей (2.2) и правые части (2.3) уравнений решается система (2.1) и (2.4), определяются циркуляции Г'£:+1г и 3**+1. По этим значениям определяются циркуляции остальных вихрей, проводится расчет поля скоростей внешнего обтекания крыла и находятся аэродинамические нагрузки [6].

3. Расчет трехмерного пограничного слоя. Из решения задачи о внешнем обтекании крыла в каждый расчетный момент времени находится распределение скорости на верхней поверхности крыла. Затем на ней проводится выстраивание линий тока от линии растекания до-точек отрыва.

Интегральные соотношения получаются интегрированием уравнений пограничного слоя по координате п, отсчитываемой в направлении, нормальном к обтекаемой поверхности. Два интегральных уравнения импульса, одно и в основном направлении — по оси s, связанной с линией тока внешнего течения, и другое в поперечном направлении г, могут быть сведены к системе двух дифференциальных уравнений в частных производных первого порядка '[7], если задать профиль скорости в направлении основного потока в виде степенной зависимости, профиль скорости вторичного течения — в виде профиля Магера, а коэф-циент трения — по формуле, которая получена на основе экспериментальных данных и является обобщением формулы Людвига—Тилмана на трехмерные пограничные слои с градиентом давления.

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

В результате решения этих уравнений определяются толщина потери импульса 8и, формпараметр пограничного слоя Я, равный отношению толщины вытеснения к толщине потери импульса, и тангенс угла а№ предельной линии тока (eN = tgaw).

Граничные условия для расчета пограничного слоя задаются на линии растекания и в плоскости симметрии крыла. На линии растекания для определения 8ц и Н может быть использован профиль скорости в пограничном слое, полученный точным решением уравнений Навье—Стокса для течения в окрестности критической точки {8], а еш равно нулю, так как профиль скорости на линии растекания еще не успел деформироваться в поперечном направлении. Граничные условия в плоскости симметрии для решения уравнений трехмерного пограничного слоя получаются из расчета двумерного пограничного слоя вдоль корневой хорды, для чего используются уравнения, записанные для двумерного пограничного слоя.

Расчет осуществляется вдоль всех построенных линий тока до точки схода потока с задней или боковой кромок крыла или до наступления отрыва пограничного слоя. Отрыв трехмерного пограничного слоя существует двух типов {9]. К первому типу относится отрыв, когда становится предельно малым коэффициент поверхностного трения (С/—0) или предельно большой одна из производных дЯ/^ или дбц/дз, несмотря на то, что с^О. В численном решении это соответствует достижению значениям формпараметра Я = 2,2+-2,6 [10]. Другой тип отрыва имеет место при конечных С/, д8ц/д5, дН/дз в случае,, когда угол аи+Хе>90°—Хм (Хе — угол между направлением линии.

тока невязкого течения и осью Ох, /м — местный угол стреловидносте участка крыла), т. е. когда пограничный слой становится направленным вдоль размаха крыла.

4. Определение положения линии отрыва и моделирование оторвавшегося пограничного слоя. Отрыв пограничного слоя моделируется сходом. с поверхности крыла вихревой пелены П. Она состоит из дискретных вихревых отрезков, которые уходят в поток из узловых точек вихревой схемы крыла, ближайших к огибающей точек отрыва, полученных из расчета пограничного слоя. Интенсивность свободной вихревой пелены II считается равной потоку завихренности в сечении отрыва пограничного слоя [3, 6], т. е.

V = 5 V (- ) *,

о ' '

где V — модуль скорости в пограничном слое, V — скорость по коорди-

нате п, и — по координате s.

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

8(2)7==Ах 5 (-£ —4г )уЧп.

о 4 '

Так как в интегральных методах расчета пограничного слоя скорость V не вычисляется, и известно, что на большей части пограничного слоя она много меньше продольной скорости и, а ее производная ^/дя меньше производной ди/дп, то

о ди . ме 1 «5 и ~дп ЧП = -+ . о

Но при приближении к точке отрыва скорость V и производная ди/дя резко возрастают и в точке отрыва значение V становится соизмеримым с продольной скоростью и. Как показали расчеты пограничного слоя конечно-разностными методами [6], в сечении отрыва величина скорости V составляет 10—15%' от продольной скорости и. Следовательно, интенсивность у вихревой пелены получается несколько. меньшей, что учитывается введением коэффициента е, а '

Коэффициент е учитывает также наличие зоны обратного течения за местом отрыва и выбирается в диапазоне 0,7—0,8.

Достаточно большое значение скорости V в сечении отрыва приводит к движению свободного вихря 8(2) г не по касательной к поверхно-

• и

сти, " , что суще-

ствено влияет на параметры потока в окрестности точки отрыва и ее положение.

Чтобы определить скорость движения дискретного вихря 0(2) Г В момент отрыва, необходимо найти положение центра завихренности, который, как известно, отстоит от поверхности крыла на толщину вы-

правлении основного течения в пограничном слое и зная из расчета значение формпараметра Н, можно определить скорость оторвавшегося вихря 8(2)

Дальнейшее движение вихрей системы II так же, как и дискретных вихрей кормовой и боковой пелены, осуществляется по вектору местной скорости потока, вычисляемой в каждый расчетный момент времени r во всех узловых точках свободной вихревой пелены [4J.

5. Примеры расчета. С целью проверки разработанной методики расчета трехмерного турбулентного пограничного слоя было произведено сравнение с экспериментальными данными Ван ден Берга—Эль-сенаара и результатами расчета других авторов различными конечноразностными и нтегральными методами [11]. Из этого сопоставления можно сделать вывод о том, что описанный в п. 3 метод расчета пограничного слоя обладает удовлетворительной точностью как по вычислению его параметров, так и по определению положения линии отрыва на крыле.

На рис. 2, 3 приведены суммарные аэродинамические характеристики прямоугольного крыла Л = 4, полученные в эксперименте и в расчете по разработанной методике. Сплошной линией везде обозначены результаты расчета, кружочками — экспериментальные данные, полученные при M = 0,15, Re=10^105. Наблюдается хорошее согласование в зависимости коэффициента подъемной силы от угла атаки а (см. рис. 2). Моделирование отрыва потока с поверхности крыла позволяет оценивать как максимальный коэффициент подъемной силы сУа шах , так и критический угол атаки СХ;.:., а также выявить такое сложное явление как «выполаживание» зависимости коэффициента момента

теснения 87 пограничного слоя. Задавшись профилем скорости в на-

Рис. 2

тангажа тг от угла атаки. Качественно правильное согласование расчетных данных с экспериментальными наблюдается и по коэффициентам силы лобового сопротивления Сха и продольной силы Сх (рис. 3). Характерным для прямоугольного крыла является резкое падение подсасывающей силы и увеличение силы лобового сопротивления на закритических углах атаки. По этой причине у такого крыла обратная ветвь поляры второго рода идет выше прямой ветви и в расчете и в эксперименте. Обе расчетные поляры качественно правильно согласуются с экспериментальными во всем диапазоне углов атаки. Построение зависимости аэродинамического качества К от угла атаки п позволяет найти наивыгоднейший угол атаки (а = 8°) и оценить значение максимального аэродинамического качества Ктах=12—13.

Изменение коэффициента нормальной силы Су и циркуляций первых вихрей, сходящих с задней кромки 8Ш г и верхней поверхности

крыла 8\|) г в каждый расчетный момент времени г для первой расчетной полосы (р =1), иллюстрирует рис. 4. Здесь ". — безразмерное время:

где 1 — время, Ьо — хорда крыла. В конце переходного процесса нагрузка на крыло изменяется около какого-то среднего значения. Это происходит тогда, когда циркуляции 8(1) г и 8<2) г устанавливаются примерно равными по величине и противоположными по знаку. Это означает, что циркуляции присоединенных вихрей крыла тоже перестают изменяться, а следовательно, устанавливаются и коэффициенты аэродинамических сил. Изменение координат точек отрыва Хотр = Хотр/Ьо в трех различных сечениях по размаху крыла видно из рис. 5. Здесь zl=2z/I, где I — размах крыла. В первый момент времени, когда еще пелена, сходящая с верхней поверхности крыла, не образовалась, точка отрыва занимает наиболее удаленное от передней кромки положение. По мере накопления над верхней поверхностью крыла дискретных вихрей положительной циркуляции, моделирующих оторвавшуюся

пелену, скорость обтекания верхней поверхности в области сильного влияния пелены заметно падает, т. е. возрастает положительный градиент давления, что приводит к более раннему отрыву пограничного слоя. К моменту времени т^3,5 положение линии отрыва стабилизируется. В переходном процессе можно отметить моменты времени, в которые точки отрыва скачкообразно перемещаются вперед. Это связано с образованием и подходом к поверхности крыла вихревых сгустков в оторвавшейся пелене. Видно, что положение точки отрыва на ближайшей к торцу крыла линии тока (2; = = 0,9625) изменяется во времени в широких пределах. Это объясняется влиянием сложной вихревой структуры над концевой частью крыла, которая образуется вследствие взаимодействия вихревой пелены с поверхности и боковой пелены.

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

ЛИТЕРАТУРА

1. Х ирш е л ь Э., К о р д у л л а В. Сдвиговое течение сжимаемой жидкости. Численный расчет пограничного слоя. Пер. с англ. — М.: Мир,

1987 г.

2. Ш е н г Д. С. Обзор численных методов решения уравнений Навье—Стокса для течений сжимаемого газа. — Аэрокосмическая техника. — М.: Т. 52, 1986.

3. Б е л о ц е р к о в с к и й С. М., К о т о в с к и й В. Н., Н и ш т М. И.,

Ф е д о р о в Р. М. Математическое моделирование нестационарного обтекания кругового цилиндра.— Изв. АН СССР, МЖГ, 1983, N2 4.

4. Б е л о ц е р к о в с к и й С. М., Н и ш т М. И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. — М.: Наука,

1978.

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

6. Б е л о ц е р к о в с к и й С. М., К о т о в с к и й В. Н., Н и ш т М. И.,

Ф е д о р о в Р. М. Математическое моделирование плоскопараллельного отрывного обтекания тел. —М.: Наука, 1988.

7. Л а к ш м и н а р а я н а Б., Г о в и н д а н ного пограничного слоя на профилях плоских роторов турбомашин. — РТК, 1981, т. 19, N2 11.

8. Ш л и х т и н г Г. Теория пограничного

9. К и с е л е в А. Ф., Х о н ь к и н А. Д.,

Турбулентный пограничный слой несжимаемой ном крыле.— В сб.: Турбулентные течения и Труды ЦАГИ, 1985, вып. 2265.

10. Ч ж е н Р. Отрывные течения, т. 1. — М.:

11. Щек ин Г. А. Расчет трехмерного турбулентного пограничного слоя на крыле при наличии области отрыва. — В сб.: Современные задачи аэродинамики летательных аппаратов. — М.: МАИ, 1985.

Т. Р. Расчет турбулент-решеток и на лопатках

слоя. — М.: Наука, 1974.

В о р о т н и к о в П. П. жидкости на стреловид-сопротивление трения. —

Мир, 1972.

Рукопись поступила 14/У1 1989 г.