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

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

Том XXXI 2 000

Ко 1—2

УДК 629.735.33.015.4:533.6.013.425

ОПРЕДЕЛЕНИЕ АЭРОДИНАМИЧЕСКИХ СИЛ, ДЕЙСТВУЮЩИХ В СВЕРХЗВУКОВОМ ПОТОКЕ НА КОЛЕБЛЮЩИЕСЯ НЕСУЩИЕ ПОВЕРХНОСТИ, РАСПОЛОЖЕННЫЕ В РАЗНЫХ ПЛОСКОСТЯХ

В. П. Кузьмин, С. И. Кузьмина, В. А. Мосунов

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

Разработана гибкая вычислительная процедура для обеспечения высокой точности результатов при различных сочетаниях числа Маха М, числа Струха-ля к и размеров панели. Разработанная программа включена в комплекс КС-2 расчета на флаттер, причем организация исходных данных одинакова как при М > 1, так и при М < 1.

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

ОБОЗНАЧЕНИЯ

координаты точек на несущей поверхности (панели); координаты контрольной точки, в которой вычисляется скос;

координаты панели относительно контрольной точки;

вектор единичной нормали к панели;

вектор единичной нормали в контрольной точке;

х,у,г — ЛЬ г0 ~

^ = х0-х

п = уо->'

С = 20-г

П=\Пу,П2

«0 =|«у0’иг0

гиперболический радиус;

число Маха; число Струхаля;

частота колебаний несущей поверхности; скорость набегающего потока;

характерный линейный размер, к которому отнесены все линейные размеры;

скоростной напор;

безразмерная разность давлений на несущей поверхности;

коэффициенты обобщенных сил; размах панели.

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

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

В случае полета со скоростью, меньшей чем скорость звука, широко применяется метод дипольной решетки [1]. В ЦАГИ уже в течение 20 лет в задачах динамической и статической аэроупругости с успехом используется оригинальная версия этого метода и созданная на ее основе программа расчета [2]. При сверхзвуковых скоростях аналогичного надежного, универсального метода до сих пор не было, хотя публикаций, касающихся решения этой проблемы, достаточно много [3] — [7]. Описанные линейные методы можно условно разделить на три группы, основанные на использовании: 1 — потенциала скорости; 2 — градиента потенциала скорости; 3 — потенциала ускорения. Многие из алгоритмов, применяемых в случае компланарных несущих поверхностей и требующих учета вихревого следа за задней кромкой, становятся неприемлемыми при определении аэродинамических нагрузок на поверхностях, расположенных в разных плоскостях [7], [8].

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

г2=п2^2. я = \^2-3 2г2 З2 =м2 -1,

м

со/

у!

к =

к' ■■

Ш

ч=-

2

А с,

Отп

1п

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

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

1. Основные соотношения. Для вычисления элементов матрицы скосов от панели несущей поверхности используется интегральное соотношение для потенциала скорости в произвольной точке от единичной безразмерной разности давлений [3]

п сД[/(л)

Ф(дго,7о^о) = -Т | (1)

2 * ■> дп

П£ ^(П)

где

\

Р = ехр(-Ш;) | ехр (Зг

' щ;' "со^'Л')’

Ь2; Л'

<3^ .

В формуле (1) и последующих предполагается, что оси %,г| лежат в плоскости панели, а начало координат находится в контрольной точке.

Нормальная составляющая возмущенной скорости в контрольной точке (скос) определяется производной

м;^дФ(хо,У(),*о) (2)

дп0

Для вычисления производной вдоль нормали используется соотношение

дР 1 дР. ...

(3)

дп г дг 7 Подставив соотношение (3) в (1), получим

1 цЧ^и г(т|) 1 ир 1

Щх0,у0,г0) = (пуУ] + п2О- ] ] ^(^^(ПуЦ + ПгО^ . (4)

Л/, 4/.(П)Г Г

Интегрирование в (4) ведется по поверхности одной панели, где сомножитель {пуу\ + п2С,) сохраняет постоянное значение, поскольку различные точки на панели отличаются на величину вектора, перпендикулярного нормали п.

В работе [3] двойной интеграл в (4) приведен к виду

1 ~ Л ---2~“П >

где

ехр(-/^\Й;) зт(£7?)

Я

Ъи(Ю

5*,(ч) Ыл)

+ /М | ехр(-г^'М^) 8т(ЛТ?)<зй;

4Л(П) ?£(Л)

+

+ 2 | ехр(-г^)

Чь(Ю

+ •< 1-ехр

б'! ехр

/Щ-МЯ)

-п+

т=\

ехр

&(§ + МЛ) > +

1_ Р2 ] 1 )

тЬ(^ + МЯ) /*(§ + МД)

р 2г Р2

(6)

-1

О,

1-ехр

+ 8-

р 2г

(7)

бдв^-МЛ), />! =2 тЬ + Исг, и02=ЫтЬ + И<г.

Коэффициенты ряда в (7) взяты из работы [3]:

Ь = 0,0350039075, ^=0,00432952, а2 =0,00160137, а3 =0.033195, а4 =0,0986923, а5 =0,37673986, а6 = 0,822464, а7 =-0,3802627, а8 = 0,04340039.

Производная (2) вычисляется аналитически [4]:

П£/

1

*Хл)

г2

<Й1 +

+ | (пу1) + п2О(ну0Т\ + пг0О

^Р-г2-2гР{ц)

дг

-сЙ1

(8)

Я[Г(г\\

Производная —-— в (8) вычисляется численно. дг

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

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

И/Аср = 2ъ2, (9)

где IV — комплексная матрица скосов, Аср — столбец комплексных амплитуд безразмерных разностей давлений на панелях, Z — столбец комплексных амплитуд скоростей в контрольных точках, нормальных к несущим поверхностям.

Если задана форма колебаний поверхности ^(£,т},0 = /(^, т|) ехр(/А:/),

„ д/ то компоненты вектора 2 определяются соотношением г : =—{с, ,■ ,г| .• ) +

J

+ ikf(£)j,^r]j), где7 — номер контрольной точки.

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

0,тп

2. Численные процедуры. Ниже описаны наиболее существенные элементы численной процедуры, от точности которой зависит точность вычисления аэродинамических сил.

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

Маха лежит точно в плоскости панели или расстояние от нее до панели значительно меньше размаха панели.

2.2. Определение области интегрирования. При вычислении скоса интегрирование в (6) и (8) производится по части панели, попадающей в обратный конус Маха с вершиной в контрольной точке (рис. 1).

Рис.

При численном интегрировании для каждой панели вводится местная переменная 8, которая меняется в диапазоне [0, 1], и интегрирование по плоскости панели ведется в переменных (4, е), при этом

г = {[тц+8(т12-г11)]2+[^1+8(С2-С1)]2}1/2 и <Ь\ =1р&,

где 1р = д/(г|2 -Л1)2 +(^2 —размах панели.

Таким образом, основной интеграл (5) записывается в виде

Л[/

1

м> = — 2

Г(е)

.2 Р

lr.dE +

П(/

+ | (пуг\ + п2О(пу0Ч + пг0О-

Лх

дПг)г2 _2г7г(8) дг

-1рс1г

(И)

Часть панели, попадающая в обратный конус Маха, определяется очевидным неравенством £>Рг (рис. 1). В данном методе рассматриваются только панели с боковыми кромками, параллельными оси Для таких панелей попадание любой их части внутрь конуса возможно лишь в том случае, если туда попадает часть передней кромки. Следовательно, для определения условий пересечения панели конусом достаточно определить две точки пересечения отрезка передней кромки и конуса, т. е. величины е/, и

Zu . Очевидно, что если один или оба конца передней кромки лежат внутри конуса, ТО следует ПОЛОЖИТЬ Е^-Ои (или) Ъц - 1.

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

Л1(Л2 -Л])+?1(С2 _?i)

sc =-----—------—■ - ---—, определяющая точку минимума величины г.

I

1Р

2.3. Численное определение интегралов с помощью квадратурных формул. Для численного вычисления интегралов в данной работе используются квадратурные формулы Гаусса двух типов [12]: 1 тип — стандартная формула Гаусса, дающая точные значения интегралов для функций

2 3

\,v, v , v ,... на отрезке [0,1] и 2 тип — квадратурная формула, которая может быть получена из стандартной заменой переменной v - sin S'. В этом случае новые значения узлов и весов квадратурной формулы связаны со стандартными очевидными соотношениями

~ л

vn =sin(7TU„ /2), ап = an—cos(nvn/2).

Использование квадратурной формулы означает замену интеграла конечной суммой

Sи м

f/(^*(^-^)£*M/(U).

\L т=Х

где = +(Ли

Узлы и веса для обеих квадратурных формул приведены в табл. 1 для семи узлов.

Использование квадратурных формул двух типов связано с видом подынтегральных функций в интегралах (6) и (11). Так, если в интегралах (6) пределы интегрирования определяются границами панелей, то используется квадратурная формула первого типа, а если один из пределов интегрирования находится на конусе Маха, то — второго.

Таблица 1

Узлы и веса для квадратурных формул

п vn ап vn ап

] 0,974554 0,06474241 0,999201 0,00406381

2 0,870766 0,139853 0,979466 0,0442896

3 0,702923 0,190915 0,893081 0,134918

4 0,500000 0,208980 0,707107 0,232118

5 0,297077 0,190915 0,449895 0,267825

6 0,129234 0,139853 0,201610 0,215169

7 0,0254460 0,06474241 0,0399599 0,101616

При использовании квадратурных формул интервал интегрирования по переменным £ или £ приводится к отрезку [0,1 ] для переменной V, причем, при использовании квадратурной формулы второго типа значение и = 1 должно соответствовать пределу интегрирования, лежащему на конусе Маха.

2.4. Вычисление функции /'’(е). Функция Р(е) вычисляется с помощью квадратурных формул Гаусса, описанных выше. Заданное значение переменной 8 однозначно определяет величину г и пределы интегрирования по переменной 2,

£1/(е) = £1+(£з-^)е,

^ £ (£) = шах{эг, %£ },

где '^1 = %2 +(^4 - £,2 )8 координата задней кромки, соответствующая данному значению в.

Если рг>£>1, то нижний предел интегрирования в (6) лежит на конусе Маха, подынтегральная функция в окрестности нижнего предела пропорциональна в этом случае используется квадратурная формула

второго типа. Если (Зг < ^, то нижний предел интегрирования в (6) лежит на границе панели, и в этом случае используется квадратурная формула первого типа.

Трудоемкость численного интегрирования (6) связана, в основном, с

л

наличием в подынтегральном выражении сомножителя ехр(-г£'/р £,), ко-

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

для каждой из которых интегралы (6) вычисляются по квадратурным формулам. Количество делений определяется по эмпирической формуле

О,7—у——(4{у -£і) + 1 М2-1

, где іпІ(.) — целая часть числа. Очевид-

но, что необходимость в дополнительном делении интервала интегрирования на части возникает лишь при больших числах Струхаля и близких к единице числах Маха.

Отметим также, что если нижний предел интегрирования лежит на конусе Маха и > 1, то квадратурная формула второго типа используется

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

г п * - 5^(8) гт 5^(е)

2.5. Вычисление производной---------. Производная -------- определяет-

дг дг

ся численным дифференцированием по переменной г с помощью центральной разности. Стандартная величина приращения по переменной г определяется по формуле

* л-д/^- Р2'2 ГПЧ

Аг = Аг-------------, (13)

где AF = 10-5.

При вычислении функции F(e) по формуле (6) нижний предел интегрирования может лежать либо на задней кромке панели, либо на конусе

5F(£)

Маха, и в точке перехода производная -------терпит разрыв второго рода.

дг

Поэтому при окончательном выборе величины приращения Аг контролируются два условия: 1) — R(r + Аг) > 0, R(r - Аг) > 0 для передней и задней кромок и 2) — при изменении величины г не должен меняться тип нижнего предела интегрирования. Если одно из этих условий не выполняется при стандартном значении приращения (13), приращение А г уменьшается до приемлемой величины. Описанные условия выполняются, если величины г + Аг и г - Аг лежат в пределах от rmax до rmin, которые определяются формулами

О при Cz. ^Р>% при %L<$r.

r£z/p при %L>$г,

Ду/р при %L<$r.

(14)

3^(є) л

При г = 0 принимается-- = 0.

дг

2.6. Интегрирование по 'размаху панели. При вычислении интеграла (11) существенное значение имеет расстояние с/ от оси конуса Маха до плоскости панели

іі = пуг\ + п(15)

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

Вся панель делится на области по всем точкам пересечения передней и задней кромок панели конусом Маха. Интегрирование внутри каждой об-

ласти осуществляется либо по квадратурным формулам, либо аналитически с использованием квадратичной аппроксимации функций Р(е) и ее произ-5^(8)

водной

дг

Количество областей, на которые может быть поделена панель, зависит от количества точек пересечения передней и задней кромок панели и конуса Маха и может изменяться от 1 до 4. На рис. 1 приведен пример с максимальным возможным числом делений. Каждая область интегрирования характеризуется граничными значениями переменной е. Область, внутри которой расположена точка ес, будем называть центральной. Во всех областях, кроме центральной, интегрирование ведется численно по квадратурным формулам одного из двух приведенных выше типов. Выбор типа квадратурной формулы определяется поведением функции Р(е) в окрестности точек пересечения панели и конуса Маха. На рис. 2 приведены характерные зависимости функции Р(г) для нескольких случаев пересечения панели и конуса Маха.

На рис. 2 изображены зависимости действительной (сплошные линии) и мнимой (штриховые линии) частей функции /г(8) от переменной е. На правой части рисунков показаны: расположение панели относительно контрольной точки, линия пересечения конуса и плоскости панели (штриховая линия), значения числа Маха (М), числа Струхаля (к) и высоты контрольной точки над панелью (^).

На приведенных рисунках видны изломы функции F(e) в точках пересечения передней и задней кромок панели конусом Маха. В этих точках 3F(e)

производная -------- имеет разрыв. Описанное выше деление панели на

дг

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

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

Обозначим пределы интегрирования для одной из таких области, как S) и г2 ■ Тогда интегрирование по квадратурным формулам ведется для областей, внутри которых не лежит центральная линия, т. е. если

(ei -ес)(в2 -ес)>0,

или в том случае, если расстояние от оси конуса Маха до панели достаточно велико — > Aemjn, где d -пуц + п2С,, a Aemin - min(je] — sc|,je2 ~EC|)-

Ip

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

Предусмотрено максимальное число делений, равное 9, что обеспечивает возможность точного определения скоса при изменении величины г

■з

более, чем в 10 раз, в пределах одной панели.

Если средняя линия конуса Маха лежит внутри отрезка интегрирования, т. е., если (s] -sc)(s2 -ес)<0 и на расстоянии, удовлетворяющем неравенству

то интеграл (11) также вычисляется с помощью квадратурных формул, но независимо для двух областей [е^,ес] и [вс,Б2], причем это деление на две области производится до проверки условий, определяющих необходимость дополнительного деления отрезка интегрирования из-за значительного изменения величины г.

2.8. Аналитическое интегрирование с помощью квадратичной аппроксимации. Если контрольная точка лежит в плоскости панели или очень близко к ней и координата ес центральной линии конуса находится внутри отрезка интегрирования, то для вычисления интеграла (11) используется

л гг \ - 9^(е)

квадратичная аппроксимация функции г (е) и ее производной ----------------,

дг

а также аналитические формулы.

Если (е] -ес )(в2 -£с )< 0 и — Лет;п > — , то отрезок интегрирова-

30 Iр

ния [б|,82] делится натри части (см. рис. 1) [е],ес -Де], [ес -Ле,ес +Де] И [ес + Де,82 ], где Д8 * ^-Д8т;п .

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

что внутри этого отрезка производная —■ — ^ не имеет особенностей, что

дг

важно для аппроксимации ее полиномом.

г, \ ^(е)

Функция г (г) и ее производная---------представляются в виде

дг

Р(е) * ?(е) = а0+а1{£-Ес)+а2(е~ес)2 (16)

дР(е) дР(е) І(є —єс)+62(е — ес)2при £<ес>

дг дг , А. (е-О иА-їс-о Ї2

г>0+Аі(є-єс)+г»2(є-єс) при є>єс.

(17)

Для определения коэффициентов квадратичной аппроксимации используются центральная точка єс и точки, отстоящие от центральной, на расстоянии « 0,10,2 от длины отрезка интегрирования. Для повышения

. Э^( 8)

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

дг

стояния й используются две формулы (17).

Подставив соотношения (16) и (17) в (11), получим аналитическое выражение для скоса

а2

(Л с12+12

+

С2а]

М

+ Сл

Ьп1

Я2

с/2+12 (?|~»|)

г +

+/■

-аХ ап —

А

Г 12 (I2 л-12

+

-И

+ Ь\ )+ С\с12 ('р2 ^ ^2 ) 2

Л

2 +12

+ 1ов

/ +

+/■

где Сх = (и-и0), С2

1(^2 “Л Ко+(22-21)«го)

(18)

и

, / = Д 8, й? =----------------.

1Р

-4'

При значениях величины с/ , близких к нулю {<111 <10 ), соотноше

ние (18) заменяется своим предельным значением для й = 0:

С,

Н' = -

/

21 п

(19)

После вычисления приближенного значения скоса (>у) по формулам

(18) или (19), соответствующего квадратичной аппроксимации функций ар(е)

г (е) и -------, производится дополнительное вычисление интеграла (11)

дг

(Ам>) по квадратурной формуле для разностей /’(е) - /г(е) и

ЭДе) дР(е) _ ,

---------------. В этом случае используются квадратурная формула перво-

дг дг

го типа с четным числом узлов (N=6), чтобы исключить вычисления в средней точке при с1 = 0. Окончательное значение скоса определяется суммой

м? = ц> + Атф' .

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

На рис. 3, 4 приведены примеры зависимости скоса от высоты контрольной точки над панелью (С,) для различных панелей, различных значе-

Рис. 3

й

И=1,1В А = 2,00 Л о' < •

50 -0,25 (Л ,00 0,25 0

Т~

Н>

. о ф- 90 Л Г\ 1

«о зг?-- ч V Н ,0 /Т. 5 1 V

* ф-0

,0 'Ы Г 1

\ 1

\ /

К

Рис. 4

ний числа Маха, числа Струхаля и различных углов между нормалями к панели и в точке скоса (ф = агссоз(и ■ «о) )• Сплошные линии на рис. 3, 4 соответствуют действительной части скоса, а штриховые — мнимой.

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

Рассчитывались две стандартные конфигурации АОАШЗ, представляющие интерферирующие между собой поверхности, расположенные в разных плоскостях, параллельных вектору скорости набегающего потока. Исходные данные для несущих поверхностей и результаты расчетов обобщенных аэродинамических сил, используемые для сравнения, взяты из работы [4].

Вариант X У z

^г.о = ^ а 0,00 0 0

Ь 2,25 0 0

с 2,75 0 0

d 3,70 0 0

е 2,70 0 0

/ 4,00 0 0

g 3,90 1,00 0

h 4,25 1,00 0

Кг. о=30° g 3,90 0,866 0,50

h 4,25 0,866 0,50

Форма крыло стабилизатор

кручение крыла j**-2,25|>|-0,85) 0

изгиб крыла у\у\ д 0

тангаж стабилизатора 0 х-3,35

крен стабилизатора 0 У

.X

Рис. 5

Пример 1. Крыло — стабилизатор ( AGARD wing-tail configuration ). Рассматривались два варианта: а) — крыло и стабилизатор лежат в одной плоскости и б) — угол наклона стабилизатора к плоскости крыла составля-

ет Vr 0 =30°. На рис. 5 показаны вид в плане, геометрические размеры и формулы для четырех форм антисимметричных колебаний.

В расчетах число панелей по размаху было равно 12, число панелей по хорде равно 10 на крыле и 8 на стабилизаторе.

При использовании линейных методов определения аэродинамических нагрузок на сверхзвуковых скоростях нет устоявшегося правила для определения положения контрольных точек, поэтому расчеты по данной программе проводились при двух различных правилах, определяющих положение точки удовлетворения граничных условий^ Контрольная точка лежит на средней хорде панели, и правила касаются только ее положения на ЭТОЙ хорде.

Согласно первому правилу контрольная точка располагается на одинаковом относительном расстоянии от передней кромки каждой панели.

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

— 50% хорды панели, если передняя кромка панели сверхзвуковая (<т = 0,5);

— 85% хорды панели, если передняя кромка панели дозвуковая (ст — 0,85);

— 92,5% хорды панели, если задняя кромка панели совпадает с дозвуковой задней кромкой несущей поверхности (ст = 0,925).

В табл. 2 приведены действительные (Яе(0) и мнимые (1ш(0) части коэффициентов обобщенных аэродинамических сил, полученные по программе, и результаты работ [4] и [9] при М= 1,2Д= 1,5.

Таблица 2

Сравнение коэффициентов обобщенных аэродинамических сил для конфигурации крыло — стабилизатор, М = 1,2, к = 1,5

Плоская конфигурация (УГ 0 =0)

Результаты работы Результаты работы Результаты данной работы

[4] [9] а = 0,95 переменная а

У Яе(0) 1т(0) Яе(0) 1т(0) Яе(0) 1т (0) Яе(б) 1«п(0)

11 -0,13 0,20 -0,09 0,27 -0,10 0,21 -0,10 0,25

21 0,30 0,40 0,40 0,42 0,33 0,39 0,38 0,41

12 -0,13 -0,02 -0,15 0,00 -0,14 -0,00 -0,15 -0,00

22 -0,31 0,24 -0,30 0,28 -0,28 0,25 -0,30 0,27

31 -0,26 0,19 -0,24 0,20 -0,24 0,20 -0,25 0,20

41 -0,41 0,31 -0,42 0,32 -0,39 0,32 -0,40 0,33

32 -0,30 -0,02 -0,31 -0,01 -0,30 -0,01 -0,32 -0,01

42 -0,45 -0,01 -0,46 0,00 -0,45 0,00 -0,47 0,00

33 0,48 0,63 0,50 0,64 0,51 0,63 0,52 0,65

43 0,95 0,57 1,0 0,57 0,99 0,57 1,0 0,57

34 -0,15 0,29 -0,15 0,29 -0,14 0,31 -0,15 0.31

44 -0,097 0,471 -0,098 0,483 0,08 0,48 -0,09 0,49

Пространственная конфигурация (Уг о =30°)

Результаты работы Результаты данной работы

[4] а = 0,95 переменная а

и 1*е(0 1т(0) Яе(® Ьп(б) МС?) 1т(0

11 -0,09 0,30 -0,10 0,21 -0,10 0,32

21 0,30 0,55 0,32 0,39 0,33 0,58

12 -0,12 -0,01 -0,14 -0,00 -0,13 -0,00

22 -0,27 0,36 -0,28 0,25 -0,28 0,38

31 -0,03 -0,08 -0,02 0,06 -0,02 0,07

41 -0,09 0,14 -0,10 0,10 -0,07 0,12

32 -0,07 0,01 -0,07 0,01 -0,06 0,02

42 -0,12 0,02 0,12 0,02 -0,09 0,02

33 0,33 0,69 0,48 0,63 0,37 0,71

43 0,67 0,65 0,95 0,59 0,62 0,57

34 -0,11 0,32 -0,15 0,30 -0,09 0,30

44 -0,08 0,52 -0,10 0,49 -0,05 0,41

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

\

\

\

X У Z

а 3,55' 0 1,20

b 4,85 0 1,20

с 4,75 1,00 1,20

d 5,1 1,00 1,20

е 2,5 0 0

/ 4,2 0 0

g 3,70 0 1,20

h 4,6 0 1,20

Форма киль стабилизатор

изгиб киля z2 0

кручение киля |z|(*-0,875|z|-3) 0

крен стабилизатора 0 У

Рис. 6

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

Пример 2. Стабилизатор — киль (AGARD tail-fm configuration). На рис. 6 показаны вид сбоку и вид сверху комбинации стабилизатор -— киль, их координаты и приведены формулы для трех форм антисимметричных колебаний, а именно: изгиб киля, кручение киля и крен стабилизатора. Каждая поверхность разбивалась на 150 панелей, 15 по размаху и 10 по хорде; положение контрольной точки было переменным и определялось по второму правилу. В табл. 3 приведено сравнение результатов, полученных по представленной программе, с данными других работ при М= 1,6, к=0 и к= 1,5. В данном примере сравниваются модули и аргументы комплексных коэффициентов обобщенных аэродинамических сил, вычисленных для половины конфигурации стабилизатор — киль. Наиболее близкими оказываются результаты данной работы и результаты работы [4]. Следует считать удовлетворительным согласование между представленными результатами, принимая во внимание сложность рассчитываемой пространственной компоновки.

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

стабилизатор — киль, М = 1,6

Результаты к- = 0 к = 1,5

работы V MOD(0 MOD(0 ARG(0, градус

10 0,82 89,4

11 11 0,78 90,9

4 0,80 92,3

данная 0,76 88,6

10 0,09 128,6

11 12 0,09 131,9

4 0,11 129,8

данная 0,12 116,7

10 0,11 -63,3

11 13 0,19 72,8

4 0,32 15,6

данная 0,33 9,7

10 0,79 0,71 13,2

11 21 0,81 0,69 18,2

4 0,78 0,71 18,7

данная 0,76 0,66 14,7

10 0,081 0,26 63,4

11 22 0,097 0,27 67,3.

4 0,098 0,26 63,2

данная 0,109 0,27 56,8

10 0,18 0,05 ■ -167,0

11 23 0,23 0,14 177,1

4 0,43 0,30 -69,4

данная 0,45 0,31 -76,7

10 0,095 61,9

И 34 0,087 64,9

4 0,051 59,2

данная 0,052 49,9

10 0,03 48,2

11 32 0,18 7,9

4 0,02 51,3

данная 0,02 40,0

10 0,72 88,5

11 33 0,67 89,8

4 0,70 88,4

данная 0,66 86,8

ЛИТЕРАТУРА

1.Albano Е., Rod den W. A doublet lattice method for calculating lift distribution on oscillating surfaces in subsonic flows // AIAA Paper. № 69-73.

2. Мосунов В. А., Набиуллин Э. H. Определение аэродинамических сил, действующих в дозвуковом потоке на упругие колеблющиеся поверхности, расположенные в разных плоскостях // Труды ЦАГИ. — 1981. Вып. 2118.

3. Арра К. Constant pressure panel method for supersonic unsteady airload analysis // Journal of Aircraft. — Oct. 1987. Vol. 24.

4. Арра K., Smith M. J. C. Evaluation of the constant pressure panel method for supersonic unsteady airloads prediction // Journal of Aircraft. — Nov. 1988 . Vol. 26.

5. Lо11ati I., Nissim E. Nonplanar, supersonic, three-dimensional, oscillatory, piecewise continuous kernel function method // Journal of Aircraft. — Jan. 1987. Vol. 24

6. L i u D. D., J a m e s t D. К., С h e n P. С., P о t о t z к у A. S. Further studies of harmonic gradient method for supersonic aeroelastic applications // Journal of Aircraft. — Sep. 1990. Vol. 28.

7. Кузьмина С. И. Расчет нестационарных аэродинамических характеристик тонких крыльев конечного удлинения при сверхзвуковых скоростях полета // Труды ЦАГИ. — 1981. Вып. 2118.

8. Б у н ь к о в В. Г. Комбинированный метод расчета аэродинамических сил на колеблющемся летательном аппарате в сверхзвуковом потоке // Ученые Записки ЦАГИ. — 1984. T. XV, JV» 3.

9. Р о 11 о k S. J., Huttsell L. J. Application of three unsteady aerodynamic load prediction methods // Air Force Flight Dynamics Laboratory. — May 1974. TR-73-147.

10. J о n e s W. P., A p p a K. Unsteady supersonic aerodynamic theory using a fmite-element doublet representation // AIAA Journal. — Jan. 1977. Vol. 15.

11. H о u n j e t M. N. L. Improved potentional gradient method to calculate airloads on oscillatory supersonic interfering surfaces // Journal of Aircraft. — May 1982. Vol. 19.

12. Корн Г. и Корн Т. Справочник по математике. — М.: Наука. —

1974.

Рукопись поступила 10/XI1998 г.