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

УЧЕНЫЕ ЗАПИСКИ НАГИ Т о м VII 19 7 6

№ 2

УДК 629.735.33.015.4:533.6.013.425

ОБ УЧЕТЕ ВЛИЯНИЯ УПРУГОСТИ КОНСТРУКЦИИ НА ПРОДОЛЬНОЕ КОРОТКОПЕРИОДИЧЕСКОЕ ДВИЖЕНИЕ САМОЛЕТА

Ю. А. Кублин, В. И. Чубарое

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

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

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

Предлагаемая методика является дальнейшим развитием работ [2, 4].

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

кости симметрии самолета. Ось Ог направлена в сторону правого полукрыла.

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

; Н,;' Г : /(*, г, Ц = У 9;(х, г)%ф), (1)

. ! ', ■ ■ (=1 ■ ■ '

где л:, г — координаты точки срединной поверхности; ?,•(*, г)— заданные линейно независимые координатные функции; щ(1) — обобщенные координаты; п — число обобщенных координат.

Функция ^ (х, г) = 1 соответствует поступательному смещению самолета как твердого тела, а величина ихУ) равна изменению высоты полета. Вторая координатная функция <?о(х, г) =—х соответствует повороту самолета как твердого тела вокруг оси Ог, а величина и2 (?) равна изменению угла тангажа. Функции г) с индексами /> 3 представляют изгибные и крутильные деформации несущих поверхностей, изгибные деформации фюзеляжа и хорд крыла, так что относительное смещение срединной поверхности вычисляется по формуле

/о (х, г, *) = 2?|(-*\ г) щ (*). (2)

1=3

В качестве координатных функций ?;(.х, г) при г!>3 могут использоваться как формы собственных колебаний в пустоте, так и любые функции, имеющие непрерывные вторые производные, удовлетворяющие существенным граничным условиям и выбираемые в соответствии со спецификой распределения жесткостных и массовых характеристик конструкции.

Распределение местных углов атаки при нулевом угле атаки и нулевом угле отклонения рулевой поверхности для стапельной формы самолета будем задавать в виде функции <х0(х, г).

Так как отклонение рулевой поверхности является частным случаем деформации срединной поверхности самолета, то закон распределения смещений зададим в виде функции <?ъ(х, г) 3 (^), где 5(0 — угол отклонения рулевой поверхности. Вне рулевой поверхности <р6 (х, г) = 0.

При расчетах конструкций с крылом малого удлинения „бес-хвостой1* схемы (фиг. 1) с достаточно равномерным распределением массовых и жесткостных параметров в качестве координатных функций с индексами г';>3 удобно применять произведения полиномов

т‘(л'Ъ Ш' (3)

где рь <7; — номера полиномов; — расстояние от носка фюзеляжа до начала координат; £к — полуразмах крыла.

: Полиномы //>дуу для носка и кормы фюзеляжа удовлетво-

ряют балочным граничным условиям на свободном конце балки

<1~Х*

ач,

й X3

^0,

(4)

где л

/ • — г

Функции / (——) ПРИ 2 —— = 0 удовлетворяют кинематичес-

' < \ £к / 1-к

кому условию

V dfЧ.

</ г

-=о,

(5)

а при г= 1 — условиям типа (4) на свободном конце балки.

Фиг. 1

Для Х1=2 = 0 все функции г) с индексами г>3 удовлетворяют условиям

ср.(0, 0) = -^- = -^- = 0. (6)

д х

д г

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

Фиг. 2

ласти располагается фюзеляж, во второй — крыло, в третьей — горизонтальное оперение.

Смещение срединной поверхности в первой области представим функциями вида

где функции аналогичны функциям формуле (3).

Для точек второй и третьей областей смещение, соответствующее функциям (7), определим по формуле линейной экстраполяции.

dv\(xn, z0)

?] (X, z) = fj (х0, Zq) н-------(х - х0), (8)

1 ‘ д х

где х0, z0 — координаты начала оси жесткости крыла или горизонтального оперения (см. фиг. 2).

Считая хорды крыла недеформируемыми, зададим деформации во второй области функциями вида

(9)

где Ьк — средняя хорда крыла; LK — длина оси жесткости крыла;

х2=(х — хк) cos — — Zk) sin Хк5

z2 = (x- xK) sin Хк + (Z — ZK) cos yK;

zK — координаты начала оси жесткости крыла;

Хк — угол стреловидности оси жесткости.

Функция j равна 1 или

Если = 1' то фУнкции

вые вторую и третью производные, а при = 0 обращаются в

Lk

нуль и имеют нулевые первые производные.

Если то функции имеют при = 1

1 \ 1 \

нулевые производные, а при -^- = 0 обращаются в нуль.

Деформации горизонтального оперения представим функциями <pf (х, z), определенными в области 3 (см. фиг. 2) и аналогичными функциям (9) 2).

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

С и = (D - Dk) й + (В - Bk) и + С* Ъ + + В*Ъ + GB В0, (10)

где и, и, и — «-мерный вектор обобщенных координат и его производные; 8 — угол отклонения рулевой поверхности; С—симметричная матрица коэффициентов обобщенных масс; D — матрица коэффициентов аэродинамического демпфирования; /^ — симметричная матрица коэффициентов конструкционного демпфирования; В — матрица коэффициентов „аэродинамической жесткости"; Bk — симметричная матрица коэффициентов обобщенных жесткостей

при —J— = 1 имеют нуле-Lit

конструкции; С1-—вектор обобщенных сил инерции при колебаниях рулевой поверхности; /)5 — вектор обобщенных аэродинамических сил, вызванных угловой скоростью вращения рулевой поверхности; Въ — вектор обобщенных аэродинамических сил, возникающих при отклонении рулевой поверхности на единичный угол; {?в — вектор обобщенных сил веса; В0 — вектор обобщенных аэродинамических сил, определяемых стапельной формой самолета при нулевом угле атаки и нулевом угле отклонения рулевой поверхности.

Элементы матриц й, В и векторов й0, В'', В0 являются функциями числа М, скорости полета V и плотности воздуха р.

Элементы матриц В11г В12, й удовлетворяют соотношениям

В 1\ 0, 1 V— В/ 2-

При вычислении аэродинамических сил заменим крыло, горизонтальное оперение и фюзеляж тонкой несущей поверхностью, имеющей форму проекции самолета на плоскость Охг, и воспользуемся хорошо зарекомендовавшей себя при исследованиях флаттера гипотезой гармоничности [1]. Для дозвуковых скоростей (М<^ <<0,9) применим апробированный на практике алгоритм [1], основанный на замене вихревого слоя системой П-образных вихрей. Для сверхзвуковых скоростей (М< 1,1) воспользуемся алгоритмом [3], основанным на замене вихревого слоя прямоугольными ячейками с постоянной напряженностью вихрей (при этом несущая поверхность аппроксимируется поверхностью с пилообразными кромками).

Для расчета элементов инерционной матрицы СI] распределение масс конструкции и топлива схематизируем системой дискретных грузов. Центр инерции каждого груза расположен на

оси прямолинейного недеформируемого стержня, основание которого жестко связано со срединной поверхностью конструкции (фиг. 3). Одна из главных центральных осей инерции груза направлена по оси стержня. Кинетическая энергия груза Ет с учетом разложения (1) и трех степеней свободы груза (смещение, параллельное оси Оу, и два поворота вокруг главных центральных осей инерции груза) вычисляется по формуле

(И)

£г = — мтД Си,

2

где „т“ — индекс транспортирования матриц; ДС — симметричная матрица обобщенных масс, соответствующая т-иу грузу, элементы которой вычисляются в виде квадратичной формы векторов обобщенных смещений груза

= (12) где 1т — диагональная матрица третьего порядка с элементами

гт = м

*11 /Г1Г>

22

1т = I *

133 Л’

МТ, /5, /л~ масса и моменты инерции тп-го груза; \£г;, векторы

обобщенных смещений для г-й и /-й координатных функций.

Компоненты векторов при А5 = /,у для функций <?к{х, г) имеют вид

= *) + ?*,(*. 2)5> = ?**(*. 2). 8 = ?*„(■*. 2), (13)

где х, г — координаты основания кронштейна; ср^(х, г), ?к„(х, г) — производные функции по направлению оси кронштейна и нормали к ней соответственно; а — длина жесткого кронштейна (см. фиг. 3).

Элементы матрицы С с учетом (12) определяются в виде суммы

Су =2 о4)

/71 = 1

где Л/,.— общее число сосредоточенных грузов расчетной схемы.

Элементы вектора С“ вычисляются по формуле

N . ■

С\= (15)

т=1

где Л^р. п — число грузов на рулевой поверхности; — вектор обобщенных смещений, соответствующих функции срв (л:, г).

В соответствии с выбором функций <р,(.т, г) и ?2(х, г), с1? =

— с21 = 0; си, с22 — масса Мс и главный центральный момент инерции самолета /гс относительно оси, параллельной оси Ог.

Элементы вектора бв пропорциональны элементам С,, матрицы С

Qвl—gClX, (16)

где £ = 9,81 м/с2 — ускорение свободного падения.

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

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

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

П = — ит Вк и, 2

(17)

- _ _ ' ’ где Вк = V Ъ.Вкт — симметричная матрица коэффициентов обоб-

т=1

щенных жесткостей, имеющая нулевыми две первые строки и два первых столбца; &Вкт — матрица коэффициентов обобщенных жесткостей для т-го конечного упругого элемента; Ny — общее число конечных элементов силовой схемы.

Элементы матрицы ДВкт представляются в виде квадратичной формы от обобщенных смещений fk при (& = /, /) граничных точек то-го конечного элемента

(18)

где R — матрица коэффициентов жесткости конечного элемента. Для элемента балки (см. фиг. 4) принято, что в пределах его

длины I величины —-—--------- --- изменяются линейно. Тогда мат-

£/(5)’ С/р(5)

рица жесткости элемента балки вычисляется по обычным правилам строительной механики и имеет вид

Яб =

82, -*12 0

1 °12 0

®11®22 ~ 8,2 о 1 0 °11®22 ~ ®|2

®33

(19)

где

и,1 ------

822 =

Е Л

1

+

Е /,

Р

12

» I 1

о,., = -------

и А

Е /, Е [, ) 2 ’

О/,

Ч"

+■

0,5

Е1,

1

р і

О/,

/? 2

з

і

~2~

В1и Е12, С/р1, б/р2 — значения жесткостей изгиба и кручения в концевых сечениях элемента балки.

Компоненты векторов fк при к = г, / для элемента балки выражаются в виде

«*) — ?*(•*!, г,)—:'*?*,(*!»

Л* = 9*5 (-«2, *2)“ ?*,(•*!. 2,), /*з =?*„(*!, г,) —?*„(х„22), і

(20)

где л:,, 2и х2, г2 — координаты концевых точек оси элемента балки; (х, 2), <рАя(лг, 2) — производные функции % (х, 2) по направлению оси балки и по нормали к ней соответственно.

Для трапециевидного в плане кессона, боковые стенки которого приняты абсолютно жесткими на сдвиг, матрица жесткости Як имеет вид

Ль

о

— г

- г~‘

7-ї" — 1 1

1 ; — 1 — 1 1

г — г

— г г

(21)

где г = 1/7;

ЗС5во„с(1 + |*) (/ Я, + А /Уз)2

(5В + 8Н) (Л + /) [1,5 (1 + |*) с2 + (й - /)2 +(<*-/)(<*-£) + (<*- 1)2] ’

8в) °н — редуцированные толщины обшивки верхней и нижней стенок кессона; (7 — модуль упругости материала на сдвиг; ц — коэффициент Пуассона; I, I, с, й — линейные размеры проекции кессона на плоскость (см. фиг. 5); Ни Нъ — длины боковых ребер кессона, соответствующих вершинам 1 и 3 трапеции.

Векторы смещений /к при & = г, /

где ?/(-*„,, гт), <ру (хт, гт) — векторы линейных смещений граничных точек кессона для г'-й и у'-й координатных функций.

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

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

Как показывают расчеты, при колебаниях рулей с малыми частотами влияние составляющих С‘о и Ь18 в уравнениях (10) мало, поэтому при выводе упрощенных уравнений короткопериодического движения этими составляющими пренебрегаем.

Деформации конструкции в системе дифференциальных уравнений (10) представляются обобщенными координатами н; (г> 3). В соответствии с предположением о квазистатическом изменении упругих деформаций при короткопериодических колебаниях самолета примем

Тогда система я-линейных дифференциальных уравнений сводится к системе двух дифференциальных уравнений с переменными м, и и2:

1=3

и линейной алгебраической системе (я — 2)-го порядка с восемью слагаемыми в правой части

// = ?/(*„.. гт), /} = 4}(хт, гт) (т = 1, 2, 3, 4), (22)

(23)

П

/ = 1 п

(24)

С22 Но — ^21 "Ь 22 4-2 В

22 ^2

-р В2 й б/в 2 ~Ь В02 4“ ^ В21

(.В — Вк) и =— са' «1 — с“а «2 + /)“■ и, + /)из и> + Ви- иг 4+ В*Ъ + Ов + В0,

где и — (п — 2)-мерный вектор, составленный из обобщенных координат «;(г>-3), (В — Вк) — матрица, образованная из матрицы (В —

— Вк) вычеркиванием двух первых строк и столбцов.

Векторы-столбцы в правой части уравнений (25) формируются из коэффициентов системы (10) по правилам:

Определив из (25) обобщенные координаты и1 (г>-3) и подставив их в правые части дифференциальных уравнений (24), получим уравнения короткопериодического движения самолета с учетом упругих деформаций

Все приращения коэффициентов в системе уравнений (27) являются функциями числа М и скоростного напора <?. Коэффициенты Дс12 и Дс2и определяющие инерционное взаимодействие степеней свободы щ и и2 и зависящие от закона распределения масс по конструкции, не равны друг другу.

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

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

где V — скорость полета; а — угол атаки самолета; д — угол тангажа; ю., —угловая скорость самолета; Дпу — приращение нормальной перегрузки; <7 = 9,81 м/с2 — ускорение свободного падения.

Учитывая соотношения (11) и (16), преобразуем систему (27) к виду

(26)

О в і — О в іт2 , В0 1 = Во І-(_2.

С]] = (ґ)\\ 4~ Д/>И) и1 -(- (/)]2 + Д/)(2) Ч2 4~ (/?|2 "I- Дв,2) ^2 4~

4~ + А#і)3 (Ов 5 -(-Дйв і)-{-(/?0і 4- ДЄ01) — ^си и\ ~ Дсг> и> > ^7)

С22 и2 — (/?21 4~ Д ^2і) М1 4~ (^22 4~ Д^2ї) М2 4" (Во2 4" ^В22) И2 +

4- (Во -|- Д/?г) § + (Св 2 + Д£/в 2) 4~ (В02 4- ДВй2) — Д<?21 щ —Дс22 и >.

«1 = V (& — “)> и2 = &- «2 = “а, «1 = А/гу ё,

(28)

Мс V(<»г — я) = (-0)2 + д^іг) “г 4~ (В\г + ДЄіг) я-)-(Єі 4- Д/?і')84-

4- б?в і + /?0і -|- ДОв і (1 -}- Дпу) -|- ДВ01 — Дс12 шг;

1гс шг — (-^22 + Д^гг) “г 4~ (Є22 4~ ДВ22) <* ~Ь (^2 + Дф>) 0 4" С?В 2 4" 4- -®02 4“ Д0В 2 (1 + ^Яу) 4“ ДД)2 — Д^гг^г-

(29)

9— Ученые записки ЦАГИ № 2

129

В соответствии с правыми частями уравнений (29) приращения коэффициентов подъемной силы Дсу и продольного момента Ат2, вызванные упругими деформациями, представим в форме

ДСУ = ДC^z (q, М) ш2 + ДСу (q, М)а+ Дс°у (q, М) 8 + Дсу0 (q, М) +

+ Дell (q, М) nv + Дс»г (q, М) шг,

Дягг = Дm™2(q, М) -f Дm\{q, М)л-{- Д/ге*(^, Л4)8 -f-

+ Дотг0 (9, А*) + bmj {q, М) пу + Дт“г (?, Ж) ш„ где <7 —скоростной напор; пу — величина нормальной перегрузки самолета; юг = гуА— безразмерная угловая скорость самолета;

— <oz b\

ш =---------безразмерное угловое ускорение; Ьк — средняя аэро-

- уг -

динамическая хорда.

В результате учет влияния упругих деформаций на короткопериодическое движение самолета сводится к вычислению 12 функций двух переменных: скоростного напора q и числа М. Функции Дc*y(q,M), Дc\(q,M), Дcy0(q, М), Дm^q, М), Дm\{q,M), Arnz0(q, М) могут быть также определены по весовым испытаниям упругоподобных моделей в аэродинамических трубах.

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

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

ЛИТЕРАТУРА

1. Белоцерковский С. М., Скрипач Б. К., Табачников В. Г. Крыло в нестационарном потоке газа. М., „Наука*, 197!.

2. Бисплингхофф Р. Л., Эшли X., Халфмэн Р. Аэроупругость. М., Изд. иностр. лит., 1958.

3. Ч у б а р о в В. И. Применение вихревого метода расчета аэродинамических характеристик тонкого крыла произвольной формы в плане в установившемся сверхзвуковом потоке. „Ученые записки ЦАГИ-, т. VI, № 1, 1975.

4. W у k е s J. and Lawrence R. Aerothermoelasticitysome recent studies of the impact on stability and conrol of winged aerospace vehicles.

A1AA Paper, N 64-489, 1964.

Рукопись поступила ///AY 1974 г.