Исследование нелинейных аэродинамических характеристик сверхзвуковых летательных аппаратов мeтoдoм локальной линеаризации Текст научной статьи по специальности «Физика»

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

Том XXV 199 4 №3-4

УДК 629. 735. 33. 015. 3

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

С. С. Граськин

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

1. Расширение диапазона эксплуатационных режимов полета сверхзвуковых летательных аппаратов (ЛА) потребовало получения аэродинамических характеристик на умеренных и больших углах атаки. Для решения этой задачи наряду с экспериментальными методами исследования используются численные методы, среди которых видное место занимают методы расчета, основанные на интегральном представлении решения нелинейного уравнения для потенциала возмущенной скорости [1 — 3]. Однако практическое использование этих методов связано с необходимостью применения итерационных процедур как для построения вихревых структур, так и для решения самого уравнения (см., например, [1, 3]). Все это приводит к большим вычислительным затратам и исключает возможность их реализации на ПЭВМ и в системах автоматизированного проектирования, а также при проведении различных параметрических исследований. Поэтому идея разработки быстродействующих методов и алгоритмов, несмотря на рост производительности вычислительной техники, по-прежнему сохраняет свою актуальность. С этой точки зрения в работе [4] для тонких крыльев и дозвукового потока был предложен способ замены в каждой расчетной точке общего нелинейного уравнения его линейным приближением, для решения которого затем применен метод вихревой решетки или метод дискретных вихрей [5].

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

2. Постановка задачи. Основные положения метода расчета. Пусть ЛА, состоящий из системы тонких несущих поверхностей (рис. 1), движется в идеальном газе с постоянной дозвуковой или сверхзвуковой скоростью . Будем считать, что в процессе этого движения массовые

силы и теплообмен отсутствуют, завихренность потока П сосредоточена только на несущих поверхностях 5 и вихревых пеленах ст, сходящих с фиксированных мест отрыва X, и на сверхзвуковых скоростях; возникшие на кромках элементов Я скачки уплотнения и характеристики эту завихренность не изменяют. Тогда в произвольной точке М пространства вне 5 и ст существует потенциал ср (А/) возмущенной скорости У'{М), который в скоростной подвижной системе координат 0ХС¥С2С удовлетворяет нелинейному волновому уравнению [6]:

где Мм — число Маха невозмущенного потока, у — отношение удельных теплоемкостей. Индексы х, у, г означают первые, а хх, уу, тх, ху, XI, уг — вторые производные потенциала возмущений ср(А/). В соответствии с [7] интегральное представление решения уравнения (1) включает в себя поверхностные интегралы по несущим элементам £ и вихревым пеленам с и объемный интеграл по некоторой области

Рис. 1. Схема обтекания ЛА

(1)

+^2~ <?1 +2м^[ф^,фу(1 + фх) + сругфуфг+ +ФхгФг(1 + Фх)],

У +1 -Л

пространства Т. В численном методе эти интегралы определены через размещенные на Б и ст вихревые рамки или П-образные вихри и объемные источники, интенсивности которых необходимо определять в процессе решения задачи на основе итераций. Однако, согласно работе [4], нелинейное уравнение (1) можно упростить

(1 - + <руу + фк = М^фхе(у + 1)фх +

^ “ !)фх + 2М» (ф^фу + ФхгФг)

и после алгебраических преобразований записать в виде

I1 - Мср)ф** + Ф*у + Фгг = °> (2)

Ы? =

х 2 ФхФхк + фуфду + ФгФхг Фхк

ср 1-М2(у-1)фд

(3)

Таким образом, уравнение (2) с учетом (3) в локальном смысле становится линейным и для его решения теперь применим численный метод, в котором основные положения дискретных вихрей распространены на сверхзвуковые скорости потока [8]. Следовательно, при таком подходе для получения требуемого решения нет необходимости в использовании объемных источников, а вихревые особенности размещаются только на несущих поверхностях Я и вихревых пеленах ст. Граничные условия при этом являются традиционными для МДО [5].

Алгоритм расчета следующий. В соответствии с принятой схемой разбиения каждый несущий элемент ЛА (крыло, оперение и т. п.), геометрия которых задана в связанных осях ОХК2Т, разбивается на ряд элементарных вихревых ячеек с постоянной циркуляцией Гу. Для определения Гу используем условие непротекания на несущих элементах, выполнять которое будем в совокупности контрольных точек, размещенных в центре ячеек. В результате получим

[£]{Г} = {Г},

где [£] — матрица аэродинамического влияния (скосов), определяемая, согласно [8], при Мм = Мср; {г} — вектор-столбец искомых циркуляций Гу; {V} — вектор-столбец, каждый элемент которого V} равен нормальной скорости, создаваемой в /-й контрольной точке набегающим потоком. Начальное приближение циркуляций на пелене

Г^0) =0, а величина =М00. В процессе формирования матрицы [X] следует учесть, что в основе построения дискретных схем для рассматриваемых ЛА лежит прямолинейный вихревой отрезок, который В сверхзвуковом потоке В произвольной точке А/о(хо, Уо> %о) индуцирует возмущенную скорость [8], равную:

(Ф1 - ф2) +

Ь = (*1 - Х0)(х2 - *!) + (1 - М2)[(л - У0)(У2 ~ Ух) + +(^1 - ^0 )(^2 “«!)],

с = (Х1- х0)2 +(1 - м2п)[(л - у0)2 + & - го)2}

(5)

О , л/а + Ь + с = 0, 0 , у[с = 0.

Здесь А*, В*, С* определены согласно [5]; ху> *1,2, Ух ^, Zlt2 — угол стреловидности и координаты концов вихревого отрезка; А Xj — средняя хорда ячейки, которой принадлежит рассматриваемый вихревой отрезок. Второе слагаемое в формуле (4) есть результат выделения особенностей интегрального решения для уравнения (2) на скачке уплотнения [9, 10]. Это слагаемое имеет смысл лишь тогда, когда подкорен-

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

Помимо этого, применение формулы (4) имеет особенности, связанные с взаимным положением вихревого отрезка и характеристического конуса с вершиной в контрольной точке. Вихревой отрезок может находиться внутри, вне или пересекаться этим конусом. Если отрезок находится внутри конуса, то при расчете возмущенной скорости V' имеется особенность лишь тогда, когда контрольная точка лежит на самом вихревом отрезке или на его продолжении. Этот случай соответствует обращению в нуль выражения ? = 4ас - А2. Если отрезок не попадает в конус, то возмущения от него не достигают расчетной точки и скорость V' принимается равной нулю. При пересечении отрезка характеристическим конусом имеют место особенности, связанные с обращением в нуль выражений л]а + Ь + с и 4с . В этом случае справедливы соотношения (5).

положительно, т. е. кромки ячеики

несущих элементах £ и на к-й итерации;

В процессе расчета итерации содержат в себе построение вихревых структур и одновременно уточнение Мср в соответствии с выражением (3). Построение вихревых структур проходит по формуле

{41)}=м_1(м -[4*‘1)]{г«*'1)}} <«>

где |г^ | — вектор-столбец искомых циркуляций вихревых ячеек на

Г^_1)}, [4*_1)] - вектор-

столбец циркуляций вихревых ячеек на пеленах ст и матрица их аэродинамического влияния на (к — 1)-й итерации. При этом произвольной

„ « ~ ~(к) ц-и рамке этой поверхности с циркуляциеи присваивается значение циркуляции у-и вихревои рамки, получаемой на

предыдущем расчетном шаге и примыкающей к ст со стороны поверхности Я. Граничные условия на ст удовлетворяются тем, что вихревые рамки перемещаются в пространстве по траекториям жидких частиц. Из формулы (6) следует, что на первой итерации вихревые пелены ст отсутствуют, а на второй — пускаются по вектору скорости набегающего потока Ут или лежат в плоскости несущих элементов £. На

последующих итерациях по известным циркуляциям |гу^| и |г^_1^|

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

влияния [х£] и |4*_1)}. Следует отметить, что здесь в отличие от традиционного МДО требуется расчет элементов матрицы [Ху] на каждой

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

V — V V' л — V'

гхц+1 уц+1 уц

=-------гг---------> Уху

Ах тл' Ах

Фхг =

^гц+1

Ах

где Ах, Ух, Уу, У£ — расстояние и составляющие возмущенной скорости У' в ц-й и (ц + 1)-й контрольных точках.

После проверки итерационного процесса на сходимость определяем аэродинамическую нагрузку

Ар = ——— = р2

Я +

где р+, р_ — возмущенное давление на верхней и нижней поверхностях несущих элементов .У, — скоростной напор невозмущенного по-

тока. Для расчета Ьр воспользуемся интегралом Коши—Лагранжа, который для стационарного потока имеет вид [5]:

*+,-

V2

г

Y-1

Можно показать, что

-у1 = ±—[гх (cosa cosР + V¿) + ^(sina cosр + ) +

rz(-smp + V¿)l Тогда, с учетом (7), получим:

1 т м» ^-^{rx(cosacosp + V¿) +

(7)

_ 2

уК

. т

(sin a cos р + ) + Гг (-sin 0 + V¿ )}|т -1.

(8)

Переход к суммарным характеристикам произведем путем интегрирования распределенной нагрузки по всем несущим элементам.

3. Примеры расчетов. По изложенному выше алгоритму был проведен ряд расчетов аэродинамических характеристик и вихревых структур, формируемых тонкими несущими поверхностями и сложными пространственными конфигурациями. Некоторые результаты этих расчетов показаны на рис. 2—4. Так, например, на рис. 2 приводятся зависимости коэффициентов подъемной силы суа от угла атаки а, момента тангажа т^, и числа Мю полета для трапециевидного крыла с канардой. Результаты расчетов (черные кружки и треугольники) сравниваются с данными эксперимента (сплошная и штриховая линии) [11]. Здесь же светлыми кружками и треугольниками показаны результаты, полученные по методу [8]. Видно, что на малых и умеренных углах атаки оба метода практически совпадают. Однако на больших углах атаки метод локальной линеаризации по сравнению с методом [8] на больших дозвуковых скоростях увеличивает суа примерно на 10%, а на сверхзвуковых, наоборот, на такую же величину уменьшает, приближаясь к экспериментальным данным. Это связано прежде всего с тем, что по мере роста а возрастает по сравнению с числом Мм невозмущенного потока величина Мср. Следовательно, при больших а несущие поверхности (крыло и канарда) исследуются при большем значении числа М потока. Это, в свою очередь, означает, что при М<1 значение коэффициента подъемной силы суа возрастает, а при

го момента для трапециевидного крыла с канардой М^-1,75

■/Й—8° . 10° 12° ос.

о» ------------

• •

О

О

т2

6)

♦ ♦ \ \ ^' нЩ1“«"

^\\\! | \ * 1

Сечение х=3,52 М^-0,9; а.“10а: р=Ч°

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

М > 1 — уменьшается. Момент тангажа при этом изменяется незначительно. Это говорит о том, что распределение аэродинамической нагрузки качественно по МДО и методу локальной линеаризации совпадает. В процессе моделирования непрерывный вихревой слой на крыле заменялся 72, а на канарце — 12 вихревыми рамками. При этом пелены сходили с задних и боковых кромок крыла и канарды.

Аналогичные результаты, но для аэродинамической компоновки, приведенной на рис. 1, показаны на рис. 3. Здесь рассматриваются зависимости коэффициента нормальной силы сп и момента тангажа т1 от угла атаки а для числа М«, = 1,75 (рис. 3, а) и Мда =3,0 (рис. 3, б). Полученные коэффициенты отнесены к характерной площади «Ухр = 0,002 м2 и линейному размеру / = 1 м. Данные расчетов по предлагаемому методу (черные кружки) сравниваются с экспериментом (сплошная линия) и с результатами, полученными по методу [8]. Видно, что существенное различие расчетных данных зависит не только от угла атаки а, но и числа Мте набегающего потока. Так, например, чем больше число М«,, тем больше при фиксированном а отличие результатов расчета друг от друга и тем ярче проявляются положительные качества метода локальной линеаризации.

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

значении угла скольжения р = 4° и числе Мю = 0,9 для аэродинамической компоновки, схема которой показана на рис. 4. Отличительной чертой этой компоновки является наличие на крыле наплыва большой стреловидности, который на умеренных и больших углах атаки генерирует систему вихрей, существенно влияющих на характер обтекания несущих и управляющих поверхностей, что приводит к неблагоприятному изменению характеристик устойчивости и управляемости в боковом движении. В настоящее время это обстоятельство ограничивает диапазон эксплуатационных углов атаки. Из анализа приведенных результатов расчета следует, что при увеличении угла атаки в диапазоне от 0 ...15° коэффициент момента крена тх изменяется от положительных (крен на правое крыло) до отрицательных значений (крен на левое крыло). Этот факт объясняется перераспределением аэродинамической нагрузки между горизонтальным оперением, крылом, наплывом, фюзеляжем и вертикальным оперением, в результате чего при небольшом

скольжении на правое крыло (р = 4е) при а > 6° создается момент крена, стремящийся поднять это крыло вверх. При дальнейшем росте а вклад вертикального оперения в создание момента тх уменьшается, но вместе с тем возрастает доля момента, создаваемая наплывом и крылом. В результате при а = 18° знак момента тх снова изменяется на противоположный. Аналогичным образом объясняется и характер изменения коэффициента момента рыскания ту. Однако здесь следует

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

Таким образом, исходя из анализа приведенных материалов, можно сделать вывод о целесообразности применения метода локальной линеаризации при исследовании JIA на дозвуковых скоростях при а > 16 4-18°, а на сверхзвуковых — при а > 10 +12°.

ЛИТЕРАТУРА

1. Белоцерковский С. М., Коржнев В. Н., Шипи-лов С. Д. Метод расчета отрывного обтекания крыльев дозвуковым потоком газа // Изв. АН СССР, Механика жидкости и газа.— 1984, N° 4.

2. X е м ш М., Нилсен Дж. Аэродинамика ракет. Кн. 2. Методы аэродинамического расчета. — М.: Мир, 1989.

3. Г р а с ь к и н С. С. Численное исследование стационарного отрывного обтекания пространственных несущих систем сверхзвуковым потоком //

Изв. АН РФ, МЖГ,- 1993, № 3.

4. Ye Zheng-Yin, Zhao Ling-Cheng. Localized linearization method for wings at high angle of attack // AIAA J.— 1990. Vol. 28,

N 10.

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

1978.

6. С и р с У. Р. Общая теория аэродинамики больших скоростей/ Пер. с англ. — М.: Воениздат, 1962.

7. А д а м а р Ж. Задача Коши для линейных уравнений с частными производными гиперболического типа / Пер. с франц. — М.: Наука, 1978.

8. Г р а с ь к и н С. С. Расчет отрывного обтекания треугольных крыльев сжимаемым потоком // Ученые записки ЦАГИ.— 1991. Т. 22, № 4.

9. Алексеев Г. Ю., Еремин В. Ю. Дискретный вихрь в сверхзвуковом потоке // Труды ЦАГИ.— 1987. Вып. 2321.

10. Г р а с ь к и н С. С. Математическое обоснование метода дискретных вихрей решения линейных задач сверхзвуковой аэродинамики // Ученые записки ЦАГИ.- 1992. Т. 23, № 4.

11. Е г i с s о n G. Е., S с h г е i п е г J. A., Rogers Z. W. Canard-wing vortex interactions at subsonic through supersonic speeds // AIAA Paper.—

1990, N 2814.

Рукопись поступила 7/V1992 г.