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

УЧЕНЫЕ ЗАПИСКИ Ц АГ И Том XIII 1982

№ 1

УДК 533.6.013.2.011.35:629.7.025.73

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

Ю. П. Нуштаев

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

Задача об автоколебаниях органов управления самолета в трансзвуковом потоке известна уже приблизительно 30 лет. В 50-х годах с появлением трансзвуковых аэродинамических труб получили развитие экспериментальные методы исследования этого явления [1]. Было показано, в частности, что потеря устойчивости колебаний элерона при трансзвуковых скоростях может быть трех видов. Два из них наблюдаются в случае, когда скачки уплотнения размещаются за осью вращения органа управления или не-' посредственно на задней кромке профиля. Однако дальнейшее исследование, особенно исследование флаттера в системе с несколькими степенями свободы, тормозилось из-за отсутствия теоретических методов. Разработанные в последнее время численные методы расчета позволяют рассчитать необходимые для решения задачи о флаттере нестационарные силы, действующие на профиль при наличии на его поверхности колеблющихся скачков уплотнения, а также изучить особенности нестационарных трансзвуковых течений. Первой отечественной работой в этой области была работа [2], где для решения нестационарных уравнений Эйлера использован конечно-разностный метод Годунова. Однако при этом потребовались очень большие затраты машинного времени и поэтому полученные численные результаты служат в ос-

новном для проверки более приближенных быстродействующих методов расчета.

С целью упрощения рассматриваются безвихревые изоэнтропи-ческие течения, что позволяет свести задачу к решению уравнения для потенциала скорости возмущения. При этом исследуются течения только со слабыми скачками уплотнения. В данной работе решается уравнение для потенциала малых возмущений, получаемое в предположении, что угол атаки и толщина профиля, а также амплитуды колебательного движения профиля малы. Это упрощение приводит к уменьшению точности расчета в области резкого изменения наклона поверхности, например, в области носка профиля или элерона. Для низкочастотного диапазона колебательного движения профиля задача сводится к решению более простого нелинейного уравнения Линя — Рейснера — Тзяна [3]. Интерес к движениям с малой частотой обусловлен тем, что автоколебания органов управления обычно характеризуются относительно небольшими приведенными частотами (^^0,1). В данной работе на основе численного интегрирования по времени 'уравнения Линя—Рейснера— Тзяна исследованы аэродинамические характеристики профиля NACA 64А—006 при его колебаниях относительно оси, проходящей через середину хорды и при колебаниях элерона, занимающего 1/4 хорды данного профиля. Рассчитанные при этом нестационарные аэродинамические коэффициенты могут быть использованы для решения модельной задачи о крутильноэлеронном флаттере. Полученные результаты дополняют результаты работ [2, 4].

1. Рассмотрим уравнение Линя — Рейснера — Тзяна для потенциала малых возмущений, которое описывает двумерные нестационарные трансзвуковые течения идеального газа [3]:

2В<рх/ = С<рхх 4- <оууу (1)

где

Ш* 1—м* ,

Я—6 = --^-(Т4-1)М1<р„

•е — относительная толщина профиля, ?— показатель адиабаты,

Л Л

принятый в расчетах равным 1,4; приведенная частота

00

<о — физическая частота колебательного движения профиля;, Ь — длина полухорды профиля. Уравнение (1) записано в безразмерных величинах; х, у, <р, Ь— отнесены соответственно к величинам

Ь, й/т1-'3, дт2'3 К», Уравнение (1) выводится в рамках приближения:

А — — (1 _М1)« 1. (2)

Уравнение характеристик для уравнения (1) имеет вид

Сі2 — В2 у2 + 25^=0.

Имея в виду, что в уравнении (1), записанном в размерных вели-

М2 *

чинах, В > С= 1—М«, и заменяя М<» на V^ja, получим

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

того же порядка, что и характерное время распространения

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

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

где /^(х, *) — форма поверхности профиля. В обоих выражениях в низкочастотном приближении отброшены члены, пропорциональные малой частоте к (величина коэффициента давления в звуковых точках ср выбрана равной 0). Таким образом, скачок потенциала через вихревую пелену, как и в стационарном случае, принимается постоянным вдоль оси ;с. В силу этого, в качестве наиболее простого граничного условия на большом удалении вниз по течению от профиля можно использовать следующее: ¥*(+00, _у) = 0, а для остальных границ потока, например, <р=0. В качестве начального решения уравнения (1) в работе использовано значение потенциала возмущения, полученное при решении итерационным методом релаксации стационарного уравнения малых возмущений, совпадающего с уравнением (1) при £ — 0. Для уменьшения влияния границ расчетного поля на результаты расчета в данной работе использовалось преобразование координат 1 = Г1~Г1(У)> переводящее бесконечную физическую область течения в область внутри квадрата |^| <11. При этом неравномерная

в физической плоскости (х, у) расчетная сетка с концентрацией узлов в направлении у в окрестности профиля, а в направлении х — в окрестности кромок профиля, преобразуется в равномерную

или для трансзвукового потока, где V^^a:

у2 — 2at [х + (а — 1/<») /].

(х, 0, t)

ср = ~2 Нг'3?х

расчетную сетку в плоскости (5, 73). Использованное преобразование координат имеет вид:

*

1

е«>- 1

е-*г * _ 1

аі

1

при 1!>.х;>0,

при — 1 -<л: <0,

1 —(1 - У при *>1,

(1 — £0) £“»(!+■*) — 1 при х < 1,

Ч = + (1 — е^а»*) при у^О,

где 5= + ^о, Т1 = 0 координаты кромок профиля в новой системе. Параметры аи а2, а3 определяют степень концентрации узлов расчетной сетки при фиксированном количестве узлов всей сетки. При этом параметры а2 и ах связаны между собой условием гладкости функции £ (х) в точке Е = £0. Уравнение (1), записанное в дивергентной форме для системы координат (?, -г]), принимает вид:

4т\

(Т + О м'

іІт,

Чу

МІ — і

(7+1) г2/3 М

Лх

?1

=0.(3)

Граничные условия записываются аналогично. На луче расчетной сетки тг} = 0 вводятся два слоя расчетных точек (/ — у'0 4- 0), в которых разностная аппроксимация производной ущ записывается с учетом граничного условия на теле и следе. Использующаяся при этом величина циркуляции Г рассчитывается по формуле:

Г = ^ 1ТЕ> уо_0 — (Р/ге, /0+о,

где индекс 1ТЕ — обозначает номер расчетного узла на профиле, ближайшего к задней кромке. Эта запись представляет собой разностную аппроксимацию условия Жуковского для острой задней кромки профиля.

Численное интегрирование уравнения (3) в данной работе проведено с помощью неявной конечно-разностной схемы переменных направлений, описанной в [4]. Составлена программа расчета на ЭВМ для расчетной сетки, имеющей 81 и 42 узла, соответственно в направлении € и % Шаг интегрирования для всех расчетов был выбран равным 3°, что соответствует 120 шагам по времени в течение одного периода. При этом для расчета одного периода требуется 13 минут машинного времени на ЭВМ БЭСМ-6. Результаты расчета сходятся за два—три периода.

2. Рассмотрим симметричный профиль ЫАСА 64А—006 с относительной толщиной 6%. Он имеет тупую переднюю кромку и поэтому его кромки при расчетах располагаются между узлами сетки. Ближайший к кромке расчетный узел на поверхности профиля располагался на расстоянии 1% хорды, а количество расчетных точек на всей поверхности профиля равнялось 50. Угол атаки профиля изменялся по закону а = а05т©^, где амплитуда колебания а0 принималась равной 1/2° и 3/2°. Считалось, что элерон

колеблется относительно своего носка по закону S == £0 sin с амплитудой 50, равной 1° и 3/2°. Все расчеты проведены для величины 6 = 0,05 в диапазоне изменения чисел М« от 0,8 до 0,95.

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

в виде: тг = \тг\ sin (<ot — Ф) или тг = am* + гдеа = -^-, Ф —фаза запаздывания изменения величины т2 относительно движения профиля.

Для иллюстрации точности используемого численного алгоритма на рис. 1 проведено сравнение численного решения линейного нестационарного уравнения, получающегося из уравнения (1)

при 7 = — 1 с данными точной линейной теории для пластины при колебаниях относительно середины с частотой k = 0,05. Причинами расхождения сравниваемых результатов является численная погрешность алгоритма и погрешность низкочастотного приближения (2), в рамках которого в нестационарном уравнении малых возмущений отброшен член, пропорциональный <р/Л а в граничном условии и в выражении, для коэффициента ср — члены, пропорциональные <Dt. Как видно из рис. 1, при числе Мс* 1 расхождение сравниваемых результатов уменьшается, что связано с уменьшением погрешности приближения: (1—<С 1.

На рис. 2 и 3 представлены графики распределения модуля величины \&ср\ и фазы запаздывания Ф коэффициента распределенной нагрузки Аср = срн, п — срв. п в зависимости отсечения профиля для а0=1/2°. На рис. 2 представлена как 1-я, так и 3*я гармоника разложения в ряд Фурье величины Дср (пунктирные кривые), дающая для приведенных вариантов расчета наибольший вклад из всех гармоник высокого порядка. Рассмотрим подробнее результаты, представленные на рис. 2. При числе Моо = 0,8, соответствующем бесскачковому обтеканию данного профиля, „пик“ распределения нагрузки, наблюдаемый в области тупой передней кромки, вызван резким изменением в этой области наклона поверхности. Появление скачка уплотнения на поверхности профиля (Моо = 0,85) приводит к появлению дополнительного „пика" нагрузки в области перемещения скачка. При увеличении числа Мсо этот „пик“ растет и смещается к задней кромке. При этом величина [Дс^1 в точках, располагающихся вверх по течению от области колебания скачков уплотнения, заметно уменьшается. Это показывает,

Рис. 2 Рис. 3

что чем интенсивнее скачок уплотнения, тем он сильнее препятствует распространению возмущений вверх по течению от области его колебания. Сравнение 1-й и 3-й гармоник разложения величины Дср показывает, что при колебании развитых сверхзвуковых зон вдоль поверхности профиля нагрузка в большинстве точек профиля остается линейной от угла а, т. е. изменяется по закону, близкому к синусоидальному. Исключением являются точки профиля, через которые проходит скачок уплотнения. В этих точках вклад 1-й и З‘й гармоник становится одного порядка и линейность нарушается. Расчеты показали, что по мере приближения области колебания скачков уплотнения к задней кромке профиля (Мсс^0,9) вклад гармоник высокого порядка в распределение величины Дср уменьшается. При увеличении чисел Моо^0,9 наблюдается слабое изменение величины |Дср| на большей части поверхности профиля, что является отражением закона стабилизации стационарных распределенных характеристик профиля по числам М,*,. Графики рис. 2 показывают, что при дальнейшем увеличении чисел М*» (Моо^0,91) по мере стабилизации скачков уплотнения на задней кромке профиля их интенсивность уменьшается. Аналогичные результаты получены при колебании элерона. В этом случае вблизи оси шарнира элерона наблюдается дополнительный „пик* распределения нагрузки, обусловленный резким изменением в этой области наклона поверхности п рофиля с отклоненным элероном. Результаты, представленные на рис. 2, позволяют судить о вкладе

гармоник высокого порядка в величину модуля суммарных аэродинамических коэффициентов. Если относительный вклад этих гармоник в величину коэффициента Су мал во всем диапазоне чисел Моо, то их вклад в величину тг становится определяющим для числа Моо, при котором величина J/nJ, обусловленная 1*й гармоникой разложения, близка к нулю. При увеличении или уменьшении числа Моо относительный вклад гармоник высокого порядка резко уменьшается и величина тг становится практически линейной функцией от угла а.

На рис. 3 представлены распределенные фазовые характеристики колеблющегося профиля. Если при бесскачковом обтекании профиля (Мсо = 0,8) во всех точках профиля изменение величины &ср запаздывает по фазе относительно изменения угла а(Ф>0), то при появлении скачков уплотнения, в области их перемещения но верхней и нижней поверхностям профиля наблюдается сильное изменение и смена знака величины Ф. При этом в точках профиля, расположенных вниз по течению от области колебания .скачков, величина Ф<0, т. е. изменение величины Аср, опережает по фазе изменение угла а. При увеличении числа М» рост интенсивности скачков уплотнения сопровождается более интенсивным изменением величины Ф в области колебания скачков и увеличением величины отрицательного фазового запаздывания. По мере стабилизации скачков уплотнения в районе задней кромки профиля интенсивность изменения величины Ф уменьшается. При этом величина Ф вновь становится положительной на всей поверхности профиля.

На рис. 4 показан график перемещения скачка уплотнения по верхней поверхности профиля для амплитуды а0 == 1/2°. Расчеты показали, что на поверхности данного профиля ■ при числе Mr» — 0,85 в некоторый момент времени (см. рис. 4) возникает скачок уплотнения. Перемещаясь затем вниз по течению, он достигает крайнего к задней кромке положения не в момент отклонения профиля на максимальный угол, как это имеет место в стационарном течении, а в момент времени, когда угол а уменьшается. Таким образом, скачок уплотнения перемещается с фазовым запаздыванием относительно движения профиля. При перемещении вверх по течению скачок достигает своей максимальной интенсивности. Следовательно, интенсивность скачка в свою очередь изменяется с фазовым запаздыванием относительно его перемещения. В результате последующего уменьшения интенсивности скачок уплотнения переходит в звуковую линию, после чего образовавшаяся бесскачковая местная сверхзвуковая зона постепенно сокращается и исчезает. Такая закономерность возникновения и исчезновения скачка уплотнения имеет периодический характер. При увеличении числа Моо скачок уплотнения присутствует на протяжении относительно большей части периода, а при числе Моо^0,88 для данного профиля — в течение всего периода •и колеблется по закону, близкому к синусоидальному. В силу симметричности рассматриваемого профиля и колебания его относительно нулевого угла атаки перемещение скачка уплотнения на нижней поверхности совпадает с перемещением скачка на верхней поверхности с фазовым сдвигом, составляющим 180°. При увеличении числа Моо, как уже отмечалось, область колебания скачков уплотнения на поверхности профиля смещается к задней

si'acvtt

Рис. 4

Рис. 5

Рис. 6

з^омке, а амплитуда колебания скачков уменьшается (см. рис. 4). При числе Моо^0,92 скачок уплотнения практически стабилизируется в области задней кромки профиля. При этом наблюдается уменьшение фазы запаздывания перемещения и интенсивности скачка практически до нуля. Аналогичные результаты получены хля профиля с колеблющимся элероном. Рассмотренные типы движения скачка уплотнения соответствуют двум типам движения сыачка, наблюдавшимся экспериментально на поверхности профиля : колеблющимся элероном [5]. Расчеты показали, что увеличение амплитуды колебательного движения профиля приводит к увели* чению амплитуды колебания скачков уплотнения и смещению области их колебания к задней кромке профиля, что влечет за со-, бой соответствующие изменения в распределении нестационарной нагрузки Саср. При увеличении амплитуд а0 или 50 фаза запаздывания перемещения и интенсивности скачка уплотнения относительно движения профиля увеличивается. По мере приближения области колебания скачков к задней кромке влияние амплитуды или 80 на движение скачков, а, следовательно, и на распределен |Дср!

ление величин ; Ф уменьшается.

На рис. 5 и б представлены результаты расчета нестационарных аэродинамических коэффициентов (/я^бшФ и [ тщ( sin Ф. Величины этих коэффициентов пропорциональны работе, совершаемой силами давления, действующими на профиль и элерон, за период колебания. Положительный знак этой работы или коэффициентов аэродинамического демпфирования |m£|sln<D и (т^эшФ означает, что колебания системы с одной из рассмотренных степеней свободы будут затухать. Характер изменения величины аэродинамического демпфирования в зависимости от перемещения по поверхности профиля скачков уплотнения можно выяснить из рассмотрения распределения работы сил давления А (х) вдоль профиля (см. рис. 5). Площади заштрихованных областей обозначают отрицательный вклад в величину демпфирования профиля. Положение и ширина этих областей соответствуют области колебания скачков уплотнения на поверхности профиля. Таким образом, отрицательное демпфирование- формируется на той части профиля, вдоль которой колеблются скачки уплотнения, при условии, что эта область располагается за осью вращения. Расчет величины демпфирования данного профиля |m*|sin® относительно оси вращения, проходящей через середину хорды (см. рис. 5), показал, что при увеличении числа Моо величина аэродинамического демпфирования растет и достигает своего максимального значения при возникновении на поверхности профиля скачка уплотнения. При этом скачок колеблется вблизи оси вращения и его вклад в величину отрицательного демпфирования невелик. При дальнейшем увеличении числа М,» область колебания скачков уплотнения на нижней и верхней поверхности смещается к задней кромке профиля, а вклад их в величину отрицательного демпфирования увеличивается. Величина \maz\ sin Ф при этом уменьшается до нуля и становится отрицательной. По мере приближения скачков уплотнения к задней кромке и их стабилизации (см. рис. 2) интенсивность скачков уменьшается, что приводит сначала к замедлению роста отрицательного демпфирования, а затем к его резкому уменьшению. При числе Мсо = 0,92 для данного профиля наблю-

дается появление положительного демпфирования профиля. Таким образом, область неустойчивости колебаний профиля ио числам Моо составляет приблизительно 0,04. Результаты расчета для числа Моо <0,88 согласуются с результатами [4]. В настоящей работе расчет проведен В более широком диапазоне изменения чисел Моо» что позволяет оценить область потери устойчивости профиля и элерона, а также показать, что „выход“ из области неустойчивости обусловлен стабилизацией скачков уплотнения в районе задней кромки профиля. Сравнение результатов расчета по нелинейной и линейной (пунктирная кривая), теориям, приведенное на рис. 5, показывает, что при появлении скачков уплотнения линейная теория становится непригодной. Расчет по нелинейной теории в случае, когда ось вращения расположена между серединой и передней кромкой профиля, например, на расстоянии 1/4 хорды от носка, показывает, что в этом случае отрицательное демпфирование Профиля появляется При меньшем числе Моо-Характер изменения величины \тьш\ sin Ф при изменении числа М<х обусловлен .перемещением скачков уплотнения на поверхности элерона и аналогичен изменению коэффициента \ml\ sin® (см. рис. 6). Это показывает, что причина возникновения отрицательного демпфирования профиля и элерона одна и та же. Полученная численным расчетом зона неустойчивости колебания элерона соответствует зоне „В“ в классификации автоколебаний элерона в трансзвуковом потоке, данной в [!]• Таким образом, этот вид неустойчивости описывается в рамках теории идеального газа. Как уже отмечалось, увеличение амплитуды колебания а0 или приводит к смещению области колебания скачков уплотнения к задней кромке профиля. Это, в свою очередь, приводит к появления отрицательного демпфирования профиля или элерона для меньших чисел Моо. Однако при увеличении амплитуды а0 и 80 вели чина максимального отрицательного демпфирования уменьшается Кроме приведенных нестационарных аэродинамических коэффициентов т“ и т*ш, рассчитаны остальные шесть коэффициентов

т<г, тш> тш’ т\' <> также необходимых при решении задач!

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

Автор благодарит Ю. Б. Лифшица за руководство работой

ЛИТЕРАТУРА

1. Lambourne N. С. Control — surface buzz, Aeronautical Research Council Reports and Memoranda N 3364, 1962.

2. Кузьмина С. И. Численное исследование колебаний профиля в трансзвуковом потоке газа Труды ЦАГИ, вып. 1862, 1977.

3. Lin С. С., Reissner Е., TsienH. S. „On two-dimensional nonsteady motion of a slender body in a compressible fluid". J. „Mathematics and Physics*, vol. 27, N 3, 1948.

4. Ballhaus W. F. and G о о r j i a n P. M. „Implicit finite difference computations of unsteady transonic flows about airfoils, including the treatment of irregular shock wave motions*, AIAA Paper 77—205,

3977.

5. Tijdeman H. On the motion of shock waves on an airfoils with oscillating flap, Symposium Transsonlcum II, 1975.

Рукопись поступила 14j VIII 1980 г