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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том V 1 9 7 4 Мб

УДК 629.735.33.015.4 :533.6.013.43

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

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

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

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

1. Упруго-массовая схематизация самолета. Фюзеляж, крыло, горизонтальное оперение представим балками переменной жесткости, расположенными в одной плоскости, работающими на изгиб и кручение и нагруженными сосредоточенными грузами (фиг. 1). Оси балок крыла и горизонтального оперения в общем случае представляются ломаными линиями, состоящими из прямолинейных отрезков. Чем больше отрезков, тем точнее аппроксимируется форма криволинейной оси жесткости.

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

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

1) V — смещение центра инерции груза 0Г параллельно оси Оу;

2) <р— угол поворота груза вокруг главной центральной оси инерции груза Ог хт (фиг. 2);

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

117

N..

3) 6 —угол поворота груза вокруг главной центральной оси инерции груза Ог zt (см. фиг. 2).

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

В принятой схематизации число упругих элементов равно числу грузов.

Каждый элемент балки вместе со своим грузом (см. фиг. 2) описывается 12 параметрами:

1) Dz— проекция элемента балки на ось Ог (всегда Dz >0);

2) Dx — проекция элемента балки на ось Ох (при Dx > 0);

3) Gr— вес груза;

4), 5) 1Х г, /2Г — главные центральные весовые моменты инерции груза;

6), 7) Е1У, £/2—величины из-гибной жесткости в концевых сечениях элемента балки;

8), 9) G/j, G/2 —величины крутильной жесткости в концевых сечениях элемента балки;

10) а — расстояние от центра инерции груза до оси элемента

О

Яроншгпейн

Г Центр инерции груза.

Фиг. 1

Фиг. 2

балки (о>0 при положительных координатах центра инерции груза в системе осей элемента Оэхэгэу,

11) т — расстояние от проекции центра инерции груза на ось элемента

балки до основания жесткого кронштейна (т^-О, если в системе осей элемента Оэ хэ гэ.координата второго сечения гэ2 меньше координаты центра инерции груза гэ г ); .

12) е — угол между осью Оэхэ и главной центральной осью инерции груза 0ГхГ (е>0, если угол стреловидности оси Ог-*т меньше угла стреловидности оси Оэ г3).

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

Все элементы расчетной схемы пронумерованы по порядку. Первый элемент» расположенный в начале координат, имеет нулевую длину Бгх = Охх == 0. Далее следуют элементы носовой части фюзеляжа с горизонтальным оперением (в случае самолета схемы „утка*), хвостовой части фюзеляжа с горизонтальным оперением (в случае самолетов нормальной схемы) и крыла. Для всех элементов балки фюзеляжа о; — = 0. Для первого груза, кроме того, принято ^1=0. Эле-

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

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

і <

■*< = 2 °хр = 2 °г) (' = 1............... ЛГ12);

У=1 /=1

горизонтальное оперение в носовой части фюзеляжа

N,1 і

*< = 2 2 Вх)'

/=1 /=ЛГ12+1

Л'п *

= 2 Ог2 Е*г/ (г' = -^12 +*>•■•• Л^із);

/=і /=ЛГ12+1

хвостовая часть фюзеляжа

І І

~ ^ = {і = ^13^-1» * •» А^зз);

/=ЛГ 18+1 , /=Л^18+1

горизонтальное оперение в хвостовой части фюзеляжа

хі= 2 °х1+ 2 °хл

/ = Л^ 13+1 /— Л^22+1

ЛГ21 г

'*/ = 2 °*і+ 2 °гі (г' = ^я+1.---.^я);

/ = Л/І3+1 /“/V22 "Ь 1

крыло

і '

Хі= 2 Вхі' гі

І — іїаз+1 , /=Л^3+1

Величины А^ц, Л^,2, Л^13, Л/'з!, Л^22. Л^23. Л^ЗІ. Л^зз являются порядковыми

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

Угол стреловидности г'-го элемента Хэ/ вычисляется через величины Бгі к Бхі

Охі

2 Dzj {і — + 1.......Ni3).

Хэ і = arcsin

V Dx]+Dz)

Если У Dx] + Dzj = 0, то хэ і = Хэ

Для элементов балки носовой части фюзеляжа Хэ/= —%12. для элементов балки хвостовой части фюзеляжа Хэ/ = я/2, для элементов балок крыла и горизонтального оперения

— я/2 < Хэ і < л/2-

Координаты центра инерции г-го груза хГ ¡, zr ¡ и угол стреловидности оси Or i zT і іт і определяются через координаты г-го узла оси жесткости и угол стреловидности элемента по формулам

хг i = xi + a¡ cos хэ / + */ sin хэ/.

¿г і = ¿i - Ч sin Хэ і + У cos хэ і,

Хг і = Хэ і — Ч-

При колебаниях положение упругого самолета в момент времени t будет определяться вектором обобщенных координат размерности n = SNS3

t>i (о

Ті {і)

0i W

S(0 =

*ч,(*>

Ч,<*>

где 1<; (¿), <р; (¿), 6/ (<) — обобщенные координаты г-го груза в моменты времени

(1)

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

Cl 0 0

с = 0 «2 0

0 0 с3

(2)

где С], с2, с3 — квадратные симметричные матрицы коэффициентов влияния балок носовой части фюзеляжа, хвостовой части фюзеляжа и крыла. Порядки пи п2, соответственно равны

Пх = ЗЛ^з, = 3 (ЛЛ.з — ЛГ13), «з = 3 (Л33 — ^зз)- (3)

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

грузов. Стрелками в точках стыков элементов показаны главные центральные оси. инерции грузов.

Для принятых координат матрица коэффициентов влияния разветвленной балки су по методу Аргироса [1] представляется в виде произведения матриц

(4>

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

Обе матрицы Вр Пу порядка Иу состоят из квадратных блоков третьего порядка. Квадратный блок б*; матрицы усилий Ву, соответствующий обобщенным усилиям в к-м стыке при действии обобщенных сил в г'-м стыке (г>£), опреде-

ляется из условия равновесия и имеет вид

1 0 0

Вы = (xí — xk) sin х* + + (*i - Zk) cos Xft cos(x¡— Xk) sin (Xt - ~Xk)

- (x¡ - xk) cos x* + + (*/ - *k) sin Xk - sin (x* -Xk) cos (хг ^ ~Xk)

(4)

где XI, 21, Хк, — координаты стыков элементов,

■ Хг> Хк — Углы стреловидности главных осей инерции Ог/гг(- и Ог*гг*. Так как при к = 1 квадратные блоки Вц являются единичными матрица*!» третьего порядка, а при £>г элементы = 0, то матрица усилий бу является верхней треугольной матрицей.

Для силовой схемы, представленной на фиг. 3, матрица усилий имеет вид

В,

Bn 6] 2 Bi3 Bu Вхь В\ъ Bi7

0 B22 B23 624 Bro #26 627

0 0 B33 B.u Взъ 0 0

0 0 0 Bu В,5 0 0

0 0 0 0 Въь 0 0

0 0 0 0 0 •See B&j

0 0 0 0 0 0 Bn

(5)

Блоки В*;(3^;й<;5, 6 7) в матрице Ву равны нулю, так как силы, дей-

ствующие на стыки 6 и 7, не вызывают усилий в стыках 3, 4 и 5.

Блок Пн третьего порядка матрицы Пу является матрицей податливости куска балки, расположенного между стыками с номерами (* — 1) и і по направлениям обобщенных смещений г-го груза и вычисляется по формуле

(6)

где

COS s - sin Є 8U 8,2

COS є

О

(7)

^21 ^22 О

О 0 833

Величины т, о и е для г-го груза определены выше. Коэффициенты податливости куска балки 8И, 6]2, 822, 533 вычисляются при предположении, что величины

—1— и — вдоль прямолинейного упругого звена изменяются линейно:

ЕЦБ) О/ (5)

Е1Х + £/;

1

Р

12

512 — 621 — í

+

2 El,

JL

З

, 1 1

= (ж +

Eh

°33^(g/1 +

1

G/a

(8)

(9)

■где EIU GIU EI2, G/г — величины изгибной и крутильной жесткостей куска балкй в его концевых сечениях, I — длина куска балки.

Перемножив матрицы в формуле (6), получим коэффициенты матрицы податливости П« в виде

хц = ^11 Ч" 2 6]2 т -[- й22 т2 "Ь 833 а2,

*12 = *21 ” (812 ~Ъ 822 т) сое — 833 а sin £,

*13 = *31 = — (®12 “Ь 822 z) sin £ — ®33 0 COS в,

*22 = 822 “Ь (833 — 822) Sin2 £,

7^23 = *32 = (633 — В22) COS £ Sin £,

*33 = 833 — (533 — 822) sin2 s.

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

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

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

V V г2 0 1 cos т] 1 + + Г2 ф 1 sini) 1 г ¥ sin i¡ г в eos t)

ф r<f sin r¡ <p 0

0 rd cos r¡ 0 0

P

• Mx , (10)

Mz

где V, ср, 0 — обобщенные координаты; р, Мх, Мг — обобщенные силы; V, ?, 6 — коэффициенты податливости; г — расстояние от центра инерции груза до центра вращения груза при парциальных колебаниях (всегда 0) (фиг. 4); -(¡ — угол между вектором г и главной центральной осью инерции Огхг.

Величина 1) может принимать только четыре значения —я/2, 0, я/2, к. Различные положения вектора г представлены на фиг. 5.

4. Уравнения колебаний самолета в пустоте, итерационный метод расчета низших тонов собственных симметричных колебаний. При свободных колебаниях незакрепленного самолета вектор деформаций 5Д(£) в системе координат Охуг, связанной с первым элементом схемы (см. фиг. 1), запишем в виде

Sí(t) = S(t)~'£íaí(t)sl,

(П>

¡ = 1

где 5 (<) — вектор обобщенных координат (1); 5,-нат, соответствующий поступательному смещению самолета как твердого тела, параллельному оси Оу; 52 — вектор обобщенных координат, соответствующий повороту вокруг оси Осгс, параллельной оси Ог и проходящей через центр тяжести самолета; <*¿(7) — функции времени.

•вектор обобщенных коорди-

Центр инерции^ £

(т,1)

Центр дращения

Фиг, 4

N r(7¡ = rr) in-f) >ч <-=* /%

°r \ г>!?

Центр инерции

Г (7¡-0) груза

X,

Фиг. 5

Векторы 51( соответствующие собственным колебаниям с нулевыми частотами, отнормированы так, что

STAÍS/=1 (/=1,2),

(12>

где М — диагональная матрица удвоенных масс и моментов инерции. ■

При известной матрице коэффициентов влияния с (2) напишем уравнения симметричных свободных колебаний конструкции в пустоте:

5 (<) - £ (t) S¡ = А [— MS (<)],

í = 1

(¡3>

где S(¿)— вектор обобщенных ускорений; А — матрица коэффициентов влияния, элементы которой вдвое меньше элементов матрицы с.

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

S (¿)TAfS(-= 0 (/ = 1, 2), (14)

исключим из уравнения (13) функции a¡(t):

2

S (t) = - AMS (It) + £{(/Ш5'(0]т MS i) S¡. (I5>

¿=i

Представим закон движения при собственных колебаниях в виде

S(¿) = Ssinco¿, (16)

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

После подстановки (16) в дифференциальное уравнение (15) получим характеристическое уравнение

2

/.5 = АМБ - £ [(;Ш5)Т уМ5(] 5/,

(17)

<=]

где

Введя преобразование координат по формуле

5 = М,/2 5,

где

Д*1/2

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

(18)

1-і

где матрица В = М АМ ,как и матрица А, симметричная. Векторы ¿’¿(г=1,2) единичной длины и попарно ортогональные.

Уравнение (18) имеет п — 2 положительных собственных значений

^ (( — 3....и), которым соответствуют попарно ортогональные собственные

векторы 5;. Будем считать, что все собственные числа различны

(19)

Обозначим матричный оператор правой части уравнения (18) через Ь (Б) и воспользуемся при расчете А3 и 53 методом последовательных приближений [3].

з

(20)

где к — номер приближения.

В качестве начального вектора ^ можно взять вектор

(21)

Итерационный процесс (20) в силу неравенств (19) сходится к наибольшему собственному значению Л3 и собственному вектору 53.

Сходимость оценивается по норме

тах|(53А+1),-(5^|,

(22)

где (5|)у — у'-й элемент вектора 5* .

При вычислении последующих собственных значений и собственных векторов воспользуемся ортогональностью собственных векторов. В результате для расчета /-го тона (У = 4, 5,...) оператор 1(5) в формулах итерационного процесса (20) нужно заменить на оператор

І, (5) = і ($)-£[£ (5)т5/]5г,

(23)

г=з

где 5; — собственные векторы с номерами г' = 3, 4.....(у—1).

_ В соответствии с алгоритмом расчета все вычисленные собственные векторы 5/, включая нулевые, имеют единичную длину и попарно ортогональны

Каждому собственному значению А; (г = 3, 4,...) и его собственному вектору 5,- соответствует собственное колебание с частотой юг-= 1/]^Хг- и формой

1. Аргирос Дж. Современные методы расчета сложных статически неопределимых систем. М., Судпромгиз, 1971.

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

3. Д е м и д о в и ч Б. П., М а р о н И. А. Основы вычислительной математики. М., Физматгиз, 1963.

4. Кандидов В. П., К и м Л. П. Дискретная модель для исследования колебаний балок, Вестник МГУ, сер. „Физика и астрономия“, № 5, 1968.

і

¿¡ = М 2 5|.

ЛИТЕРАТУРА

Рукопись поступила 17/1 1973 г.