УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том 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 г.