Документ взят из кэша поисковой машины. Адрес оригинального документа : http://belyakov.imec.msu.ru/papers/Diser.pdf Дата изменения: Mon Sep 8 14:32:51 2008 Дата индексирования: Mon Oct 1 19:36:38 2012 Кодировка: Windows-1251

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ им. М.В. ЛОМОНОСОВА Механикоматематический факультет На правах рукописи

Беляков Антон Олегович

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

Диссертация на соискание ученой степени кандидата физикоматематических наук

Научный руководитель доктор физикоматематических наук А.П. Сейранян

Москва 2005 г.


Содержание
Введение Глава I. Описание метода измерений 1 Существующие методы измерения моментов инерции 2 Конструкция измерительного стенда и процесс измерений 3 Математическая модель стенда и постановка задачи Глава II. Способы возбуждения колебаний и алгоритмы определения моментов инерции 1 Возбуждение известной силой, приложенной до начала движения 2 Нахождение моментов инерции без информации о способе возбуждения колебаний
2.1 2.2 Пружины одинаковы, известна масса . . . . . . . . . . . . Известна масса и положение центра масс тела . . . . . . .

5 11

11

19 22

32

36

39
39 40

3 Идентификация линейной колебательной системы
2

42


3.1

Определение частот, декрементов затухания и амплитуд сигнала методом Прони . . . . . . . . . . . . . . . . . . . . . . 43 47

3.2

Идентификация методами пространства состояний . . . . .

4 Оценка минимального времени возбуждения многомерной колебательной системы с управлением
4.1 4.2 Синтез управления . . . . . . . . . . . . . . . . . . . . . . . Синтез управления при малом управлении и неизвестных параметрах системы . . . . . . . . . . . . . . . . . . . . . . 4.3 4.4 Максимальное время возбуждения системы . . . . . . . . . Оценка минимального времени возбуждения многомерной колебательной системы . . . . . . . . . . . . . . . . . . . . 4.5 Пример расчета . . . . . . . . . . . . . . . . . . . . . . . . . 68 70 64 64

55
59

Глава III. Численный эксперимент и анализ погрешностей 1 Численное моделирование процесса колебаний 2 Пример вычислений динамических параметров тела 3 Анализ погрешности определения динамических параметров
3.1 Анализ чувствительности динамических параметров к ошибке идентификации . . . . . . . . . . . . . . . . . . . . . . .

73 73

78

82
82

3


3.2

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

4 Влияние демпфирования на частоты и формы колебаний системы Выводы Список литературы Приложение 90 98 100 106

4


Введение
Актуальность темы
Для обеспечения требуемых маневренных характеристик самолетов и быстроходных морских судов конструкторам требуется знать моменты инерции их массивных деталей. Но из-за сложности конструкции некоторых элементов, таких как силовые установки, аналитически определить их моменты инерции не представляется возможным. Возникает задача измерения моментов инерции массивных крупногабаритных тел. В Центральный аэрогидродинамический институт им. проф. Н. Е. Жуковского (ЦАГИ) поступил заказ от Калужского турбинного завода на разработку стенда для измерения моментов инерции крупногабаритных тел. По мнению специалистов ЦАГИ существующие методы измерения моментов инерции [1, 2, 3, 5, 6, 7] или трудно применимы для крупногабаритных тел или не дают требуемой точности. В. В. Богданов [9] предложил новый метод измерения моментов инерции, где помещенное на четыре пружины тело совершает свободные колебания по трем степеням свободы. Моменты инерции тела (включая центробежные) содержатся в инерционной матрице колебательной системы. Инерционная матрица определяется по сигналу с датчиков, измеряющих смещение тела. При этом не требуется знание жесткостей пружин. На основе предложенного метода в ЦАГИ спроектирован стенд для измерения моментов инерции крупногабаритных тел [9].

Научная новизна
5


Метод измерения моментов инерции, на основе которого в ЦАГИ разработан стенд [9], является новым, следовательно, возникла необходимость разработки соответствующих математических алгоритмов. Настоящая диссертация посвящена созданию этих алгоритмов и анализу погрешности определения моментов инерции тела. В работе показано, что разработанный математический аппарат совместим с различными методами идентификации линейных систем [13, 14, 15]. Представлены три варианта решения задачи определения моментов инерции тела в зависимости от сведений о способе возбуждения колебаний или о параметрах тела и жесткостях пружин стенда: 1. Начальное смещение системы от положения равновесия вызывается при помощи известной силы. При этом разработанный автором алгоритм определяет не только моменты инерции тела, но также его массу и положение центра масс. Все динамические параметры тела определяются без информации о жесткостях пружин, кроме предположения о постоянстве жесткостей пружин в процессе колебаний. Другими словами, жесткости пружин косвенно определяются в процессе измерений. Таким образом учитывается изменение жесткостей пружин при их деформации под весом тела, так как эта деформация на порядок больше амплитуды колебаний в процессе измерений. Это избавляет от необходимости производить калибровку пружин перед измерениями. 2. Способ приведения системы в движение неизвестен, но известна мас6

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

Практическая значимость
Полученные алгоритмы предполагается использовать на разработанном в ЦАГИ стенде для определения моментов инерции крупногабаритных тел. Данные автором рекомендации могут быть использованы при конструировании аналогичных измерительных стендов.

Метод исследования
При разработке алгоритмов определения моментов инерции использовались методы теории колебаний [23, 38, 39, 40, 41, 42], теории идентифика7


ции линейных систем [13, 14, 15, 18, 19, 21], теории оптимального управления и теории возмущения матричных операторов.

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

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

8

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

Апробация
По теме диссертации подготовлены публикации [30, 31, 32, 33]. Основные результаты были доложены на следующих конференциях: 1. Восьмой Всероссийский съезд по теоретической и прикладной механике. Пермь, 2001. Численное моделирование процесса измерения моментов инерции крупногабаритных тел методом свободных колебаний. 2. Научная конференция МФТИ, 2001. Определение инерционной матрицы по формам свободных колебаний. 3. Международная конференция Математические идеи П. Л. Чебышева и их приложения к современным проблемам естествознания. Обнинск, 2002. Моделирование процесса измерения динамических параметров массивных тел по формам колебаний.

9


4. Научная конференция МФТИ, 2002. О способах определения моментов инерции по упругим колебаниям, (первая премия). 5. Вторая международная научно-техническая конференция молодых специалистов Современные проблемы аэрокосмической науки и техники. Жуковский, 2002. Определение инерционных характеристик крупногабаритных тел по колебаниям в упругом подвесе. 6. Международная научная конференция по механике Третьи Поляховские Чтения. С.-Петербург, 2003. Способы определения динамических параметров тел по колебаниям в упругом подвесе. 7. Международная конференция Физика и управление. С.-Петербург, 2003. Optimal excitation of oscillations by a limited control force. Автор выражает глубокую признательность своему научному руководителю А. П. Сейраняну за руководство и за переданное умение самостоятельно отличать полезный научный результат. Автор благодарит Ю. В. Болотина, И. В. Новожилова, И. Л. Антонова, В. В. Александрова, Ю. Г. Мартыненко за ценные замечания по тексту диссертации, а также О. Н. Кириллова и А. А. Майлыбаева за полезные советы. Автор особо благодарит В. В. Богданова за постановку задачи и консультации и Л. Ю. Блаженнову-Микулич за совместную работу и за указание на литературу по методам пространства состояний.

10


Глава I. Описание метода измерений 1 Существующие методы измерения моментов инерции
Все методы экспериментального определения моментов инерции используют вращение с ускорением или крутильные колебания тел. Конструкции измерительных устройств можно разделить на две группы: конструкции, фиксирующие ось вращения, и конструкции, не фиксирующие ось вращения. Простой пример конструкции с фиксированной осью враще-

Рис. 1: Определение момента инерции махового колеса

ния показан на рисунке 1 из [5]. Здесь момент инерции махового колеса определяется по времени, за которое опустится груз, прикрепленный к тросу, намотанному на это колесо с радиусом R. Для учета диссипации измерения производят с двумя грузами разных масс M1 и M2 , которые
11


опускаются с высоты h за времена T1 и T2 . Получается следующее выражение для определения момента инерции [5]

I=R

g 2h 2

(M1 - M2 ) -
1 2 T1

M1 2 T1

-

M2 2 T2

-

1 2 T2

,

где g ускорение свободного падения. В современных устройствах такого типа исследуемое тело закрепляется на валу электродвигателя [7]. Момент инерции определяется по работе

A1 электродвигателя при разгоне до угловой скорости 1 и по работе A

2

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

I=

A1 - A2 2 2 - J, 1 - 0

где J приведенный момент устройства. Очень точным считается метод измерения на монофилярном подвесе [2], конструкция которого представляет из себя крутильные весы. Тело висит на упругом стержне и совершает крутильные колебания вокруг его оси. Другой простой способ определения моментов инерции тел метод физического маятника (рис. 2). Известно расстояние r от оси вращения до центра масс тела, масса тела m и величина ускорения свободного падения

g . Момент инерции относительно оси, проходящей через центр масс и
12


параллельной оси вращения, находится по формуле

IC = m

gr - m r2, 2

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

c l2 I = 2 - m d2 - J,
где m масса тела, d расстояние от центра масс тела до оси платформы, l расстояние от оси платформы до пружины, J момент инерции платформы, c суммарная жесткость пружин, измеренная частота малых колебаний платформы с телом. К группе конструкций с фиксированными осями относятся также разнообразные варианты подвесов на нерастяжимых тросах [3, 4], где тросы расположены под углом друг к другу. Например, трифилярный подвес (рис. 3) и так называемые мультифилярные подвесы (рис. 5). Схемы измерительных устройств могут не содержать колебательных контуров. Например в [8], синусоидальные колебания тела на платформе производятся в одной плоскости тремя гидравлическими силовозбудителями, (рис. 5). Масса, положение центра масс и момент инерции тела определяется исходя из измеренных датчиками сил, ускорений и смещений в точках крепления силовозбудителей к платформе. К сожалению,
13


устройство, представленное в [8], обеспечивает относительную погрешность определения момента инерции лишь в пределах 10%. В качестве примера конструкции, не фиксирующей ось вращения, можно привести бифилярный подвес, в котором тросы расположены параллельно друг к другу. Если ось вращения находится далеко от центра масс тела (рис. 2, 3, 4, 7), то это вызывает большую погрешность при вычислении момента инерции относительно оси, проходящей через центр масс. Из теоремы Гюйгенса-Штейнера следует [25]

2 = 2 - r2 , m
где m радиус инерции относительно оси, проходящей через центр масс тела, радиус инерции относительно оси вращения, r расстояние от центра масс до оси вращения. При большом r в правой части формулы возникает разность больших чисел, из-за чего относительная погрешность

m становится большой, даже если и r измерены с небольшими относительными погрешностями. Одни конструкции измерительных устройств имеют одну степень свободы, другие несколько, что предпочтительнее, если целью является определение всего тензора инерции тела. Чтобы вычислить тензор инерции тела, используя установку с одной степенью свободы, необходимо определить моменты инерции относительно как минимум шести осей, так как тензор инерции определяется шестью величинами. Эти величины связаны с определяемыми моментами инерции линейными уравнениями типа
14


Рис. 2: Физический маятник.

[25]

Iu = Ix 2 + Iy 2 + Iz 2 - 2Ixy - 2Ixz - 2Iyz ,
где Iu момент инерции относительно оси u,

, , косинусы углов, образуемых осью u с осями Ox, Oy , Oz , Ix , Iy , Iz осевые моменты инерции, Ixy , Ixz , Iyz центробежные моменты инерции. Осевые и центробежные
моменты инерции составляют тензор инерции тела относительно точки пересечения осей. Можно сделать вывод, что для определения тензора инерции тел хорошо подошел бы метод, в котором крутильные колебания тела происходят по нескольким степеням свободы вокруг осей, проходящих вблизи центра масс.

15


Рис. 3: Трифилярный подвес.

Рис. 4: Мультифилярный подвес.

16


Рис. 5: Определение момента инерции по колебаниям тела в плоскости.

Рис. 6: Бифилярный подвес.

17


I r J mg с

l

Рис. 7: Определение момента инерции на качающейся платформе.

18


2 Конструкция измерительного стенда и процесс измерений
В ЦАГИ В. В. Богданов [9] предложил метод измерения моментов инерции, в котором колебательная система создается посредством упругих элементов, при этом колебания происходят по трем степеням свободы и две оси крутильных колебаний проходят недалеко от центра масс. Как сказано в первой главе, эти обстоятельства являются преимуществами данного метода. Принципиальная схема спроектированного в ЦАГИ стенда для измерения моментов инерции крупногабаритных тел показана на рисунке 8. Стенд предназначен для измерения инерционных параметров тел цилиндрической формы массой порядка 15т, длиной 6м и диаметром 2м. Тело при помощи специальных хомутов помещается на шарнирах на четыре пружины, закрепленные консольно другими концами на бетонном основании. Изгибные жесткости пружин обозначены горизонтальными пружинами на рисунке 8. Точки крепления пружин к хомутам образуют прямоугольник, также как точки крепления пружин к основанию. В точках крепления пружин к хомутам расположены датчики, измеряющие вертикальное смещение. В системе возбуждаются малые колебания по трем степеням свободы. В результате определяются два осевых и один центробежный моменты инерции тела. Затем измерения повторяются с телом, повернутым на 900 вокруг оси C x (рис. 8). Таким образом, определяется весь тензор инерции тела за исключением одного центробежного момента инерции. Другими словами, кроме значения главных моментов
19


Рис. 8: Принципиальная схема измерительного стенда.

20

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

21


3 Математическая модель стенда и постановка задачи
Введем связанную с телом ортогональную систему координат Oxy z , рис. 8 и 9. Точка O начала координат находится в точке крепления первой пружины к телу. Ось Ox проходит через точку крепления второй пружины к телу. Ось Oz проходит через точку крепления четвертой пружины к телу. Ось Oy направлена вверх перпендикулярно осям Ox и Oz . Аналогично введем неподвижную систему отсчета Oa X Y Z (рис. 8), связав е? с точками крепления соответствующих пружин к основанию стенда. Запишем уравнения движения твердого тела относительно центра масс. Для этого введем систему координат C xy z (рис. 8), начало которой расположено в центре масс тела C , а оси параллельны осям системы Oxy z . Пусть система координат Ca X Y Z неподвижна относительно системы Oa X Y Z и если тело находится в положении равновесия, то система

Ca X Y Z и связанная с телом система C xy z совпадают. Движение системы описывается законами изменения импульса и момента импульса, записанными в проекциях на оси системы C xy z

I

dw + w Ч Iw = dt d2 x m 2= dt

i Ч f
i

i

(1) (2)

f
i

i

где I тензор инерции тела в системе C xy z , w вектор угловой скорости вращения системы C xy z относительно системы Ca X Y Z в проекциях на оси системы C xy z , x вектор смещения центра масс тела от положения равновесия,
d2 x dt2

- вектор ускорения центра масс, i радиус-векторы точек
22


приложения сил fi . К уравнениям движения (1) и (2) следует добавить уравнение наблюдений

y = h(1 , . . . , 4 , x),

(3)

где h вообще говоря, нелинейная вектор-функция, y вектор вертикальных смещений точек крепления пружин относительно положения равновесия, получаемых с 4-х датчиков. Предполагается, что в процессе измерений смещение точек крепления пружин к телу относительно положения равновесия мало по сравнению с характерными размерами стенда и, следовательно, амплитуды крутильных и поступательных мод колебаний также малы. Предполагается, что скорости центра масс и угловые скорости тела малы по сравнению со своими единицами размерности. Легко показать [10], что в уравнении (1) с точностью до малых величин второго порядка можно отбросить член

w Ч I w. Рассмотрим правые части уравнений (1) и (2), оставляя там только члены первого порядка малости. Каждая пружина создает усилие, направленное противоположно смещению точки крепления к телу. При этом не создается дополнительного момента относительно осей, проходящих через точку крепления пружины к телу, так как крепление шарнирное. Ввиду того, что пружины имеют цилиндрическую форму, они обладают тремя взаимно перпендикулярными плоскостями симметрии упругих свойств [50]. Поэтому в системе координат Oa X Y Z , совпадающей с осями симметрии пружины, упругие постоянные пружины можно представить

23


в виде матрицы




(4)

c00 i Gi = 0 ki 0 , 0 0 ci

где i порядковый номер пружины, ki жесткость пружины в продольном (вертикальном) направлении, а ci жесткость пружины при смещении е? свободного конца в любом поперечном направлении. По диагонали матрицы Gi расположены коэффициенты жесткости пружины в направлениях Oa X , Oa Y и Oa Z . Эти коэффициенты считаются постоянными на всем протяжении процесса колебаний. Действующие на тело со стороны пружин силы выражаются по следующей формуле

fi = -G i , i

i = 1, 2, 3, 4,

(5)

где i вектор смещения точки крепления i-ой пружины от положения равновесия, a G матрица жесткости i-ой пружины в связанной с телом i системе C xy z . Эта матрица связана с матрицей Gi следующим образом

G = ST ST Gi SO S, O i

i = 1, 2, 3, 4,

(6)

где S матрица преобразования поворота системы C xy z относительно системы Ca X Y Z , а S
O

системы Ca X Y Z относительно системы

Oa X Y Z . Наклон системы Ca X Y Z относительно системы Oa X Y Z также будем считать малым. Положение тела в инерциальной системе отсчета Ca X Y Z описывается радиус-вектором положения центра масс x и тремя углами поворотов

x , y и z связанной с телом системы координат C xy z относительно осей
24


системы Ca X Y Z . Углы x , y и z называются углами Крылова [11]. Они удобны тем, что из малых углов Крылова можно составить вектор малого поворота , взаимнооднозначно связанный с матрицей преобразования поворота S. Матрица S состоит из произведения матриц поворотов

S = Sx Sy Sz относительно оси Ca X 1 0 0 Sx = 0 cos(x ) - sin(x ) , 0 sin(x ) cos(x )
оси Ca Y

Sy = sin(y ) 0 0 cos(y ) 1 0



- sin(y ) 0 cos(y )
и оси Ca Z



cos(z ) - sin(z ) 0 Sz = sin(z ) cos(z ) 0 . 0 0 1
Выразим смещения пружин i через векторы x и

i = x + Si - i ,

i = 1, 2, 3, 4.

(7)

Углы x , y и z малы, следовательно, матрица поворота может быть записана в виде

1 S = z - - 1
y z

-x . 1
y

(8)


25

x


Аналогично запишется матрица поворота системы Ca X Y Z относительно системы Oa X Y Z на малый угол

1 - Z Y SO = Z 1 -X -Y X 1

.

Из формулы (6) видно, что с точностью до первого порядка малости углов

и G = Gi , i i = 1, 2, 3, 4.

Теперь уравнение (5) можно записать следующим образом

fi = -Gi i ,
После введения обозначения

i = 1, 2, 3, 4.

(9)

0 = S - E = z -
выражение (7) примет вид

- 0
y z



x

-x 0



y

(10)

i = x + i ,

i = 1, 2, 3, 4.

(11)

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

ћ = Ч = - Ч = -R ћ ,
26

(12)


где матрица R имеет структуру, аналогичную матрице





0 -z y R = z 0 - -y x 0

x

.

(13)

Отсюда выражение (11) запишется линейно относительно векторов x и

i = x - Ri .

(14)

После подстановки полученного выражения в (9) получаем выражение для сил

fi = -Gi x + Gi Ri ,

i = 1, 2, 3, 4.

(15)

Учитывая, что i Ч fi = Ri fi , получим выражения для моментов сил

i Ч fi = -Ri Gi x + Ri Gi Ri ,

i = 1, 2, 3, 4.

(16)

Уравнения движения (1) и (2) теперь можно записать в матричном виде













mE 0 0 I

Е x Е

+

Gi

-

Gi R

i i

x

= 0,

Ri Gi -

Ri Gi R

(17)

где E единичная матрица. Заметим, что так как матрицы Ri кососимметрические, а матрицы Gi симметрические, то
T

Ri Gi = -
i i

Gi R

i

.

(18)

Легко показать, что матрица жесткости представима в следующем виде





G

i

-

Gi Ri Ri Gi Ri
27

= CT GC,

Ri Gi -

(19)


где

G 1 0 G= 0 0 0 G
2

0 0 , G3 0 0 G4 -R1 -R2 . -R3 -R4 0 0

0 0

(20)

E E C= E E

(21)

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

стемы. Это точка, в которой приложенная сила вызывает только поступательное смещение тела, то есть = 0. В частности, если точка C была бы центром жесткости, то ввиду (18)
T

Ri Gi = -
i i

Gi R

i

= 0.

Выразим радиус-векторы точек крепления пружин к телу i в системе

C xy z через положение центра масс r в системе Oxy z , используя геометрию стенда (см. рис. 9),

1 = -r, 2 = (Lx , 0, 0)T - r, 3 = (Lx , 0, Lz )T - r,
28


4 = (0, 0, Lz )T - r,
где Lx и Lz базы стенда, т. е. расстояния между точками крепления пружин по осям Ox и Oz .

Матрица C примет вид = C 100 010 001 100 010 001 100 010 001 100 0 r
z y

-r 0 r
x z z

-r 0 r
z

-r 0
y

-r 0 r
z

-Lx + r -r 0
z

x

-r 0

y

-Lx + r Lz - r 0 r
x z

x

0 1 0 -Lz + rz 001 -r
y

-rx 0 ry Lx - rx 0 . ry Lx - rx 0 ry -rx 0

r

y

(22)

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


тремя степенями свободы, показанными на рисунке 8. Это предположение будет проверено при помощи компьютерной модели стенда в третьей

главе. Тогда матрица G преобразуется в матрицу G размерности 4 Ч 4,
имеющую следующий вид

k 1 0 G= 0 0 0 0 0

k2 0 0 0 k3 0 0 0 k4

В матрицу C размерности 4 Ч 3 войдут столбцы 2, 6, 4 и строки 2,

5, 8, 11 матрицы C 1 -rx 1 Lx - r C= 1 Lx - r 1 -rx r
z x x

rz . -Lz + rz -Lz + rz

(23)

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

Mq + Kq = 0 Е , y = Cq
30

q(0) = q0 ,

q(0) = v0 ,

(24)


где q вектор размерности 3, первой компонентой которого является вертикальное смещение центра масс тела, второй и третьей углы поворота тела вокруг осей Ca Z и Ca X (рис. 8), y наблюдаемый вектор вертикальных смещений, получаемых с 4-х датчиков, M матрица инерции, K матрица жесткости, C матрица наблюдения размерности 4 Ч 3. Матрицы M и K симметрические и положительно определенные размерности

3 Ч 3.
Запишем такое же как (19) выражение для матрицы жесткости

K = CT GC.

(25)

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

m0 0 M = 0 Izz -Izx 0 -Ixz Ixx

,
(26)

где m масса тела, а Iz z , Ixx , Ixz , Iz x элементы тензора инерции тела. Это легко показать, выбрав вторые, шестые и четвертые столбцы и строки из инерционной матрицы системы (17). Задача состоит в определении матрицы M по последовательности значений вектора наблюдаемого сигнала y.

31


Глава II. Способы возбуждения колебаний и алгоритмы определения моментов инерции
В этой главе будут рассмотрены три разные постановки задачи, позволяющие определить моменты инерции тела по его колебаниям в упругом подвесе. Вначале получим выражение для матрицы инерции M через матрицы

K иA =M

-1

K , где матрицы K и M - матрицы жесткости и инерции

системы (24), записанной в нормальных координатах

Mq +Kq = 0 Е , y = Cq

q (0) = q0 ,

q (0) = v0 .

(27)

В нормальных координатах матрицы K и M диагональны по определению. Элементами диагональной матрицы A являются квадраты частот колебаний системы. Матрица S
-1

преобразования координат системы (24)

к нормальным координатам (27) состоит из собственных векторов матрицы

A = M-1 K.
Следовательно выполняется следующее соотношение

(28)

A = S-1 A S.
Заметим, что матрица A может быть диагональной не только в нормальных координатах системы (27), (например, если в (28) M = K, где K не диагональна). Известно [23], что такое возможно только в случае
32

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

M = ST M S.

(29)

Матрицу преобразования S, переводящую систему (24) из нормальных координат в исходные физические координаты, можно определить если известны матрицы наблюдения C и C систем (24) и (27). Они связаны следующим образом

C = C S.

(30)

Так как матрицы C и C имеют размерность 4 Ч 3, выразим матрицу S из уравнения (30) по методу наименьших квадратов

S = C + C,
где C
+

(31)

= (C T C )-1 C T псевдообратная матрица. Подставив выражение

(31) в формулу преобразования матрицы инерции (29) и учтя, что M =

KA

-1

, получим выражение для матрицы M в следующем виде

M = CT C

+T

KA

-1

C + C.

(32)

33


Это основное выражение, из которого будут получены три алгоритма определения моментов инерции тела. Полагаем, что оценки матриц A и

C будут получены в результате процедуры идентификации, базы стенда Lx и Lz , которые входят в матрицу C согласно (23), известны заранее.
Из выражения (32), учитывая структуру матрицы инерции (26), получим шесть уравнений

m 0 0 I zz I xx Ixz

= M11 (K11 , K22 , K33 ), = M12 (rx , rz , K11 , K22 , K33 ), = M13 (rx , rz , K11 , K22 , K33 ), = M22 (rx , rz , K11 , K22 , K33 ), = M33 (rx , rz , K11 , K22 , K33 ), = -M23 (rx , rz , K11 , K22 , K33 ),
(33)

где неизвестными считаются координаты центра масс rx и rz , диагональные элементы матрицы жесткости в нормальных координатах K11 , K22 , K33 , масса тела m, два главных момента инерции тела Iz z , I
xx

и один центро-

бежный момент инерции Ixz . Всего 9 неизвестных и 6 уравнений. Чтобы решить эту систему и получить моменты инерции, необходимы еще три независимых условия. Заметим, что благодаря структуре (23) матрицы наблюдения C, в первое уравнение системы (33) координаты центра масс rx и rz не входят. Во второе и третье уравнения rx и rz входят линейно, а в остальные три уравнения rx и rz входят квадратично. Это легко показать, рассматривая покомпонентно произведения матриц в выражении (32). Недостающие три независимых условия можно найти из дополнитель34

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

~~~ ~ ки матриц A, C и векторов q0 , v0 . Обозначим их A, C, q0 и v0 . Эти
оценки могут оказаться комплексными. Базис, в котором получены оценки параметров системы, вообще говоря, неизвестен.

~ 2. Найдя преобразование T, диагонализирующее матрицу A, получаем
оценку матриц A и C

~ ~ A = T-1 A T, ~ ~ C = CT ~ ~ 3. Заменяем комплексные элементы матриц A и C их модулями. Обоснование этого действия дано в 4-ом пункте третьей главы.

~ 4. Проверяем, нет ли кратных диагональных элементов матрицы A с
фиксированной относительной погрешностью. Если таковые имеются, прекращаем выполнение алгоритма, выдав сообщение об ошибке. Далее алгоритмы определения моментов инерции отличаются дополнительными данными для решения системы уравнений (33). Получим выражения для заключительного пятого этапа алгоритмов.
35


1 Возбуждение известной силой, приложенной до начала движения
Получим дополнительные условия для нахождения матрицы инерции, возбуждая колебания определенным способом. Пусть перед началом движения в точке крепления к телу одной из пружин приложена известная сила p в вертикальном направлении (см. рис. 10). При этом тело не движется. Затем силу мгновенно убирают, и тело начинает колебательное движение по трем степеням свободы. Пусть f = (0, -p, 0, 0)
T

век-

тор значений сил, приложенных в вертикальном направлении в точках крепления пружин к телу. Тогда вектор обобщенных сил получается из умножения некой матрицы размерности 3 Ч 4 на вектор f . Эта матрица зависит от точек приложения сил и положения центра масс. Она совпадает с транспонированной матрицей наблюдения CT , так как точки приложения сил совпадают с точками, в которых находятся датчики смещения. Примем за дополнительные условия к системе уравнений (33) векторное соотношение в нормальных координатах, связывающее начальное смещение тела от положения равновесия q0 с приложенной перед началом движения силой C T f

C T f = K q0 .

(34)

Из векторного уравнения (34) находим элементы диагональной матрицы K покомпонентным делением вектора обобщенных сил C T f на вектор начального смещения q0 . Затем подставляем диагональные элементы

36


r
O

x

1 r
z

x

2 L
C x
z

z

4 L

z
x

3

Рис. 9: Пружины и датчики. Вид сверху.

y Z O X

200кг

p

Рис. 10: Схема приложения силы.

37


матрицы K в систему (33) и получаем следующую систему уравнений

m 0 0 I zz I xx Ixz

= M11 , = M12 (rx , rz ), = M13 (rx , rz ), = M22 (rx , rz ), = M33 (rx , rz ), = -M23 (rx , rz ),
(35)

в которой на 6 уравнений имеется 6 неизвестных. Опишем решение системы (35). Первое уравнение системы (35) дает нам массу тела. Из второго и третьего уравнений

0 = M12 (rx , rz ) 0 = M (r , r ) 13 x z

выражаются координаты центра масс, так как rx и rz входят в эти уравнения линейно. Полученные выражения подставляются в последние три уравнения для определения моментов инерции. Явные выражения для определения m, rx , rz , Iz z , Ixx и Ixz приводятся в приложении 1. Выражения были получены с использованием пакета MAPLE. Последним пятым этапом алгоритма является вычисление этих выражений, используя оценки матриц A и C , а также заранее известные значения баз стенда Lx , Lz и приложенной силы p. Отметим, что этот алгоритм, кроме моментов инерции, дает массу и положение центра масс тела.

38


2 Нахождение моментов инерции без информации о способе возбуждения колебаний
Рассмотрим варианты нахождения моментов инерции, когда способ приведения системы в движение неизвестен.

2.1 Пружины одинаковы, известна масса
Получим дополнительные условия для нахождения моментов инерции, зная массу тела m и предполагая, что все четыре пружины имеют одинаковую жесткость k . Матрица жесткости в нормальных координатах выражается через жесткости пружин k1 , k2 , k3 , k4 аналогично (25)

K = C T GC ,
где G = diag(k1 , k2 , k3 , k4 ), C матрица, совпадающая с матрицей наблюдения, так как пружины крепятся к телу там же, где расположены датчики. Если жесткости одинаковы и равны k , то

K = kC T C .
Подставив это выражение в (32), получим систему уравнений

m 0 0 I zz I xx Ixz

= M11 (k ), = M12 (k , rx , rz ), = M13 (k , rx , rz ), = M22 (k , rx , rz ), = M33 (k , rx , rz ), = -M23 (k , rx , rz ),
39

(36)


куда неизвестная жесткость пружин k входит линейно. Опишем решение системы (36). Первое уравнение системы (36)

m = M11 (k )
дает нам k , так как масса тела m заранее известна. Выражение для k подставляется в остальные уравнения. Из второго и третьего уравнений

0 = M12 (k , rx , rz ) 0 = M (k , r , r ) 13 x z

выражаются координаты центра масс rx и rz , так как rx и rz входят в эти уравнения линейно. Полученные выражения подставляются в последние три уравнения для определения моментов инерции. Последним пятым этапом алгоритма является вычисление выражений для rx , rz , Iz z , Ixx и Ixz , используя оценки матриц A и C , а также заранее известные значения баз стенда Lx , Lz и массу тела m. Отметим, что данный алгоритм, кроме моментов инерции, дает положение центра масс тела и не использует информацию о способе возбуждения колебаний.

2.2 Известна масса и положение центра масс тела
Получим дополнительные условия для нахождения матрицы инерции из знания массы тела m и координат центра масс rx и rz .

40


Система уравнений (33) примет следующий вид

m 0 0 I zz I xx Ixz

= M11 (K11 , K22 , K33 ), = M12 (K11 , K22 , K33 ), = M13 (K11 , K22 , K33 ), = M22 (K11 , K22 , K33 ), = M33 (K11 , K22 , K33 ), = -M23 (K11 , K22 , K33 ).
xz

(37)

Шесть неизвестных K11 , K22 , K33 , Iz z , Ixx и I линейно.

входят в уравнения (37)

Последним пятым этапом алгоритма является вычисление выражений для моментов инерции Iz z , Ixx и Ixz , используя оценки матриц A и C , а также заранее известные значения баз стенда Lx , Lz , массу m и координаты центра масс rx и rz тела. Отметим, что данный алгоритм применим при любом способе возбуждения колебаний и пружинах различных неизвестных жесткостей.

41


3 Идентификация линейной колебательной системы
Существует множество методов идентификации линейных систем [15, 18, 19, 20, 21]. Обзоры методов идентификации показывают, что точность идентификации можно повысить используя априорную информацию. В нашем случае, мы не знаем функцию распределения помехи сигнала, следовательно методы максимального правдоподобия не подходят. Зато известно количество колебательных мод в сигнале. В настоящей работе рассмотрены два известных метода идентификации, обладающие высокой точностью и удобством. Первый классический метод Прони [26], в котором используется метод наименьших квадратов. Второй метод Кунга [13, 14] редукции системы при помощи сингулярного разложения, относящийся к активно развивающимся в последнее время методам пространства состояний [15, 19, 20]. Сравнения с другими методам [15] показывают, что они или мало отличаются или уступают в точности выбранным методам. Отметим, что методами идентификации нельзя непосредственно определить параметры системы (24), поскольку существует множество систем вида (24), имеющих одинаковые решения с различными матрицами M,

K и C, записанными в одной системе координат. Такие системы называются изоспектральными [12]. Поэтому будем идентифицировать систему следующего вида

q + Aq = 0 Е , y = Cq

q(0) = q0 ,
42

q(0) = v0 ,

(38)


где A = M-1 K. Результатом идентификации является определение оценок параметров A, C, q0 и v0 системы (38) по последовательности векторов y сигнала, измеренного датчиками. Оценки параметров обозначаются теми же символами с волной. Предполагается, что матрица A не имеет кратных собственных значений, так как в противном случае одной последовательности значений вектора y будут соответствовать множество систем (38) с разными матрицами наблюдения C в одном и том же базисе. Во втором параграфе этой главы были приведены способы определения оценок параметров системы (24) по оценкам параметров системы (38), используя различные дополнительные условия. Эти способы являются ключевым результатом данной работы.

3.1 Определение частот, декрементов затухания и амплитуд сигнала методом Прони
Одним из классических методов определения частот и декрементов затухания сигнала является метод Прони [26]. Метод Прони предполагает, что число мод колебаний в сигнале известно заранее. Определив частоты колебаний, найдем оценку матрицы системы (38) в нормальных координатах

~ A , так как эта матрица диагональна и е? диагональными элементами являются квадраты частот колебаний. Матрица наблюдения в нормальных

~ координатах C составляется из амплитуд мод колебаний, определенных ~ для каждого датчика. Количество строк матрицы C равно количеству
датчиков, а количество столбцов равно количеству мод колебаний систе-

43


мы. Тогда, если начальные скорости системы равны нулю v0 = 0, решение системы (38) запишется в виде

q = (cos(1 t), cos(2 t), cos(3 t))T . y = Cq

Полагая t = 0, получим вектор начального смещения системы в нормальных координатах, состоящий из единиц q0 = (1, . . . , 1)T . Итак, мы показали, что найдя частоты и амплитуды мод колебаний для каждого датчика, мы идентифицируем линейную консервативную динамическую систему (38), причем параметры системы получаем в нормальных координатах. Опишем суть метода Прони и укажем особенности его практического применения. Одномерный гармонический сигнал с затуханием представляется следующим образом:
P

x [n] = A0 +
k =1

Ak [exp (ik nt - k nt + ik ) + exp (-ik nt - k nt - ik )] 2 (39)

где P - число мод колебаний в сигнале, t - период дискретизации сигнала, Ak - амплитуда k -ой моды, k - декремент затухания k -ой моды,

k - частота k -ой моды, k - фаза k -ой моды. Для того, чтобы определить все частоты и декременты затухания, находятся коэффициенты ak разностного уравнения:
2P +1

ak x [n - k ] = 0,
k =0

(40)

где 2P + 1 n N - 1 , N - число элементов в массиве x. Чтобы найти
44


коэффициенты ak , нужно решить систему из 2P + 1 линейных действительных уравнений, положив a0 = 1. Эта система получается из условия:



N -1

2P +1

2



ak (k =0,2P +1)

min

ak x[n - k ]
k =0

.

n=2P +1

Далее находятся корни следующего полинома степени 2P + 1:
2P +1

F (z ) =
k =0

ak z

2P +1-k

= 0.

(41)

Один из его корней является вещественным, а остальные составляют комплексно-сопряженные пары:

zk = exp (ik t - k t) z
k +1

= exp (-ik t - k t)

(42)

Таким образом, из корней (42) полинома (41) мы находим частоты и декременты затухания колебаний. Для определения значений hk =
Ak 2

exp (ik )

нужно решить систему из 2P + 1 линейных уравнений с комплексными коэффициентами. При применении метода автором обнаружено, что он дает очень точные ( 0, 001%, N = 100) значения частот и декрементов затухания слабо зашумленного сигнала (белый шум 1%) при

d [2, 10) max
где
max

(43)
2 t

максимальная частота составляющих сигнала и d =



частота дискретизации. Точность повышается при приближении отношения частот к левому краю интервала в соотношении (43) при постоянном
45


e
5

(%)

еш гр по

ь ст но

оп ре де ле ни яа мп ли туд ы

ы от а ст яч ен и д ел п ре о ость еш н погр

0

2

N

wd/ wmax

Рис. 11: Зависимость погрешности от дискретизации сигнала

числе точек сигнала, что отражено на рисунке 11 сплошной линией. Из теоремы Котельникова [29] также следует, что левая граница интервала не может быть меньше 2. В то же время для определения амплитуд мод колебаний частоту разбиения сигнала по времени лучше выбирать так, чтобы на период моды
2 x

, амплитуду которой мы хотим определить, при-

ходилось бы примерно столько же точек, сколько периодов этой моды укладывается в суммарное время измерения N ћ t (см. рис. 11, пунктирная линия), то есть:

d max N .

(44)

Отсюда следует, что при больших N соотношение (43) может не выполняться. В таком случае для определения частот нужно рассматривать не все точки сигнала, а каждую m-ую

m

d . 2max
46

(45)


Тогда соотношение (43) будет выполнено, так как вместо d в (43) будет
фигурировать d = d /m.

Недостатком метода Прони является его чувствительность к шуму, что требует применения фильтров.

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

^ x = Ax + Bu,

^ y = Cx + e,

x(0) = x0 ,

(46)

где y(t) l-мерный вектор измерений, e(t) ошибки измерений (стационарный случайный процесс с нулевым средним), u внешнее (управляющее) воздействие. Задача идентификации состоит в определении матриц

^ ^ A, B и C по известной реакции системы y(t) на некоторое входное воздействие u(t). Уравнения движения механической системы можно записать в виде

(46), где x(t) = (q(t), q(t)) 2n-мерный вектор обобщенных координат ^ ^ и скоростей, n число степеней свободы системы, а матрицы A, B, C
имеют блочную структуру





^ C= CC
v

0 B
f

^ A=

0

E
d

,

,

B=

.

AA

(47)

Современные методы идентификации методами пространства состояний [13, 14, 21] дают представление системы в некотором абстрактном пространстве, в котором блочная структура матриц (47), вообще говоря, не
47


сохраняется. Пространство этого представления связано с исходным фазовым пространством линейным невырожденным преобразованием, отыскание которого часто является нетривиальной в вычислительном плане задачей. Однако для консервативных механических систем (38), где Ad =

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

^ Sj = CAj -1 B,

j = 1, 2, . . . , m,

(48)

где m количество марковских параметров, A переходная матрица непрерывной по времени системы (46). Параметры Sj можно найти, измерив импульсную реакцию

^ ^ h(t) = C exp(At)B

(49)

непрерывной системы (46). Подставив (47) в (48) и положив Ad = 0, получим, что последовательность марковских параметров Sj распадается на две подпоследовательности:

S

2k -1

^^ = CA ^^ = CA

2k -2 2k -1

B= B=

Cv A CA

k -1

Bf , k = 1, 2, . . . , N ,

(50) (51)

S2
где N =

k

k -1

Bf ,

m 2

. Таким образом, применение алгоритма идентификации

непосредственно к последовательностям (50) и (51) для нахождения матриц A, C, Bf обеспечивает необходимую блочную структуру системе (46).
48


Запишем систему (24) в виде (46), положив x = (q, q), A = -M-1 K.
Будем считать, что внешнее воздействие отсутствует, т. е. u = 0, но зада-

ны начальные условия q(0) = q0 , q(0) = v0 . Тогда решение системы (46)
запишется в виде

^ x = exp(At)x0 ,

^ ^ y = C exp(At)x0 ,

Вектор y(t) имеет такой же вид, что импульсная реакция (49) системы (46). Значит для идентификации системы (24) можно использовать алгоритм, представленный для системы (46), с точностью до замены матрицы

B на вектор x0 . Соотношения аналогичные (50), (51), при Cv = 0, запишутся следующим образом:

S2

k -1 k

^^ = CA2 ^^ = CA2

k -2

x(0) = x(0) =

CA CA

k -1 k -1

q0 , v0 , k = 1, 2, . . . , N .

(52) (53)

S2

k -1

Найдем марковские параметры Sj системы (46). Импульсная реакция (49), как уже было сказано, при отсутствии управления совпадает с решением системы и может быть записана следующим образом

^ ^ h(t ћ i) = yi = C exp(A ћ t ћ i)x0 . ^ ^ Методом пространства состояний определим матрицы C, exp(A ћ t) и ^ вектор x0 в некотором базисе [15]. Чтобы найти матрицу A в этом пространстве, находим преобразование T в пространство собственных векто-

^ ров матрицы exp(A ћ t), где она диагональна D=T
-1

^ exp(A ћ t)T.
49

(54)


^ Вычислив логарифмы собственных значений матрицы exp(A ћ t), поделим их на t и преобразуем полученную диагональную матрицу в исходное пространство с помощью обратного преобразования T
-1

^ A = T(ln(D)/t)T-1 .
Теперь марковские параметры Si определим по формуле (48). Итак, известна последовательность {Sj }
2N +1 j =1

и матрица наблюдения

C. Требуется определить матрицу A, матрицу наблюдения C и векторы
начального состояния q0 , v0 . Рассмотрим только подпоследовательность (52). С помощью алгоритма идентификации Кунга [13] по известной последовательности {S2
N k +1 }k =1

~~~ найдем представление системы {A, C, q0 } в

некотором пространстве состояний. Алгоритм Кунга применительно к рассматриваемой задаче состоит в следующем. Составим матрицу Ганкеля H из векторов столбцов {S2k+1 }N=1 k



S3 H= ... S
где r =
N 2 2r+1

S1

S S

3 5

. . . S2 . . . S2 ...


r+1

... S
2r+3

r+3 , ...
r+1

. . . S4

. Пусть U и V - ортогональные матрицы, осуществляющие

сингулярное разложение матрицы Ганкеля H:

H=U

0 0

VT ,

где = diag{1 , 2 , . . . , n }, n число степеней свободы системы, а
50

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

U=

U

1

U

2

, V=

V1 V2

,

(55)

где dim(U1 ) = (r + 1)l Ч n, dim(V1 ) = (r + 1) Ч n. Построим матрицы

~ ~ наблюдаемости (r + 1)l Ч n и управляемости n Ч (r + 1)
1 ~ = U1 2 , 1 T ~ = 2 V1 .

(56)

~ Пусть 1 и 2 матрицы, полученные из вычеркиванием соответственно последней и первой блок-строк l Ч n. Тогда

~ A = T ћ 1 1

-1

T ћ 2 . 1

(57)

~ ~ ~ Матрица C вычисляется как первая блок-строка матрицы , а q0 как ~ первый столбец матрицы .
Таким образом, идентификация проводится дважды: первый раз для

^ определения матрицы exp(A ћ t), второй раз для определения матриц C
и A. Теперь можно выразить по ним матрицу инерции системы способами, указанными в начале главы. Сингулярные числа 1 , 2 , . . . , n являются мерой оцениваемости. Они обратнопропорциональны погрешности определения соответствующих величин.

Идентификация рассматриваемой системы
Рассмотрим преимущества и недостатки представленных методов идентификации и разработаем алгоритм определения оценок параметров A, C,
51


Погрешность определения частоты (%)

6 5 4 3 2 1 0 0 20 40 60 80 100 120 Уровень шума относительно амплит уд мод к олебаний (%) Кунг Прони

Рис. 12: Сравнение погрешности определения частот методами Прони и Кунга

q0 и v0 системы (38). Метод Кунга удобнее метода Прони. Метод Кунга
существенно менее чувствителен к шуму (см. рис. 12), то есть не требует дополнительной фильтрации сигнала, так как модель системы в методе Кунга учитывает шум. В отличие от метода Прони метод Кунга работает с многомерным сигналом. Следовательно, он нормально функционирует, даже если одна из мод колебаний на одном датчике не наблюдается. Для определения частот и декрементов затухания сигнала будем использовать метод Кунга. Частоты и декременты затухания сигнала определя-

^ ются находятся из диагональной матрицы D, см. (54). Матрица A имеет ^ ^ размерность 7, а матрица C 4 Ч 7. Первый столбец матрицы C отражает
среднее значение измеряемого сигнала. Первым диагональным элементом

52


матрицы D является единица, остальные элементы равны:

D D

2k 2k

=e =e

(-k +ik )t (-k -ik )t

, ,

k = 1, 2, 3.

2k +1 2k +1

Оценка матрицы системы (38) запишется следующим образом

ln(D2k 2k ) ~ Ak k = t

2

k = 1, 2, 3.

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

^ yj = C ћ T ћ D

j -1

T- 1 x 0 ,

j = 1, ..., N .

Вектор T-1 x0 можно представить как произведение диагональной матрицы diag(T-1 x0 ) на вектор i = (1, . . . , 1)T . Учитывая, что матрица D также диагональна, а произведение диагональных матриц коммутативно, преобразуем уравнение сигнала

^ yj = C ћ T ћ diag(T-1 x0 ) ћ D

j -1

i,

j = 1, ..., N .

^ Отсюда получаем матрицу Q = C ћ T ћ diag(T-1 x0 ) размерности 4 Ч 7,
первый действительный столбец которой содержит значения положения равновесия на четырех датчиках, а последующие столбцы комплексносопряженные векторы форм колебаний. Матрицу наблюдения составим следующим образом

~ Cj k = sign(Re(Qj 2k ))|Qj 2k |,

j = 1, 2, 3, 4,

k = 1, 2, 3.

Такая запись матрицы наблюдения означает, что мы положили вектор начального смещения q0 = (1 1 1)T , и вектор начальной скорости v0 =

(0 0 0)T .
53

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

54

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

Mq + Kq = CT u Е , y = Cq |uj | Uj ,

q(0) = q0 , j = 1, . . . , 4 ,

q(0) = v0 ,

(58) (59)

где u вектор управления размерности 4 содержит значения создаваемых силовозбудителями сил, CT матрица управления размерности 3 Ч 4 совпадает с транспонированной матрицей наблюдения, так как управляющие силы прилагаются в точках расположения датчиков. Заметим, что в силу последнего обстоятельства система (58) наблюдаема и управляема одновременно [35]. В случае, когда для возбуждения колебаний системы (58) используются алгоритмы управления, возникает задача минимизации времени возбуждения. Для многомерной системы эта задача очень сложна. Автору не удалось найти в литературе решение в общем случае. Для оценки эффективности (в смысле быстродействия) возможных алгоритмов управ55

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

x1 + 2 x1 = m-1 u , Е

|u | u , m

(60)

где x1 скалярная переменная, m масса тела, частота свободных колебаний. Требуется найти функцию управления u , приводящую систему (60) из любого начального положения x1 (0), x1 (0) к уровню энергии с амплитудой колебаний a за кратчайшее время. Запишем уравнение системы (60) в форме Коши









x=
где вектор x = (x1 ,

0



x +

0 u

,

- 0
-1

x(0) = x0 ,

(61)

x1 )T и функция управления u = m-1 -2 u . Усло

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

(x) = a2 - x2 - x2 = 0. 1 2

(62)

Решим задачу отыскания оптимального управления, используя принцип максимума Понтрягина [24, 35]. Для этого рассмотрим систему, сопряженную с системой (61)

0 - 0
56

.
(63)

=


Эта система имеет решение



1 2

= C sin( t + ), = C cos( t + ).
(64)

Гамильтониан задачи имеет следующий вид

H = 0 + u2 + x2 1 - x1 2 0,

(65)

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

u дается принципом максимума Понтрягина max H
u

(66)

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

i (T ) = c

xi

t=T

,

(67)

где T конечный момент времени, а c некоторая константа. В данной задаче условие максимума (66) записывается следующим образом

u = um sign(2 ),
где

(68)

um = m-1 -2 u . m
Из условий трансверсальности (67) следуют выражения

(69)



1 2

= cx1 , = cx2 ,
57

(70)


связывающие координаты системы (61) с координатами сопряженной системы (63) в конечный момент времени T . Теперь можно выразить функцию оптимального управления через фазовые координаты, используя условия максимума (68) и трансверсальности (70) вместе с решением (64) сопряженной системы (63). Это означает, что мы построим синтез оптимального управления (см. следующий параграф).

58


4.1 Синтез управления
Функция синтезированного управления u не зависит явно от времени, а только от фазовых переменных системы (координаты и скорости). Следовательно, такая функция управления при любых начальных условиях приводит систему в требуемое положение вне зависимости от возмущений. Поэтому на практике стремятся построить синтез управления. Функция оптимального управления принимает значения um или -um ввиду условия (68). Таким образом, оптимальная траектория системы (61) на фазовой плоскости состоит из дуг концентрических окружностей с центрами в точках (um , 0) или (-um , 0). Множество точек, в которых траектория переходит из дуги окружности с одним центром в дугу окружности с другим центром, называется линией переключения. Синтезировать оптимальное управление означает найти линию переключения и определить знак функции управления u с обоих сторон этой линии. Случай, когда a = 0 классический пример задачи быстродействия, рассмотренный D.W. Bushaw [34, 35]. Линия переключения оптимальной траектории при a = 0 показана на рисунке 13. Также известна задача выхода на фазовую окружность a > 0 снаружи [36]. Однако автору неизвестны публикации, где решается задача быстродействия для выхода на фазовую терминальную окружность изнутри. Построим линии переключения задач быстродействия изнутри и снаружи терминальной фазовой окружности различных радиусов. Для того, чтобы построить линию переключения при 0 < a < u
59
m

(см. рис.


Рис. 13: Линия переключения и пример оптимальной траектории при a = 0.

14) заметим из решения (64) сопряженной системы (63), что между двумя переключениями изображающая точка должна пройти угол вокруг центра одной окружности. Таким образом, линия переключения симметрична относительно начала координат. Рассматривая процесс в обратном времени, когда движение начинается от терминальной поверхности (62), мы можем определить из условий трансверсальности (70) и решения сопряженной задачи (64), где должно произойти первое (в обратном времени) переключение в зависимости от начальной точки на терминальной поверхности. Найдем угол , по прошествии которого происходит первое (в обратном времени) переключение управления. Из решения (64) сопряженной системы (63) получим

sin( ) =

2 (T )
2 2 1 (T ) + 2 (T )

.

(71)

60


Рис. 14: Линия переключения при 0 < a < um .

x2

a x1 -um +um

Рис. 15: Линия переключения при a = um .

61


x2

-um +um

a

x1

Рис. 16: Линия переключения при a > um .

62


Используя условия трансверсальности (70), можно представить это выражение через x1 и x
2

sin( ) =

c x 2 (T ) |c| x2 (T ) + x2 (T ) 1 2

.

(72)

Приняв во внимание уравнения терминальной поверхности (62), находим, что

sin( ) = + x2 (T )/a.
Таким образом,

(73)

= + arcsin(x2 (T )/a).

(74)

Например, легко видеть, что если стартовать в обратном времени из точки

(0, -a) управление переключится по прошествии угла /2 по окружности
вокруг точки (um , 0). Аналогично строятся линии переключения при a >

um (см. рис. 16).

63

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

um входит частота. Поэтому невозможно построить линию переключений
так, как это делалось ранее. Однако, если известно, что a

um , то за линию переключений можно

принять ось абсцисс (см. рис. 16). Функция управления

um sign(x2 ) при x2 + x2 < a2 1 2 u= 0 при x2 + x2 = a2 1 2 -u sign(x ) при x2 + x2 > a2 m 2 1 2
в этом случае примет вид

u = -um sign ( (x)) sign

x2

.

(75)

Будем называть это управление квазиоптимальным.

4.3 Максимальное время возбуждения системы
Теперь мы можем оценить максимально возможное время T возбуждения системы оптимальной (в смысле быстродействия) малой (a

um ) управ-

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


Рис. 17: Линии переключения при a > um .

возбуждения колебаний будет, если x1 (0) = x1 (0) = 0. Легко видеть из рисунка 16, что за один период амплитуда колебаний системы изменяется на величину 4 um . Таким образом, деля a на 4 u
m

и умножая на 2 / ,

мы получаем время на возбуждение колебаний до амплитуды a

T=
После подстановки um = m-1

a . 2 um
в итоге получаем (76)

-2 um

T=

am . 2 u m

Определим как время возбуждения колебаний (76) отличается от времени при оптимальном управлении. На рисунке 17 показаны различия траекторий на фазовой плоскости при оптимальном управлении и при
65


квазиоптимальном управлении (75). Время пропорционально сумме углов, которые проходит изображающая точка вокруг центров концентрических окружностей, по дугам которых она движется. Сравним эти суммы углов для оптимального и квазиоптимального управлений. Без потери общности рассмотрим процесс с начальной точкой A на оси абсцисс под нечетной полуокружностью линии переключений, считая от терминальной окружности. При оптимальном управлении изображающая точка проходит угол вокруг точки (um , 0), затем по n оборотов на попеременно вокруг -(um , 0) и (um , 0) и наконец угол от C до точки D на терминальной поверхности.

+ 2 n + - оптимальное управление, 2 n + - квазиоптимальное управление,
где n количество переходов через линию переключения слева от начала координат после точек B и A. Так как сумма углов при оптимальном управлении наименьшая, можем записать следующее неравенство



+

(77)

Кроме того, на рисунке 17 видно, что на терминальной окружности точка

D находится правее точки F , следовательно справедливо неравенство - .

Действительно, в силу симметрии точки E и C располагаются на пунктирной полуокружности так же как точки A и B . Дуга AB образована вокруг центра (um , 0) из-за чего расстояние от B до (-um , 0) больше чем
66


от A до (-um , 0). Аналогично, расстояние от C до (-um , 0) больше чем от E до (-um , 0) и дуга C D имеет больший радиус чем дуга E F , из-за чего на терминальной окружности точка D находится правее точки F . Добавив к левой части неравенства неотрицательный угол и учитывая (77) получаем


При a

+

- .

um угол um /a, следовательно, можно записать = + + O(um /a).

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

T = T0 +

O(um /a) ,

где T0 время оптимального процесса. Значит, время (76) можно считать оптимальным с точностью до величины

O(um /a) .
Такая же оценка получается, если начальная точка лежит внутри терминальной окружности, но не близко к началу координат. Если начальная точка находится близко к началу координат, разница квазиоптимального и оптимального времени не превышает / . По сравнению с возбуждением изменяющейся по синусу силой такой же амплитуды u m

x + 2 x = u sin( t)/m, Е m

x(0) = x(0) = 0

(78)

возбуждение колебаний силой (75) происходит в 4/ раз быстрее.
67


4.4 Оценка минимального времени возбуждения многомерной колебательной системы
Многомерная линейная консервативная система описывается следующим уравнением

Е Mq + Kq = B u, |uj | Uj ,

q(0) = q(0) = 0,

(79) (80)

j = 1, . . . , m ,

где матрицы M и K симметрические и положительно определенные,

q вектор переменных размерности n, Uj ограничение амплитуды j ой компоненты вектора управления u, имеющего размерность m, B матрица управления размерности n Ч m. Задача состоит в определении минимального времени, за которое все моды колебаний системы достигнут соответствующих уровней энергии. Известно, что существует линейная взаимнооднозначная замена координат, обозначим е? матрицу S, которая приводит матрицу M к единичной матрице и матрицу K к диагональной матрице 2 , содержащей квад2 раты собственных частот i системы

Е z + 2 z = Bu,

z(0) = z(0) = 0,

(81)

где z вектор нормальных координат системы и матрица B = ST B . В нормальных координатах система распадается на n одномерных осцилляторов, связанных только через управление. Каждое из уравнений
2 zi + i zi = (Bu)i , Е

zi (0) = zi (0) = 0,

(82)

совпадает с уравнением (60) с точностью до обозначения.
68


Рассмотрим процесс возбуждения отдельно каждого i-го осциллятора при ограничении на функцию управления
m

|(Bu)i | bi =
j =1

(|Bij |Uj ) .

(83)

Используя оценку (76) времени возбуждения одномерного осциллятора (60) оптимальной (в смысле быстродействия) управляющей силой, можно выбрать моду многомерного осциллятора с максимальным временем возбуждения

T

min

= max ai 2i

m i j =1

(|Bij |Uj ) .

(84)

Эта формула получена при условии, что задачей является возбуждение только выбранной i-ой моды. Многомерный осциллятор нельзя возбудить быстрее, чем его одну моду с самым большим минимальным временем возбуждения, при тех же ограничениях на управление. Действительно, предположим, что T0 время возбуждения многомерного осциллятора при помощи управления u0 (t), T минимальное время возбуждения самой медленной моды и T0 < T , тогда при использовании управления u0 (t) для возбуждения этой моды время возбуждения е? колебаний будет T0 , что меньше минимального времени T . Мы пришли к противоречию. Мы получили оценку минимально возможного времени возбуждения многомерного осциллятора ограниченными силовыми воздействиями. Эту оценку можно использовать при тестировании возможных алгоритмов управления, чтобы определить их эффективность в смысле скорости возбуждения колебаний.
69


4.5 Пример расчета
Расчитаем минимальное время возбуждения многомерной колебательной системы измерительного стенда, рассматриваемого в данной работе. Рассмотрим самый простой случай, когда все четыре пружины имеют одинаковые продольную k = 7.4 ћ 105 Н/м и поперечную c = 3.1 ћ 105 Н/м жесткости, центр масс находится посередине между пружинами, то есть

rx = Lx /2, rz = Lz /2, и главные оси инерции тела направлены вдоль
осей симметрии стенда. Тогда матрицы инерции и жесткости будут диагональны, и крутильные колебания тела будут происходить вокруг осей симметрии стенда. Масса тела

m = 15000кг,
главные моменты инерции

Iz z = 44000кг ћ м2 ,
базы стенда

I

xx

= 3000кг ћ м2 ,

Lx = 4м,

Lz = 1.74м.

Пусть необходимо возбудить колебания системы (24) из положения q(0) =

q(0) = 0 в состояние, в котором амплитуды мод колебаний ai = max qi (t),
t

i = 1, 2, 3

будут равны

a1 = am ,

a2 =

am , Lx

a3 =

am , Lz

(85)

70


где am = 2 ћ 10-3 м. Колебания возбуждаются устройствами, создающими силовые воздействия в точках крепления пружин в вертикальном направлении, имеющие ограничения по значению силы Um = 5Н. После подстановки значений положения центра масс в (23) матрица C примет вид

1 -Lx /2 Lz /2 1 Lx /2 Lz /2 C= 1 Lx /2 -Lz /2 1 -Lx /2 -Lz /2

.

Тогда матрица жесткости согласно (25) будет равна

. k

4k 0 0 K = k CT C = 0 L2 k 0 x 0 0 L2 z

Матрица управления равна транспонированной матрице C

B = C T ,
так как управляющие силовые воздействия прикладываются в точках крепления пружин к телу. Тогда матрица B примет вид

B=M B =-
-1

1 m Lx 2Izz

1 m Lx 2Izz Lz 2Ixx

1 m Lx 2Izz

1 m

.

- -

Lx 2I z z Lz 2Ixx

Lz 2Ixx

-

Lz 2Ixx

Частоты системы вычисляются по формуле

i =

Cii /Mii ,
71

i = 1, 2, 3,


где M

ii

и Cii диагональные элементы матриц инерции и жесткости.

Отсюда имеем выражения для частот мод колебаний

1 = 2

k /m,

2 = Lx

k /Izz ,

3 = L

z

k /Ixx .

(86)

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

Ti =
где bi = U получаем
m m j =1

ai i , 2 bi

i = 1, 2, 3,

|Bij |. Отсюда после подстановки выражений (85) и (86)


T T T

1 2 3

= = =

am 4 Um am 4 Um am 4 Um

mk =
=

33с , 14с ,

Iz z k Lx Ixx k Lz

= 8.5с .

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

72


Глава III. Численный эксперимент и анализ погрешностей 1 Численное моделирование процесса колебаний
Для исследования процесса колебаний системы и проверки алгоритмов определения моментов инерции написана программа [30], решающая численно уравнения динамики системы с менее строгими предположениями, чем в представленных во второй главе алгоритмах. Во-первых, тело может совершать движения по всем шести степеням свободы, то есть возможны малые поперечные колебания. Во-вторых, не делается предположения о малости колебаний, то есть учитываются геометрические нелинейности. Движение системы получается путем численного решения уравнений динамики

dK =M dt d2 x m 2 = F, dt

(87) (88)

где M результирующий момент внешних сил, F результирующий вектор внешних сил, m масса тела, K кинетический момент. Уравнение (88) записано в неподвижной системе координат Oa X Y Z (см. рис. 8). Уравнение (87) записано в Кениговой системе координат [25], то есть системе координат с началом в центре масс, движущейся поступательно относительно неподвижной системы отсчета. Кенигова система коорди-

73


нат удобна тем, что уравнения динамики в ней имеют тот же вид, что и в инерциальной системе отсчета. Пусть оси Кениговой системы координат совпадает с осями системы Oa X Y Z . Кенигова система координат обозначена на рисунке 8 C X Y Z . Так как силы, за исключением силы тяжести, являются силами упругости, то движение будет колебательным. Уравнения (87) и (88) можно переписать в форме Коши



dK dt d dt dv dt dx dt

=M = W S I-1 ST K = F/m =v ,
(89)

где v вектор скорости центра масс в неподвижной системе отсчета

Oa X Y Z , S матрица преобразования поворота системы C xy z вокруг
осей Кениговой системы координат C X Y Z на углы x , y и z , матрица W преобразует вектор угловой скорости в вектор производных углов Крылова [51]







1 x y = 0 0 z

sin(x ) sin(y ) cos(y )

-

cos(x ) sin(y ) cos(y )





cos(x ) -
sin(x ) cos(y )

sin(x )
cos(x ) cos(y )

x y . z

Система, описываемая уравнениями (89), имеет 6 степеней свободы, но, приложив силу параллельно оси Oa Y , мы можем возбудить колебания преимущественно по трем степеням свободы (рис. 8). Перед выполнением расчета в программу задаются все параметры колебательной системы: масса, тензор инерции, положение центра масс, положение и жесткости
74


Вертикальное смещение (мм)

2.0 1.5 1.0 0.5 0.0 0.0 -0.5 -1.0 -1.5 -2.0 -2.5 0.2 0.4 0.6 0.8 1.0 1.2

Время t (c)

Рис. 18: Вертикальные смещения, измеряемые четырьмя датчиками.

пружин, на которых установлено тело, начальное силовое воздействие, а также начальные условия. Результатом работы программы является определение собственных частот и амплитуд мод колебаний системы. Эти величины определяются при помощи метода идентификации описанного в третьем параграфе второй главы. Процесс моделирования колебаний представляет собой численное решение системы нелинейных дифференциальных уравнений (89) [27, 28]. Ниже приведена последовательность операций для определения векторов F и M в правых частях уравнений (89) на одном промежутке времени 1. Рассчитывается новое положение подвижной системы отсчета по прошествии времени 2. Находится матрица поворота S системы отсчета, связанной с телом

C xy z , относительно е? прежнего положения.
75


3. Определяются новые радиус-векторы точек крепления пружин к телу в Кениговой системе отсчета

ri = S ћ ri ,

i = 1, 2, 3, 4.

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

li = x + ri - Li ,

i = 1, 2, 3, 4,

где Li координаты закрепленных на основании концов пружин. 5. Отсюда находятся векторы сил, действующих на тело со стороны пружин

fi = -ki li 1 -

(l0i ћ li ) |li |2

- ci

(l0i ћ li ) li - l0 |li |2

i

,

i = 1, 2, 3, 4,

где ki и ci продольная и поперечная жесткости i-ой пружины (4),

l0i векторы li при недеформированных пружинах.
6. Определяется суммарный вектор внешних сил
4

F=
i=1

fi + m g,

где g - вектор ускорения свободного падения. 7. Определяется суммарный вектор моментов внешних сил
4

M=
i=1

ri Ч fi .

76

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

77


2 Пример вычислений динамических параметров тела
Для расчета задаются параметры системы (24), указанные в Таблицах 1 и 2. Положение центра масс отсчитывается от точки крепления первой пружины к телу (см. рис. 9). Начальное смещение системы от положения равновесия вызывается силой 2000н, приложенной в вертикальном направлении в точке крепления второй пружины к телу. В результате работы программы получается 12544 векторов наблюдения y, которые содержат значения вертикальных смещений четырех точек тела, к которым крепятся пружины. Период, с которым фиксируются измерения, составляет 0.002212 с. По полученному сигналу определяются частоты колебаний системы (см. таб. 3) и матрица, составленная из амплитуд мод колебаний на четырех датчиках. Эта матрица является матрицей наблюдения системы в нормальных координатах

-.000675 .000673 -.000673 -.000675 -.000719 -.000630 C = -.000675 -.000673 .000673 -.000675 .000719 .000630

,

причем вектор начального смещения системы от положения равновесия в нормальных координатах состоит из единиц q0 = (1, . . . , 1)T , что показано во второй главе. Матрица системы в нормальных координатах A диагональна. Диагональными элементами матрицы A являются квадраты ча78


Жесткость каждой из четырех пружин База стенда по оси Ox База стенда по оси Oz

740000 Н/м

Lx Lz

4м 1.74 м

Таблица 1: Параметры стенда.

Масса тела Положение центра масс по оси Ox Положение центра масс по оси Oy Положение центра масс по оси Oz

m r
x

15000 кг 2м 0м 0.87 м

ry rz

-Iy -Iz Ixx
x x

-Ixy Iy
y y

-Iz

Tензор инерции. -Ixz 3000 0 -400 -Iyz = 0 43000 0 Izz -400 0 44000



Таблица 2: Заданные параметры тела.

Частота вертикальных колебаний Частота колебаний вокруг оси близкой к Ca Z Частота колебаний вокруг оси близкой к Ca X

1 2 3

14.04726 рад/с 16.39794 рад/с 27.35160 рад/с

Таблица 3: Частоты колебаний системы.

79


Масса тела Положение центра масс по оси Ox Положение центра масс по оси Oz Момент инерции относительно оси C z Момент инерции относительно оси C x Центробежный момент инерции

m rx rz I
zz

15007 кг 2м 0.87 м 44031 кгћм 3006 кгћм 399.9 кгћм
2 2 2

Ixx Ixz

Таблица 4: Найденные динамические параметры тела

Рис. 19: Погрешность определения центральных моментов инерции из-за наклона главной оси инерции вверх от оси C x.

стот колебаний. По результатам идентификации системы и по известному начальному силовому воздействию определяются значения динамических параметров тела, представленные в таблице 4. Видно, что погрешности определения динамических параметров тела находятся в пределах 0.2%. Этот результат получен при добавленном в сигнал белом шуме в 1% от средней амплитуды мод сигнала на датчиках. Многочисленные численные эксперименты показали, что представлен80


ные алгоритмы определения моментов инерции удовлетворяют требуемой точности 2%, если соответствующая главная ось инерции отклонена вверх от оси C x не более чем на 3 градуса (см. рис. 19). Это является подтверждением сделанного в первой главе предположения о возможности описания движения системы тремя степенями свободы вместо шести. Для определения моментов инерции Iyy и Ixy измерения повторяются с телом, повернутым на 90o вокруг оси C x, при этом ось C y переходит в

C z , а момент инерции Ixx определяется повторно. В качестве окончательного значения Ixx рекомендуется брать значение при измерении с максимальным центробежным моментом Ixz , так как при этом отрицательно влияющий на точность момент I
xy

(момент I

xz

при другом измерении) ми-

нимален. В результате двух измерений определяются 5 из 6-ти моментов инерции тела (кроме Iyz ).

81


3 Анализ погрешности определения динамических параметров
3.1 Анализ чувствительности динамических параметров к ошибке идентификации
Выражения для динамических параметров тела через результаты идентификации системы во второй главе получены в явном виде. Это позволяет получить аналитически чувствительности динамических параметров к ошибке идентификации. При расчете чувствительности используются точечные оценки результатов идентификации. Предполагается, что их относительные погрешности могут быть ограничены по модулю. В результате идентификации системы мы получаем оценку матрицы наблюдения в нормальных координатах C и оценки частот колебаний системы 1 , 2 и 3 . Рассмотрим первый алгоритм определения моментов инерции. Составим вектор-функцию из выражений для динамических параметров тела (см. приложение 1)

= (M , rx , rz , Izz , Ixx , Ixz )T .
Вектор-аргумент состоит из 18 элементов

= (p, Lx , Lz , 1 , 2 , 3 , C11 , C12 , . . . , C43 )T ,
где p значение приложенной в вертикальном направлении до начала движения силы; Lx и Lz базы стенда, а также расстояния между датчиками по осям Ox и Oz . Пусть погрешность j -го элемента вектор-аргумента
82

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

i =
j =1

i ( ) j , j

i = 1, . . . , 6.

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



i i i
18 j =1

i ( ) j , j

i = 1, . . . , 6.

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

i 2 ћ , i

i = 1, . . . , 6.

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

3.2 Анализ чувствительности динамических параметров к ошибке измерения сигнала
Покажем, как изменяется относительная погрешность определения динамических параметров при добавлении в сигнал белого шума. Уровень белого шума задается в процентах от значения 0.00067 м амплитуд мод колебаний на датчиках. Численные эксперименты показывают (см. рис.
83


14
Погрешность Izz (%)

12 10 8 6 4 2 0 0 2 4 6 8 10 12
Уровень шума относительно амплитуд мод колебаний (%)

Рис. 20: Чувствительность динамических параметров к уровню шума в сигнале.

20), что для определения динамических параметров с погрешностью не более 2% требуется, чтобы погрешность измерения сигнала не превышала

6%. Относительная погрешность вычислялась по следующей формуле Iz z - Iz Iz z
z

100%,

где Iz z значение момента инерции, задаваемое при численном моделировании системы; I
zz

оценка, полученная в результате применения алго-

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

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


Погрешность

0.3 0.2 0.1 0 -0.1 -0.2 -0.3

Смещение центра масс относительно центра жесткости вдоль оси Z (м)

9% -10% 8% - 9% 7% - 8% 6% - 7% 5% - 6% 4% - 5% 3% - 4% 2% - 3% 1% - 2% 0% - 1%

-0.4

-0.3

-0.2

-0.1

0.1

0.2

0.3
85

Смещение центра масс относительно центра жесткости вдоль оси X (м)

Рис. 21: Относительная погрешность определения динамических параметров в зависимости от положения центра масс относительно центра жесткости.

0.4

0


Рис. 22: Относительная погрешность определения динамических параметров в зависимости от положения центра масс относительно центра жесткости.

86


Рис. 23: Минимальная энергия мод колебаний в зависимости от положения центра масс относительно центра жесткости.

87


Рис. 24: Минимальная энергия мод колебаний в зависимости от положения центра масс относительно центра жесткости.

88


системы это точка в которой приложенная сила вызывает только поступательное смещение тела, то есть тело не наклоняется. Так как пружины одинаковы, то центр жесткости системы находится в центре симметрии стенда. Видно, что форма области низкой погрешности совпадает с формой высокой минимальной энергии мод колебаний (см. рис. 23 и 24), что естественно. Образно говоря, если мода колебаний имеет низкую энергию, то она утонет в шуме. Графики на рисунках 23 и 24 строились по численным расчетам движения системы, которая приводилась в движение устранением известной дополнительной силы, приложенной до начала движения в точке крепления второй пружины к телу в вертикальном направлении. Можно заключить, что условием применимости такого способа возбуждения колебаний является нахождение центра масс тела внутри области, где минимальная энергия мод колебаний больше, чем половина средней энергии мод колебаний. На рисунке 24 этот уровень энергии равен 0.08Дж. Если минимальная энергия ниже этого уровня, то следует подкачать соответствующую слабую моду колебаний. При этом для определения моментов инерции тела надо использовать алгоритмы, работающие без информации о возбуждении, представленные во второй главе. Какой из предложенных алгоритмов определения моментов инерции окажется удобнее, покажут практические испытания измерительного стенда.

89


4 Влияние демпфирования на частоты и формы колебаний системы
Рассмотрим влияние демпфирования на частоты и формы колебаний системы, а следовательно, и на алгоритмы определения моментов инерции. Рассмотрим два варианта линейного демпфирования. Первый, когда матрица демпфирования принимает диагональный вид в том же базисе, что и матрицы инерции и жесткости. Второй вариант, когда демпфирование считается малым. Для алгоритмов определения моментов инерции, представленных во второй главе, важно определить частоты системы (24) и матрицу S из (31) преобразования из нормальных координат системы (24) в исходные физические координаты. Таким образом, задача состоит в определении частот и векторов форм колебаний системы (24) по собственным числам и векторам системы (24) с демпфированием, которая описывается уравнением

Е Mq + Dq + Kq = 0,

(90)

где D симметрическая положительно определенная матрица демпфирования. Пусть xарактеристическое уравнение

det(M2 + D + K) = 0

(91)

имеет только комплексные некратные корни, то есть движение колебательное и происходит по n модам колебаний, 1 , . . . , n собственные числа системы (90), u1 , . . . , un соответствующие им собственные векторы.
90


Из собственных векторов системы (24) можно составить невырожденную матрицу U0 = S-1 . Известно [23, 38, 39, 40, 41, 42], что в базисе из собственных векторов системы (24) матрицы квадратичных форм M и K имеют диагональный вид. Если в этом же базисе матрица демпфирования имеет тоже диагональный вид, то системы (24) и (90) имеют одинаковые (с точностью до множителя) собственные векторы. Покажем это и найдем выражение для определения частот системы (24) через собственные значения системы (90). Матрица U0 невырожденная, следовательно det(UT ) = det(U0 ) = 0. 0 Умножим уравнение (91) слева на det(UT ), а справа на det(U0 ). 0

det(UT ) det(M2 + D + K) det(U0 ) = 0. 0

(92)

Так как произведение детерминантов матриц равно детерминанту произведения матриц, характеристическое уравнение (92) можно записать в виде

det(UT MU0 2 + UT DU0 + UT KU0 ) = 0. 0 0 0

(93)

Пусть собственные векторы нормированы так, что UT MU0 = E, где E 0 единичная матрица, тогда

det(E2 + 2 + 2 ) = 0,

(94)

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

91


в виде

n 2 (2 + 2j + j ) = 0, j =1

(95)

где n число мод колебаний, j декремент затухания и j частота колебаний j -ой моды колебаний. Из выражения

j

,n+j

= -j + i

2 j -

2 j

(96)

для собственных значений системы (90) получим выражение для частот системы (24)

j =

(Im(j ))2 + (Re(j ))2 .

(97)

Покажем, что собственные векторы системы (90) являются собственными векторами системы (24). Запишем выражение для нахождения собственных векторов системы (90) в базисе собственных векторов системы (24)

(E2 + 2 + 2 )u = 0.

(98)

Собственными векторами u в этой системе координат являются базисные векторы, так как матрица в скобках диагональная. А базис состоит из собственных векторов системы (24), что и требовалось показать. Следовательно, для определения моментов инерции системы (90) с диагональной матрицей демпфирования в нормальных координатах системы (24) матрица преобразования S = U
-1

из (31) определяется как и в случае без

демпфирования, а частоты соответствующей системы без демпфирования находятся по формуле (97). Численный эксперимент при демпфировании в пружинах 10000Нс/м показал, что такой тип демпфирования не влияет на точность определения моментов инерции. В данном случае матрица
92


Рис. 25: Вертикальные смещения, измеряемые четырьмя датчиками в случае с демпфированием.

демпфирования пропорциональна матрице жесткости. Затухание колебаний видно невооруженным глазом (см. рис. 25). Рассмотрим случай, когда матрица демпфирования D системы (90) в базисе собственных векторов системы (24) произвольная симметрическая положительно определенная, но демпфирование можно считать малым. То есть матрица из (94) теперь недиагональная, а е? элементы малы по сравнению с диагональными компонентами матрицы . Систему (90) в этом случае можно записать в следующем виде

Е Mq + Dq + Kq = 0,

(99)

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


возмущений несамосопряженных операторов [43] эти приращения выражаются в виде рядов по целым или дробным степеням в зависимости от жордановой структуры собственного значения 0 . Пусть 0 простое собственное значение системы (99) при = 0, а u
0

соответствующий собственный вектор. В этом случае приращения для

0 и u0 принимают вид [43] = 0 + 1 + 2 2 + ћ ћ ћ u = u0 + u1 + u2 + ћ ћ ћ
Известно [44, 45, 46, 47], что
2

(100)

0 = +i
чисто мнимая величина, тогда как

uT Ku0 0 T Mu u0 0

uT Du0 1 = - 0 2uT Mu 0

0

действительная отрицательная величина. Следовательно, малое демпфирование изменяет частоту системы (24) на величину O(2 ). Пусть все собственные числа системы (90) простые. Следовательно, существует ровно n линейно независимых пар комплексно сопряженных собственных векторов. Из n собственных векторов u системы можно составить квадратную невырожденную матрицу U. Перепишем разложения (100) в матричном виде

= + i - + ћ ћ ћ U = U0 + U1 + ћ ћ ћ
94

,

(101)


где диагональная матрица, по диагонали которой расположены собственные числа , U0 действительная матрица, составленная из собственных векторов задачи без демпфирования, т. е. при = 0. Для представленных во второй главе алгоритмов определения моментов инерции нужна матрица U0 . Покажем, что в случае малого демпфирования матрица U1 чисто мнимая и поэтому можно определить матрицу U0 , зная матрицу U. Задачи на определение собственных векторов можно записать в матричном виде

MU2 + DU + KU = 0
и при = 0

(102)

KU0 - MU0 2 = 0.

(103)

Подставив в матричное уравнение (102) разложения (101), отбросив члены порядка больше и учтя уравнение (103), получим

KU1 - MU1 2 = +i(2MU0 - DU0 ).
Умножив это уравнение слева на (UT MU0 )-1 UT и учтя, что 0 0

(104)

(UT MU0 )-1 UT KU0 = 2 , 0 0
получаем уравнение Сильвестра [49] относительно матрицы U-1 U 0
1

2 U-1 U1 - U-1 U1 2 = +i(2 - (UT MU0 )-1 UT DU0 ). 0 0 0 0

(105)

Это матричное уравнение можно решить, рассматривая его по компонентам

(U-1 U1 )j k = 0

i((UT MU0 )-1 UT DU0 )j 0 0
95

k

k 2, - k
2 j

k=j

(106)


Для нахождения диагональных элементов матрицы U-1 U1 воспользуемся 0 условиями нормировки uT u1 = 0 для приращений собственных векторов. 0 Запишем эти условия, используя матрицы U0 и U
1

(UT U1 )kk = 0, 0

k = 1, . . . , n.

Квадратная матрица U0 невырожденная, следовательно можно записать предыдущую формулу следующим образом

(UT U0 U-1 U1 )kk = 0, 0 0
Рассматривая произведения матриц UT U0 и U-1 U1 по компонентам, по0 0 лучаем
n

(UT U0 )kj (U-1 U1 )j k = 0. 0 0
j =1

Из этого уравнения выражаются диагональные элементы матрицы U-1 U 0 через е? недиагональные элементы
n k j =1 j =k

1

(U U1 )k

-1 0

k

-1 = (UT U0 )k 0

(UT U0 )kj (U-1 U1 )j k . 0 0

Подставив в это уравнение выражение (106) для недиагональных элементов матрицы U-1 U1 , получаем 0

(U U1 )k

-1 0

k

+i = (UT U0 ) 0

n kk j =1 j =k

(UT U0 )kj ((UT MU0 )-1 UT DU0 )j 0 0 0

k

k 2 2. j - k

Из решения видно, что приращение действительных собственных векторов U0 является чисто мнимым порядка , то есть матрица U1 чисто мнимая. Это согласуется с известными результатами [48]. Зная вектор u, можно получить вектор u0 с точностью до действительного множителя.
96


Собственные векторы определены с точностью до комплексного множителя. Пусть

u = c(u0 + iu1 ),
где c комплексная константа. Составим вектор из модулей компонент вектора u

|uj | = |c| u2j + 2 u2j |c||u0j |. 0 1
Таким образом, в случае малого демпфирования можно получить собственные векторы задачи без демпфирования, необходимые для алгоритмов определения моментов инерции, представленных во второй главе.

97


Основные результаты и выводы
В первой главе диссертации представлен краткий обзор существующих методов измерения моментов инерции. Описана конструкция измерительного стенда, спроектированного в ЦАГИ. Дано математическое описание динамической системы стенда и поставлена задача определения моментов инерции тела. Во второй главе предложено три варианта решения задачи в зависимости от дополнительных сведений о способе возбуждения колебаний или о параметрах тела и жесткостях пружин стенда. Рассмотрены методы Прони и Кунга идентификации линейной динамической системы по последовательности наблюдаемого сигнала. Проведено сравнение этих двух методов, в результате которого предпочтение отдано методу Кунга. Получены явные выражения для определения моментов инерции тела по результатам идентификации с привлечением дополнительных условий в трех вариантах. Показано, что разработанный математический аппарат совместим с различными методами идентификации линейных систем. Дана оценка минимального времени возбуждения колебаний многомерной линейной динамической системы малой управляющей силой. В третьей главе приведены результаты численных решений полных уравнений системы, учитывающих нелинейности. При помощи этих решений проведен анализ чувствительности определяемых параметров тела к погрешности измеряемого датчиками сигнала. Проведен анализ чувствительности определяемых параметров к ошибке идентификации. Най98

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

99


Список литературы
[1] Буянов Е. В. Методика и установка для точного определения тензора инерции твердого тела // Измерительная техника. 1988. 12. С. 25 27. [2] Гернет М. М., Ротобыльский В. Ф. Определение моментов инерции. М.: "Машиностроение", 1969. [3] Orne D., Schmitz T. Analysis of a Platform for Measuring Moments and Products of Inertia of Large Vehicles //Journal of dynamic systems, measurement and control, 2, 1978. [4] Ashley S. Testing vehicle inertia //Mechanical Engineering, 117, 1995. [5] Мещерский И. В. Сборник задач по теоретической механике. М.: "Наука". Гл. ред. физ.-мат. лит., 1981. 480 с. [6] Ишлинский А. Ю., Стороженко В. А., Темченко М. Е. Вращение твердого тела на струне и смежные задачи. М.: Наука, 1991. [7] Мельников В. Г. Синтез и исследование нелинейных систем управления для параметрической идентификации тензоров инерции тел. Диссертация на соискание ученой степени канд. техн. наук, С.-Петербург, 2001. [8] Hecker F., Hahn H. Mathematical Modeling and Parameter Identication of a Planar Servo-Pneumatic Test Facility //Nonlinear Dynamics. 1997. 14, С. 269277.
100


[9] Богданов В. В., Волобуев В. С., Кудряшов А. И. Комплекс для измерения массы, координат центра масс и моментов инерции машиностроительных изделий // Измерительная техника, 2002. 2. С. 37-39. [10] Новожилов И. В. Фракционный анализ. М.: Изд-во Моск. ун-та, 1995. 224 с. [11] Ишлинский А. Ю. Механика гироскопических систем. М.: Изд-во АН СССР, 1963. [12] Gladwel l G. M. L. Isospectral vibrating beams // Procedings of Royal Society, London. A 458, pp. 2691-2703. 2002. [13] Kung S. Y. A New Identication and Model Reduction Algorithm via singular value decomposition//12th Asilomar Conf. Circuits, Syst. Comput. Pacic Grove. Calif., Nov. 1978. [14] Verhaegen M. Identication of the Deterministic Part of MIMO State Space Models Given in Innovation Forms from Input-Output Data // Automatica, V.30. 1, pp. 61-74. 1994. [15] Viberg M. Subspace-based Methods for the Identication of Linear Timeinvariant Systems // Automatica, V.31. 12, pp. 18351851. 1995. [16] Александров В. В., Садовничий В. А, Чугунов О. Д. Математические задачи динамической имитации полета. М.: Изд-во Моск. ун-та, 1986. 181 с.

101


[17] Федорова Г. А. Идентификация параметров имитационных динамических стендов. Диссертация на соискание ученой степени канд. физ.мат. наук, Москва, 1992. [18] Цыпкин Я. З. Информационная теория идентификации. М.: Наука, 1995. [19] Эйкхофф П. Основы идентификации систем управления. М.: Наука, 1975. [20] Эйкхофф П. Современные методы идентификации систем. М.: Мир, 1983. [21] Льюнг Л. Идентификация систем. Теория для пользователя: Пер. с англ. / Под ред. Я. З. Цыпкина. М.: Наука, 1991. 432 с. [22] Ефимов Н. В., Розендорн Э. Р. Линейная алгебра и многомерная геометрия. М.: Наука, 1974. С. 357361. [23] Журавлев В. Ф., Климов Д. М. Прикладные методы в теории колебаний. М.: Наука, 1983. 328 с. [24] Моисеев Н. Н. Элементы теории оптимальных систем. М.: Наука, 1975. 528 с. [25] Маркеев А. П. Теоретическая механика. М.: Наука. Гл. ред. физ.-мат. лит., 1990. 416 с. [26] Марпл-мл. С. Л. Цифровой спектральный анализ и его приложения. М.: Мир, 1990.
102


[27] Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Наука, 1978. 600 с. [28] Годунов С. К., Рябенький В. С. Разностные схемы. Введение в теорию. М.: Наука, 1977. 439 с. [29] Темников Ф. Е., Афонин В. А., Дмитриев В. И. Теоретические основы информационной техники., М.: Энергия, 1971. 423 с. [30] Беляков А. О. Численное моделирование процесса измерения моментов инерции крупногабаритных тел методом свободных колебаний // Ученые записки ЦАГИ. 2002. 12, С. 129139. [31] Беляков А. О. Определение динамических параметров массивных тел по формам колебаний. Вестник молодых ученых. Серия: Прикладная математика и механика С.-Петербург (в печати). [32] Беляков А. О. Optimal excitation of oscillations by a limited control force // Труды международной конференции Физика и управление С.-Петербург. 2003. С. 11301133. [33] Беляков А. О., Блаженнова-Микулич Л. Ю. Идентификация инерционной матрицы консервативной колебательной системы // Вестн. Моск. ун-та. сер. 1, Математика механика. 2005. 3, С. 2528. [34] Bushaw D. W. Experimental towing tank // Stevens Inst. of Technology. Reprint 169. N.Y.: Hoboken, 1953.

103


[35] Александров В. В., Болтянский В. Г., Лемак С. С., Парусни-

ков Н. А., Тихомиров В. М. Оптимизация динамики управляемых
систем. М.: Изд-во Моск. ун-та, 2000. 304 с. [36] Понтрягин Л. С., Болтянский В. Г., Гамкрелидзе Р. В., Мищен-

ко Е. Ф. Математическая теория оптимальных процессов. М.: Наука,
Гл. ред. физ.-мат. лит., 1969. 384 с. [37] Саввин А. Б. О наибыстрейшем выведении изображающей точки за пределы заданной фазовой плоскости. // Изв. АН СССР. Технич. ктберн. 1963. 4. [38] Обморшев А. Н. Введение в теорию колебаний. М.: Наука, 1965. 276 с. [39] Магнус К. Колебания. Введение в исследование колебательных систем. М.: Мир, 1982. 304 с. [40] Пановко Я. Г., Губанова И. И. Устойчивость и колебания упругих систем. Современные концепции, парадоксы и ошибки. М.: Наука, 1987. 352 с. [41] Вибрации в технике. Справочник. Том 1. Колебания линейных систем /Под редакцией В. В. Болотина. М.: Машиностроение, 1999. 504 с. [42] Трубецков Д. И., Рожнев А. Г. Линейные колебания и волны. М.: Физматлит, 2001. 416 с.

104


[43] Вишик М. И., Люстерник Л. А. Решение некоторых задач о возмущении в случае матриц и самосопряженных и несамосопряженных дифференциальных уравнений. I. УМН. 1960. Т. 15. Вып. 3. С. 380. [44] Севрюк М. Б., Сейранян А. П. Эволюция частот колебаний диссипативной системы // ПММ. 1993. Т. 57. Вып. 4. С. 2130. [45] Сейранян А. П., Шаранюк А. В. Анализ чувствительности частот колебаний механических систем // Изв. АН СССР. МТТ. 1987. 2. С. 3741. [46] Seyranian A. P., Elishako I. (Eds.) Modern Problems of Structural Stability. Wien, New York: Springer. 2002. 394 p. [47] Seyranian A. P., Mailybaev A. A. Multiparameter Stability Theory with Mechanical Applications. Singapore: World Scientic, 2003. 420 p. [48] Гантмахер Ф. Р. Лекции по аналитической механике. М.: Физматлит, 1966. [49] Гантмахер Ф. Р. Теория матриц. М.: Наука. 1988. [50] Амензаде Ю. А. Теория упругости. М.: Высшая школа. 1976. [51] Ганиев Р. Ф., Кононенко В. О. Колебания твердых тел. М.: "Наука". Гл. ред. физ.-мат. лит., 1976. 432 с.

105


Приложение

> restart; , > f:=<0,P,0,0>; 0 P f := 0 0 > C:=Matrix(<<1,1,1,1>|<-rx,Lx-rx,Lx-rx,-rx>|>); 1 rx 1 Lx rx C := 1 Lx rx 1 rx rz rz Lz rz Lz rz , ,

> Y:=matrix(4,3); , > G:=matrix(4,3); Y := array( 1 .. 4, 1 .. 3, [ ] ) Y

G := array( 1 .. 4, 1 .. 3, [ ] ) ( > MM:=Matrix(3,(i) -> sum('Y[k,i]*f[k]', 'k'=1..4)/w[i]^2,shape=diagonal); Y2, 1 P w MM :=
2 1

)

0 Y2, 2 P w2
2

0

0

0 Y2, 3 P w
2 3

0

0

> with(linalg): M:=multiply(transpose(C),G,MM,transpose(G),C):
Warning, the protected names norm and trace have been redefined and unprotected

, > mc:=solve({factor(M[1,2])=0, factor(M[1,3])=0},{rx, rz});

106


mc := { rx

Lx ( Y2, 2 G2,
2 2

2 2

w1 w3

2

2

Y2, 1 G2
2, 1 2 3, 2

2 ,1 2

w2 w
2

2

2 3 1

Y2, 3 G1, 3 w1 w2 G2
2 2

2

2

,3 ,1 3, 2

Y2, 3 G2, 3 w1 w2 G Y2, 1 G2, 1 w2 w3 G
2 2 2 2

4, 3 4, 1 3, 3

2 Y2, 1 G Y2, 2 G

w2 w3 G3,
2 2 2

Y2, 1 G1, 1 w2 w3 G3
2 ,3

w1 w3
2

Y2, 3 G2
3

w1 w2
2 3

2

2 2

Y2, 2 G
2

w1 w3 G4,

2

2

2

2 Y2, 3 G2, 3 w1 w2 G Y2, 2 G1, 2 w1 w3 G
2 2 2 2

Y2, 3 G Y2, 3 G

1, 3

w1 w2 G3,
2 2 3

Y2, 3 G3, Y2, 1 G3,
2 1

w1 w2
2 2 2

2, 2 3, 2

3, 3

w1 w2 G4,
2 2

w2 w3
2

2 Y2, 2 G2, 2 w1 w3 G Y2, 2 G1, 2 w1 w3 G Y2, 2 G1 Y2, 3 G2
2 ,2 2 ,3 2 2

Y2, 1 G Y2, 1 G
2 ,3 2 ,1

3, 1

w2 w3 G4,
2 2

1

Y2, 2 G2, 2 w1 w3 G4 ( Y2, 1 G1
2 2 ,1

,2

3, 2

1, 1

w2 w3 G2, 1 )
2 2

w2 w

2

2 3 2 2 2 2

w1 w w1 w
2 2

2

2 3 2 2 2

Y2, 3 G1 Y2, 1 G3
2, 1 2, 2 2, 3 3, 1 4, 2 4, 1

w1 w2 w2 w3
1, 1 1, 2 1, 3 2, 1 2, 3 3, 2 2

2

Y2, 1 G2 Y2, 2 G3
2 2

,1 2 ,2 1 2 3 1 3 2 2 ,3

w2 w w1 w
2

2

2 3 2 3

Y2, 2 G2 Y2, 3 G3
2

,2 2 ,3

w1 w3 w1 w2
,1 ,2 ,3 ,2 ,3 ,3 2

2 Y2, 1 G1, 1 w2 w3 G 2 Y2, 2 G1, 2 w1 w3 G 2 Y2, 3 G1, 3 w1 w2 G 2 Y2, 1 G2, 1 w2 w3 G 2 Y2, 2 G2, 2 w1 w3 G 2 Y2, 1 G3, 1 w2 w3 G Y2, 1 G4
2 ,1 2 2 2 2 2 2 2 2 2 2 2

2 Y2, 1 G 2 Y2, 2 G 2 Y2, 3 G 2 Y2, 1 G 2 Y2, 3 G 2 Y2, 2 G
2 ,2 4, 1 2

w2 w3 G3, w1 w3 G3, w1 w2 G3, w2 w3 G4, w1 w2 G3, w1 w3 G4,
2 2 2 2 2 2 2 2 2 2 2 2 2 2

2 Y2, 1 G1, 1 w2 w3 G4 2 Y2, 2 G1, 2 w1 w3 G4 2 Y2, 3 G1, 3 w1 w2 G4 2 Y2, 2 G2, 2 w1 w3 G3 2 Y2, 3 G2, 3 w1 w2 G4 2 Y2, 3 G3, 3 w1 w2 G4 w1 w2 ), rz
2 2 ,2 2 2 2 2 2 2 2 2 2 2 2 2

2

w2 w
2

2

2 3 4, 1

Y2, 2 G4

w1 w3
2 2

Y2, 3 G4

Lz ( Y2, 2 G1, 2 w1 w3 G4 Y2, 3 G
,2 2 ,1 1, 1 1, 3 2 2 ,2

Y2, 1 G2, 1 w2 w3 G
2

Y2, 1 G
4, 2 4, 3 4, 3 4, 2

w2 w3
2 2

Y2, 2 G2, 2 w1 w3 G3 Y2, 1 G3
2 ,1

Y2, 2 G2, 2 w1 w3 G Y2, 3 G2, 3 w1 w2 G
2 2 2 2

2

Y2, 2 G Y2, 3 G

4, 2 2 4, 3

w1 w3 w1 w2
2 3, 1 2 2

w2 w3
2 2

2

2

w1 w2 G4,

2

2

3

Y2, 2 G1, 2 w1 w3 G3
2 1

2 Y2, 3 G3, 3 w1 w2 G 2 Y2, 2 G3, 2 w1 w3 G Y2, 1 G2, 1 w2 w3 G Y2, 1 G1
2 ,1 2 2 2 2

2 Y2, 1 G Y2, 3 G Y2, 3 G
2 ,2

w2 w3 G4,
2 2

Y2, 1 G1, 1 w2 w3 G4
2 ,2

2

3, 3

w1 w2
2 2

Y2, 2 G3
3 2 ,3

w1 w3
2

2

2 2

Y2, 1 G

w2 w3 G3, (

2

2

1

3, 1

1, 3

w1 w2 G3,
2

Y2, 3 G2, 3 w1 w2 G3, 3 ) w1 w
2 ,1 2 2 2 2 3

w2 w
2 ,2 2 ,3

2

2 3 2 3 2 2

Y2, 2 G1

w1 w3
2 ,3 2

2

Y2, 3 G1
2 2

Y2, 1 G2

2 ,1

w2 w3
2 ,2 ,1 2 2

2

2 2

Y2, 2 G2 Y2, 3 G3

w1 w w1 w
2 2

2

Y2, 3 G2

w1 w2
2

Y2, 1 G3
,1 2

w2 w

2

Y2, 2 G3
2 2

w1 w3

2 Y2, 1 G1, 1 w2 w3 G2
2 4, 1 4, 2 4, 3

2 Y2, 1 G1, 1 w2 w3 G3
2 2 3 1

2 Y2, 1 G1, 1 w2 w3 G 2 Y2, 2 G1, 2 w1 w3 G 2 Y2, 3 G1, 3 w1 w2 G
2 2 2 2

2 Y2, 2 G 2 Y2, 3 G 2 Y2, 1 G

1, 2 1, 3 2, 1

w1 w3 G2, w1 w2 G2, w2 w3 G3,
2 2 2 2

2

2 Y2, 2 G1, 2 w1 w3 G3 2 Y2, 3 G1, 3 w1 w2 G3 2 Y2, 1 G2, 1 w2 w3 G4
2 2 2 2

,2 ,3 ,1

107


2 Y2, 2 G2, 2 w1 w3 G 2 Y2, 3 G2, 3 w1 w2 G 2 Y2, 3 G3, 3 w1 w2 G
2 2 2 2

2

2

3, 2 4, 3 4, 3

2 Y2, 2 G 2 Y2, 1 G Y2, 1 G

2, 2 3, 1 2

w1 w3 G4, w2 w3 G4,
2 2 2 2

2

2

2 1

2 Y2, 3 G2, 3 w1 w2 G3 2 Y2, 2 G3, 2 w1 w3 G4
2 ,2 2 2

2

2

,3 ,2 2

4, 1

w2 w3

Y2, 2 G4

w1 w3

2

2

Y2, 3 G

4, 3

w1 w2 ) }

2

2

> Mass:=normal(M[1,1]); Mass := P ( Y2, 1 G1 Y2, 2 G2 Y2, 3 G3
2 ,2 2 ,3 2 ,1 2

w2 w3
2 3 2 2 2

2

2

Y2, 2 G1
2 ,3

2 ,2 2

w1 w
2 2

2

2 3

Y2, 3 G1
2 ,1

2 ,3 2

w1 w2
2 3

2

2

Y2, 1 G2
2 ,2 ,1 2

2 ,1 2

w2 w
2

2

2 3

w1 w w1 w
2 2

Y2, 3 G2

w1 w2
2

Y2, 1 G3
,1 2

w2 w

Y2, 2 G3
2 2

w1 w3

2 Y2, 1 G1, 1 w2 w3 G2
4, 1 4, 2 4, 3 3, 2 4, 3 4, 3

2 Y2, 1 G1, 1 w2 w3 G3
2 2 3 1 2 1

2 Y2, 1 G1, 1 w2 w3 G 2 Y2, 2 G1, 2 w1 w3 G 2 Y2, 3 G1, 3 w1 w2 G 2 Y2, 2 G2, 2 w1 w3 G 2 Y2, 3 G2, 3 w1 w2 G 2 Y2, 3 G3, 3 w1 w2 G ( w1 w2 w3 )
2 2 2 2 2 2 2 2 2 2 2 2 2

2 Y2, 2 G 2 Y2, 3 G 2 Y2, 1 G 2 Y2, 2 G 2 Y2, 1 G Y2, 1 G

1, 2 1, 3 2, 1 2, 2 3, 1 2

w1 w3 G2, w1 w2 G2, w2 w3 G3, w1 w3 G4, w2 w3 G4,
2 2 2 2 2 2 2 2 2 2

2

2 Y2, 2 G1, 2 w1 w3 G3 2 Y2, 3 G1, 3 w1 w2 G3 2 Y2, 1 G2, 1 w2 w3 G4 2 Y2, 3 G2, 3 w1 w2 G3 2 Y2, 2 G3, 2 w1 w3 G4
2 ,2 2 2 2 2 2 2 2 2

,2 ,3 ,1 ,3 ,2 2

4, 1

w2 w3

Y2, 2 G4

w1 w3

2

2

Y2, 3 G

4, 3

w1 w2 )

2

2

> Izz:=factor(subs(mc, collect((normal(M[2,2])), {rx, rz}))); > Ixx:=factor(subs(mc, collect((normal(M[3,3])), {rx, rz}))); > Izx:=factor(subs(mc, collect((normal(M[2,3])), {rx, rz}))); Izz := P Lx2 ( 2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G1, 2 G3 2 Y2, 3 G3, 3 w1 G4, 3 Y2, 2 G2, 2 G4 2 Y2, 3 G3, 3 w2 G4, 3 Y2, 1 G1, 1 G2 2 Y2, 2 G1, 2 w3 G2, 2 Y2, 1 G3, 1 G4 2 Y2, 3 G1, 3 w2 G3, 3 Y2, 1 G3, 1 G4 2 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G1, 2 G3 2 Y2, 2 G1, 2 w3 G3, 2 Y2, 1 G1, 1 G2 2 Y2, 1 G3, 1 w3 G4, 1 Y2, 2 G2, 2 G4 2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G2, 1 G4 2 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G1, 2 G3 2 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G1, 1 G2
2 2 2 2 2 2 2 2 2 2 ,2 ,1 ,1 ,1 ,2 ,1 ,2 ,1 ,2 ,1 2 ,2

2 Y2, 3 G3, 3 w2 G4, 3 Y2, 1 G3, 1 G4
2 ,2 ,3 ,1 ,2 ,1 ,1 ,2 ,2 ,1 ,2

2

,1

2 Y2, 3 G3, 3 w1 G4, 3 Y2, 2 G1, 2 G3 2 Y2, 2 G1, 2 w1 G2, 2 Y2, 3 G3, 3 G4 2 Y2, 2 G1, 2 w3 G2, 2 Y2, 1 G1, 1 G2 2 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G2, 2 G4 2 Y2, 3 G1, 3 w2 G3, 3 Y2, 1 G1, 1 G2 2 Y2, 2 G2, 2 w3 G4, 2 Y2, 1 G1, 1 G2 2 Y2, 1 G3, 1 w3 G4, 1 Y2, 2 G1, 2 G3 2 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G2, 2 G4 2 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G3, 1 G4 2 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G3, 2 G4
2 2 2 2 2 2 2 2 2

108


2 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G1, 1 G3 2 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G1, 2 G2 2 Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G3, 3 G4 2 Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G1, 3 G3 2 Y2, 2 G3, 2 w3 G4, 2 Y2, 1 G3, 1 G4 2 Y2, 2 G3, 2 w1 G4, 2 Y2, 3 G1, 3 G3 2 Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G3, 3 G4 2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G3, 2 G4 2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G2, 2 G4 2 Y2, 2 G3 Y2, 3 G2
,3 2 ,2 2 2 2 2 2 2 2 2

2

,1 ,2 ,3 ,3 ,1 ,3 ,3 ,2 ,2

2 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G2, 1 G4 2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G1, 2 G2 2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G3, 2 G4 2 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G1, 2 G2 2 Y2, 2 G3, 2 w1 G4, 2 Y2, 3 G3, 3 G4 2 Y2, 2 G3, 2 w3 G4, 2 Y2, 1 G1, 1 G2 2 Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G1, 3 G3 2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G1, 2 G2 2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G1, 2 G3
2 3, 2 2 2 2 2 2 2 2 2 2

2

,1 ,2 ,2 ,2 ,3 ,1 ,3 ,2 ,2

w1 Y2, 3 G1, 3 G
2 2 ,1

2

4, 3

Y2, 2 G

w3 Y2, 1 G
3, 2 2

2

2 4, 1 2 4, 1 2 1, 1

Y2, 2 G3

2 ,2

w1 Y2, 3 G4

2

2 ,3

2

w2 Y2, 1 G1
2

2 Y2, 2 G2, 2 w3 G
2 ,3

Y2, 1 G

2 Y2, 2 G2, 2 w1 G3, 2 Y2, 3 G4 Y2, 3 G3 Y2, 1 G2 Y2, 3 G2
2 ,3 2 ,1 2 ,3

Y2, 3 G
2 ,3

2 3, 3 2

w2 Y2, 1 G
2 ,2

Y2, 3 G3
2 ,1

2 ,3

w1 Y2, 2 G1
2 ,2

2

2 ,2

w2 Y2, 1 G4 w2 Y2, 3 G1 w1 Y2, 2 G4
2 2 2

2

2 ,1 2 ,3 2 ,2

Y2, 3 G3

w1 Y2, 2 G4
2 2

Y2, 1 G2
,2 2 1, 1

w3 Y2, 2 G1
2 ,3 2

2

2 Y2, 1 G2

,1

w3 Y2, 2 G1, 2 G4
2 3, 2

Y2, 3 G2

w2 Y2, 1 G4

2 ,1

2 Y2, 2 G2, 2 w3 G
2 ,3

Y2, 1 G
2

2 Y2, 2 G2, 2 w1 G3, 2 Y2, 3 G1 Y2, 3 G2
2 ,3

4 Y2, 2 G
2 ,3

2, 2 2

w3 G3, 2 Y2, 1 G1, 1 G4
,1 2 1, 1

,1

w1 Y2, 2 G1
2 2

2

2 ,2

2 Y2, 3 G2
4, 2 4, 1

w2 Y2, 1 G1, 1 G4
2 3, 2

2 Y2, 3 G2 2 Y2, 2 G3

,3 2 ,2

w1 Y2, 2 G1, 2 G w3 Y2, 1 G1, 1 G
2 2

Y2, 2 G

w3 Y2, 1 G
2

2

Y2, 2 G3

2 ,2

w1 Y2, 3 G1

2

2 ,3

4 Y2, 2 G
,3

2, 2

w1 G3, 2 Y2, 3 G1, 3 G4
2

,3 2

4 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G1, 3 G4 2 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G4 Y2, 1 G3
2 ,1 2 2 ,3

2 Y2, 1 G2, 1 w3 G3, 1 Y2, 2 G4
2 2, 1

,2 ,2 2 2 ,2

4 Y2, 1 G
2 ,1 2

w3 G3, 1 Y2, 2 G1, 2 G4
2 ,3

w3 Y2, 2 G4
2

2

2 ,2

Y2, 1 G3
2 ,3 4, 3

w2 Y2, 3 G4
2 3, 1 2, 3

2 Y2, 1 G2, 1 w3 G3, 1 Y2, 2 G1
1, 2

2 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G1 2 Y2, 1 G3
2 ,1 2

2 Y2, 1 G 4 Y2, 3 G
,2

w3 Y2, 2 G
2

2

G4,

2 ,1 2 ,1 2 ,1 2 ,1 2 2 2 ,3 2 ,3 4, 1

w2 Y2, 3 G1, 3 G
2

w2 G3, 3 Y2, 1 G1, 1 G4
2

4 Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G1, 2 G4 2 Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G4 2 Y2, 1 G2 Y2, 2 G2
,2 2 ,1 2 2 2 ,2 4, 3

2 Y2, 3 G2, 3 w2 G3, 3 Y2, 1 G4
2 3, 1 2 2, 1 2

Y2, 1 G Y2, 1 G
2 ,2

w3 Y2, 2 G w3 Y2, 2 G
2 ,3 2

2

2 1, 2 2 4, 2

Y2, 1 G3 Y2, 1 G2
2 ,2

w2 Y2, 3 G1 w2 Y2, 3 G4
2

w2 Y2, 3 G1, 3 G
2 2 ,1

2

w3 Y2, 1 G1

Y2, 2 G2

w1 Y2, 3 G1

2 Y2, 2 G2

w3 Y2, 1 G1, 1 G

109


2 Y2, 2 G2

2 ,2

w1 Y2, 3 G1, 3 G
2

2

4, 3 2

Y2, 2 G

2 2, 2

w3 Y2, 1 G
2

2

2 4, 1

Y2, 2 G2
2 ,2 2

2 ,2

w1 Y2, 3 G4

2

2 ,3

2 Y2, 3 G2, 3 w2 G3, 3 Y2, 1 G1 2 Y2, 3 G3
2 ,3 2

,1 4, 1

2 Y2, 3 G 2 Y2, 3 G
,1 ,2 ,2 ,1 2

2, 3 3, 3

w1 G3, 3 Y2, 2 G1
2

w2 Y2, 1 G1, 1 G
2

w1 Y2, 2 G
2

2

1, 2

G4,

2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G3, 1 G4 2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G1, 2 G2 2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G1, 2 G3 2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G1, 1 G3 Y2, 1 G1
2 ,1 2 2 2

2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G2, 2 G4 2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G3, 2 G4 2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G1, 1 G2
2 2 2

,2 ,2 ,1

2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G2, 2 G4, 2 ) Y2, 3 G1
2 2 2 ,3

(
2 2

w2 w
2 ,2 2 ,3

2

2 3 2 3 2 2

Y2, 2 G1

2 ,2

w1 w3
2 ,3 2

2

w1 w
2 ,1

2

2 2 2 3

Y2, 1 G2

2 ,1

w2 w3
2 ,2 ,1 2 2

2

Y2, 2 G2 Y2, 3 G3

w1 w w1 w
2 2

2

Y2, 3 G2

w1 w2
2

Y2, 1 G3
,1 2

w2 w

2

Y2, 2 G3
2 2

w1 w3

2 Y2, 1 G1, 1 w2 w3 G2
2 4, 1 4, 2 4, 3 3, 2 4, 3 4, 3 2

2 Y2, 1 G1, 1 w2 w3 G3
2 2 3 1 2 1

2 Y2, 1 G1, 1 w2 w3 G 2 Y2, 2 G1, 2 w1 w3 G 2 Y2, 3 G1, 3 w1 w2 G 2 Y2, 2 G2, 2 w1 w3 G 2 Y2, 3 G2, 3 w1 w2 G 2 Y2, 3 G3, 3 w1 w2 G
2 2 2 2 2 2 2 2 2 2

2 Y2, 2 G 2 Y2, 3 G 2 Y2, 1 G 2 Y2, 2 G 2 Y2, 1 G Y2, 1 G

1, 2 1, 3 2, 1 2, 2 3, 1 2

w1 w3 G2, w1 w2 G2, w2 w3 G3, w1 w3 G4, w2 w3 G4,
2 2 2 2 2 2 2 2 2 2

2

2 Y2, 2 G1, 2 w1 w3 G3 2 Y2, 3 G1, 3 w1 w2 G3 2 Y2, 1 G2, 1 w2 w3 G4 2 Y2, 3 G2, 3 w1 w2 G3 2 Y2, 2 G3, 2 w1 w3 G4
2 ,2 2 2 2 2 2 2 2 2

,2 ,3 ,1 ,3 ,2 2

4, 1

w2 w3
,2

Y2, 2 G4

w1 w3
2

2

2

Y2, 3 G

4, 3

w1 w2 )
,1

2

2

Ixx := P Lz2 ( 2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G1, 2 G3 4 Y2, 2 G1, 2 w1 G2, 2 Y2, 3 G3, 3 G4 2 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G2, 2 G4 2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G2, 1 G4 4 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G3, 2 G4 4 Y2, 2 G3, 2 w3 G4, 2 Y2, 1 G1, 1 G2 2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G2, 2 G4 Y2, 3 G3 Y2, 3 G2 Y2, 2 G3
2 ,3 2 ,3 2 ,2 2 2 2 2 2 2 ,3 ,2 ,1 ,2 ,1 ,2 2 ,3

4 Y2, 3 G3, 3 w2 G4, 3 Y2, 1 G1, 1 G2
2 ,1 ,2 ,1 ,3 ,3 ,2 2 ,1 2

4 Y2, 2 G1, 2 w3 G2, 2 Y2, 1 G3, 1 G4 2 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G1, 2 G3 4 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G3, 1 G4 2 Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G1, 3 G3 2 Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G1, 3 G3 2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G1, 2 G3 w1 Y2, 2 G1
2 3, 2 3, 2 2 2 ,2 2 2 2 2 2

w2 Y2, 1 G1 w1 Y2, 2 G4 w1 Y2, 3 G1
2 2 2

2

2 ,1 2 ,2 2 ,3

Y2, 3 G3

Y2, 3 G2
1, 1 1, 3

2 ,3

w2 Y2, 1 G4 Y2, 2 G3,
2 2

2

2 Y2, 2 G2, 2 w3 G 2 Y2, 2 G2, 2 w1 G
,3 ,1 2

Y2, 1 G Y2, 3 G
2

G4, G4,

1 3

w3 Y2, 1 G1,

2 1

2 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G1, 3 G4 2 Y2, 3 G2, 3 w2 G3, 3 Y2, 1 G1, 1 G4
2

2 Y2, 1 G2, 1 w3 G3, 1 Y2, 2 G1, 2 G4 2 Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G1, 2 G4
2

,2 ,2

110


Y2, 1 G3 Y2, 1 G2

2 ,1 2 ,1

w3 Y2, 2 G1 w2 Y2, 3 G4
2 2

2

2 ,2 2 ,3

Y2, 1 G3 Y2, 2 G2

2 ,1 2 ,2

w2 Y2, 3 G1 w3 Y2, 1 G4
2

2

2 ,3 2 ,1 2

Y2, 1 G2 Y2, 2 G2

2 ,1 2 ,2

w3 Y2, 2 G4 w1 Y2, 3 G4
,2 ,2 ,2 2

2

2 ,2 2 ,3

2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G2, 2 G4 2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G1, 1 G3 2 Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G2, 2 G3 2 Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G1, 2 G3 2 Y2, 3 G3, 3 w1 G4, 3 Y2, 2 G1 2 Y2, 3 G2 2 Y2, 2 G3
2 ,3 2 ,2 2 2 2 ,2 4, 2 2, 1 2 ,3 2 2 2

,2 ,1 ,2 ,2

2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G1, 2 G3 2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G2, 2 G4 2 Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G2, 2 G4 2 Y2, 3 G3, 3 w2 G4, 3 Y2, 1 G1
2 2, 3 2 2, 3 2 2 ,1 2 2

2 Y2, 3 G Y2, 3 G

w2 Y2, 1 G
2

2

3, 1 2

G4,

1

w1 Y2, 2 G3, 2 G w3 Y2, 1 G1, 1 G
2 2

w2 Y2, 1 G
2

3, 1 2 ,1

2 Y2, 2 G Y2, 2 G
,2

3, 2 2

w3 G4, 2 Y2, 1 G1
2 2 2, 3

2 Y2, 2 G3, 2 w1 G4, 2 Y2, 3 G1
2

3, 2

w1 Y2, 3 G
2

2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G2, 2 G3 2 Y2, 2 G1, 2 w3 G2, 2 Y2, 1 G3
2 2 2 ,1

2 Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G2, 3 G3
2 3, 3

,3

2 Y2, 3 G
,2

w1 Y2, 2 G
2

2

1, 2

G2,

2 2 ,2 2 ,3 2 2 ,2

2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G2, 2 G3 2 Y2, 1 G3, 1 w2 G4, 1 Y2, 3 G1 2 Y2, 3 G4 2 Y2, 1 G4
2 ,3 2 ,1 2 2 2 ,3 2, 1 2, 3

2 Y2, 1 G3, 1 w3 G4, 1 Y2, 2 G1
2 4, 3

Y2, 3 G

w2 Y2, 1 G
2 2

2

2 1, 1 1, 2 2 1, 2

Y2, 3 G4 G2,
2

w1 Y2, 2 G1

w2 Y2, 1 G1, 1 G w2 Y2, 3 G1, 3 G
2 2

2 Y2, 1 G Y2, 1 G
,2 ,3 ,3 ,1 2 ,2

4, 1 2

w3 Y2, 2 G
2

4, 1

w3 Y2, 2 G
2

Y2, 1 G4

2 ,1

w2 Y2, 3 G1
,2 ,3 ,3

2

2 ,3

2 Y2, 1 G2, 1 w3 G3, 1 Y2, 2 G1, 2 G3 2 Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G2, 3 G3 2 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G1, 3 G3 2 Y2, 2 G1, 2 w3 G3, 2 Y2, 1 G1, 1 G4 Y2, 2 G4
2 ,2 2 2 2

2 Y2, 1 G2, 1 w3 G3, 1 Y2, 2 G2, 2 G4 2 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G2, 3 G3 2 Y2, 1 G1, 1 w2 G4, 1 Y2, 3 G1, 3 G3 2 Y2, 3 G4
2 2 ,3 2 2

w1 Y2, 2 G1, 2 G2
2

2

,2

w3 Y2, 1 G1
2 2

2

2 ,1

Y2, 2 G4
2, 3

w1 Y2, 3 G1
2 2, 1

,3

2 Y2, 2 G4

2 ,2

w3 Y2, 1 G1, 1 G
,2 ,1 ,1 ,1

2

2, 1

2 Y2, 2 G4

,2

w1 Y2, 3 G1, 3 G
2

2 Y2, 1 G
,2 ,2 ,1

w3 G4, 1 Y2, 2 G1, 2 G4
2

2 Y2, 1 G2, 1 w3 G3, 1 Y2, 2 G2, 2 G3 2 Y2, 3 G1, 3 w1 G4, 3 Y2, 2 G1, 2 G3 2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G2, 1 G3 2 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G3
2 2 2 ,1 2 2

2 Y2, 3 G1, 3 w2 G4, 3 Y2, 1 G1, 1 G3 2 Y2, 3 G1, 3 w2 G4, 3 Y2, 1 G1, 1 G4 2 Y2, 2 G2, 2 w3 G4, 2 Y2, 1 G1, 1 G4
1, 3 2 2

2 Y2, 3 G
,1 ,3

w1 G2, 3 Y2, 2 G3
2

2

2 ,2 ,3 ,3

2 Y2, 2 G1, 2 w3 G4, 2 Y2, 1 G1, 1 G3 2 Y2, 2 G1, 2 w1 G4, 2 Y2, 3 G1, 3 G4
2

2 Y2, 2 G1, 2 w1 G4, 2 Y2, 3 G1, 3 G3 2 Y2, 2 G1, 2 w1 G4, 2 Y2, 3 G2, 3 G4
2

111


2 Y2, 2 G1, 2 w3 G4, 2 Y2, 1 G1, 1 G4 Y2, 1 G2
2 ,1

2

,1 2 ,1

2 Y2, 1 G2
2

2 ,1

w2 Y2, 3 G3, 3 G4
2

2

,3

w2 Y2, 3 G3
2 2

2

2 ,3

Y2, 1 G2
4, 3 4, 1 2, 1

w3 Y2, 2 G3
2 2, 2 2

,2

2 Y2, 1 G2
2 3, 1

2 ,1

w3 Y2, 2 G3, 2 G
2 ,2

2

4, 2 2

2 Y2, 2 G2 2 Y2, 2 G2 2 Y2, 3 G3

,2 2 ,2 2 ,3

w1 Y2, 3 G3, 3 G w3 Y2, 1 G3, 1 G w2 Y2, 1 G1, 1 G
2 2 2

Y2, 2 G

w3 Y2, 1 G
2

Y2, 2 G2

w1 Y2, 3 G3

2

,3

2 Y2, 2 G 2 Y2, 3 G
,2

2, 2 1, 3

w1 G4, 2 Y2, 3 G1, 3 G4 w1 G3, 3 Y2, 2 G2, 2 G3
2 2

,3 ,2 ,1 2 ,2 2 ,3

2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G2, 2 G3
2

2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G1, 1 G4 ( Y2, 1 G1,
2 2 2 1 2 ,2 2 ,3 1 2 3 2 3 3

2 Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G1, 3 G4, 3 ) Y2, 3 G1 Y2, 1 G3
2 ,3 2 ,1

w2 w3 w1 w w1 w
2 2 3

2

2

Y2, 2 G1 Y2, 3 G2

w1 w
2

2

2 3 2 2

w1 w w2 w
2 2

2

2 2 2 3 2

Y2, 1 G2 Y2, 2 G3
3, 1 3, 2 3, 3 4, 1 3, 3 4, 2

2 ,1 2 ,2

w2 w3 w1 w3
1, 1 1, 2 1, 3 2, 2 2, 3 3, 3 2

2

Y2, 2 G2 Y2, 3 G3
2 2

2 2 2

w1 w2
2

2 Y2, 1 G1, 1 w2 w3 G2
2 2 ,2 ,3 ,1 ,2 ,1

,1

2 Y2, 1 G1, 1 w2 w3 G 2 Y2, 2 G1, 2 w1 w3 G 2 Y2, 3 G1, 3 w1 w2 G 2 Y2, 1 G2, 1 w2 w3 G 2 Y2, 3 G2, 3 w1 w2 G 2 Y2, 2 G3, 2 w1 w3 G Y2, 2 G4 Izx :=
2 ,2 2 2 2 2 2 2 2 2 2 2

2 Y2, 1 G 2 Y2, 2 G 2 Y2, 3 G 2 Y2, 2 G 2 Y2, 3 G 2 Y2, 3 G
2 ,3 3, 1 2

w2 w3 G4, w1 w3 G4, w1 w2 G4, w1 w3 G3, w1 w2 G4, w1 w2 G4,
2 2 2 2 2 2 2 2 2 2 2

2 Y2, 2 G1, 2 w1 w3 G2 2 Y2, 3 G1, 3 w1 w2 G2 2 Y2, 1 G2, 1 w2 w3 G3 2 Y2, 2 G2, 2 w1 w3 G4 2 Y2, 1 G3, 1 w2 w3 G4 Y2, 1 G4,
2 1 2 2 2 2 2 2 2 2

w2 w3

2

2

w1 w

2

2 3

Y2, 3 G4
2

w1 w2 ) Y2, 2 G
,2 3, 2

P Lx Lz ( Y2, 1 G2, 1 w3 G
2

G4,

2

Y2, 3 G3,
2

2 3

w1 Y2, 2 G2, 2 G4
,2

2

,2

2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G1, 2 G3 Y2, 3 G3, 3 w1 G4, 3 Y2, 2 G1, 2 G3
2 2 ,2

Y2, 3 G3, 3 w1 G4, 3 Y2, 2 G2, 2 G4 Y2, 3 G1, 3 w2 G3, 3 Y2, 1 G3, 1 G4
2 2 ,1

2 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G1, 2 G3 Y2, 2 G1, 2 w3 G3, 2 Y2, 1 G1, 1 G2 Y2, 1 G3, 1 w3 G4, 1 Y2, 2 G2, 2 G4
2 2 2 ,1 ,2

,2

Y2, 3 G1, 3 w2 G3, 3 Y2, 1 G1, 1 G2 Y2, 2 G2, 2 w3 G4, 2 Y2, 1 G1, 1 G2 Y2, 1 G3, 1 w3 G4, 1 Y2, 2 G1, 2 G3
2 2 2 ,1 ,2

,1

2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G2, 1 G4 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G1, 2 G3 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G2, 1 G4 Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G3, 3 G4
2 2 2 2 ,2 ,1 ,3

,1

Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G2, 2 G4 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G1, 1 G3 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G1, 2 G2 Y2, 1 G1, 1 w3 G3, 1 Y2, 2 G3, 2 G4
2 2 2 2 ,1 ,2 ,2

,2

2 Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G1, 3 G3 Y2, 2 G3, 2 w1 G4, 2 Y2, 3 G1, 3 G3
2 ,3

,3

Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G1, 2 G2 Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G3, 3 G4
2 ,3

,2

112


Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G3, 2 G4
2

2

,2 ,2

Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G1, 2 G2 Y2, 2 G3
2 ,2

2

,2

2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G2, 2 G4 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G1, 2 G4 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G2, 1 G3 Y2, 2 G2, 2 w1 G3, 2 Y2, 3 G4 Y2, 1 G2
2 ,1 2 2 2 ,3 4, 2 2 ,1 2 2 ,2 ,1

w1 Y2, 3 G1, 3 G4
2

2

,3 ,2

Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G2, 2 G3 Y2, 2 G2, 2 w3 G3, 2 Y2, 1 G4
2 3, 3 2 2, 3 2, 2 2 2 ,1

Y2, 3 G Y2, 3 G Y2, 2 G
,1 ,2

w2 Y2, 1 G w2 Y2, 1 G
2 2

2

2 1, 1 2 4, 1

Y2, 3 G3 Y2, 3 G2
2 ,3

2 ,3 2 ,3

w1 Y2, 2 G1 w1 Y2, 2 G4
2

2

2 ,2 2 ,2

w3 Y2, 2 G1, 2 G
2

Y2, 2 G2, 2 w3 G3, 2 Y2, 1 G1
2

w1 G3, 2 Y2, 3 G1
2

Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G1, 1 G4 Y2, 1 G1, 1 w3 G2, 1 Y2, 2 G1, 2 G4 Y2, 2 G2, 2 w1 G4, 2 Y2, 3 G1 Y2, 2 G1, 2 w1 G3, 2 Y2, 3 G4 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G3
2 2 2 2 2 ,3 2 ,3 2 ,1 2

Y2, 1 G1, 1 w2 G2, 1 Y2, 3 G1, 3 G4 Y2, 1 G1, 1 w3 G2, 1 Y2, 2 G2, 2 G3
1, 1 1, 2 2, 2 2

,3 ,2

Y2, 1 G Y2, 2 G Y2, 2 G
,3 ,2 ,2

w3 G3, 1 Y2, 2 G4 w3 G3, 2 Y2, 1 G4
2 2

2

2 ,2 2 ,1 ,1 2 ,2 ,1 ,1

w3 G3, 2 Y2, 1 G3, 1 G4
2

Y2, 1 G3, 1 w2 G4, 1 Y2, 3 G1, 3 G4 Y2, 1 G3, 1 w3 G4, 1 Y2, 2 G1, 2 G4 Y2, 3 G3, 3 w1 G4, 3 Y2, 2 G1, 2 G4 Y2, 1 G3
2 ,1 2 2

Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G3
2

Y2, 3 G3, 3 w2 G4, 3 Y2, 1 G1, 1 G4 Y2, 2 G1, 2 w3 G2, 2 Y2, 1 G1, 1 G4
1, 2 2

w3 Y2, 2 G2, 2 G
2

2

4, 2

Y2, 2 G
,2

w3 G2, 2 Y2, 1 G2, 1 G3
2

2

,1 ,3 2 2 2

Y2, 3 G3, 3 w1 G4, 3 Y2, 2 G2, 2 G3 Y2, 3 G2 Y2, 2 G3
2 ,3 2 ,2

Y2, 2 G1, 2 w1 G2, 2 Y2, 3 G2, 3 G3
2 2, 3 2

w2 Y2, 1 G1, 1 G w1 Y2, 3 G1
2 2 2 ,3

2

4, 1

Y2, 3 G
2 ,2

w1 Y2, 2 G

2

1, 2

G4,

2

Y2, 2 G3,
2 ,1

w3 Y2, 1 G1,

2

2 1 ,3

Y2, 2 G3
,3 ,3

w3 Y2, 1 G1, 1 G4
2 ,1 2

,1

Y2, 1 G2
,2

w2 Y2, 3 G1, 3 G3

Y2, 2 G1, 2 w1 G2, 2 Y2, 3 G1, 3 G4 Y2, 2 G3, 2 w1 G4, 2 Y2, 3 G1, 3 G4 Y2, 2 G2
2 ,2 2

Y2, 1 G2

w3 Y2, 2 G1, 2 G3
2

Y2, 2 G3, 2 w3 G4, 2 Y2, 1 G1, 1 G4
2 2, 2

,1

w1 Y2, 3 G1, 3 G
2

2

3, 3

Y2, 2 G
,1 ,3

w3 Y2, 1 G
2

2

1, 1

G3,

1 ,3

Y2, 3 G2, 3 w2 G3, 3 Y2, 1 G3, 1 G4 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G3, 3 G4 Y2, 3 G1, 3 w2 G3, 3 Y2, 1 G4 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G4 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G1 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G4
2 2 2 2 2 ,1 2 ,2 2 ,1 2 ,3 2

Y2, 2 G3, 2 w1 G4, 2 Y2, 3 G2, 3 G3 Y2, 3 G2
2 2, 3 2, 3 2, 1 2, 1 2 ,3 2

w1 Y2, 2 G1, 2 G3
1, 1

2

,2

Y2, 3 G Y2, 3 G Y2, 1 G Y2, 1 G

w2 Y2, 1 G
2

G3,
,2

1 2 2

w1 G4, 3 Y2, 2 G1 w3 G3, 1 Y2, 2 G4 w3 G3, 1 Y2, 2 G1
2 2

,2 2 ,2

113


Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G1 Y2, 1 G3
2 ,1 2

2

2 ,3 4, 3 2 ,2 4, 3 4, 1

Y2, 1 G Y2, 3 G Y2, 1 G Y2, 1 G Y2, 2 G

2 3, 1 2, 3 3, 1 2 2, 1 2 2, 2 2

w3 Y2, 2 G
2

2

1, 2

G4,
,1

2 2 2 ,1 2 ,1 2 2 ,3 2 ,3 2 1 2 2

w2 Y2, 3 G1, 3 G
2

w2 G3, 3 Y2, 1 G4
2

Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G4 Y2, 1 G2 Y2, 2 G2 Y2, 2 G2 Y2, 3 G3
2 ,1 2 ,2 2 ,2 2 ,3 2

w3 Y2, 2 G w3 Y2, 2 G w1 Y2, 3 G
3, 3 2 2 2

2

2 1, 2 2 4, 2 1, 3 1, 1 1, 2

Y2, 1 G3 Y2, 1 G2
3

w2 Y2, 3 G1 w2 Y2, 3 G4
2 2 2 2 2

w2 Y2, 3 G1, 3 G w3 Y2, 1 G1, 1 G w1 Y2, 3 G4
2 2 2 ,3 2

G4,
2

Y2, 2 G2,

w3 Y2, 1 G4,

Y2, 3 G2, 3 w2 G
4, 1

Y2, 1 G

Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G1,
2 ,2

w2 Y2, 1 G1, 1 G
2

Y2, 3 G
,1 ,2 ,1

2 3, 3

w1 Y2, 2 G
2

G4,

Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G3, 1 G4 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G1, 2 G2 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G1, 1 G2 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G1 Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G4 Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G1
2 2 2 2 2 ,2 2 ,3 2 ,3 2 2

2 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G2, 2 G4 Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G3, 2 G4 Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G3
2 1, 1 2, 2 2, 1 2 2 ,2 ,3 2 ,2

Y2, 1 G Y2, 2 G Y2, 1 G
,2

w2 G2, 1 Y2, 3 G2, 3 G3 w3 G4, 2 Y2, 1 G1 w2 G4, 1 Y2, 3 G3
2 2 2 2 ,1 2 ,3

Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G2, 2 G4 Y2, 3 G3, 3 w2 G4, 3 Y2, 1 G1 Y2, 3 G2 Y2, 2 G3
2 ,3 2 ,2 2 2 2 ,1 4, 1 2, 1 2 ,3

Y2, 3 G2, 3 w1 G3, 3 Y2, 2 G1, 2 G3
3, 3 2, 3 3, 2 2, 1

,2

Y2, 3 G Y2, 3 G Y2, 2 G Y2, 1 G
,3

w1 G4, 3 Y2, 2 G1
2

2

2 ,2 2 2 ,1 ,2 2 ,1 ,2

w2 Y2, 1 G3, 1 G w3 Y2, 1 G1, 1 G
2 2

w1 Y2, 2 G
2

2

3, 2

G4,

w3 G4, 2 Y2, 1 G1
2

Y2, 2 G3, 2 w1 G4, 2 Y2, 3 G1
2

w3 G4, 1 Y2, 2 G2, 2 G3
2

Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G2, 3 G3 Y2, 3 G3
2 ,3

Y2, 2 G1, 2 w3 G2, 2 Y2, 1 G3
2 1, 1 3, 1 4, 1 2, 1

w1 Y2, 2 G1, 2 G
2

2

2, 2 2

Y2, 1 G Y2, 1 G Y2, 1 G Y2, 1 G
,2 ,3 ,1

w3 G3, 1 Y2, 2 G2, 2 G3 w2 G4, 1 Y2, 3 G1
2 2 2 ,3 2

Y2, 1 G3, 1 w3 G4, 1 Y2, 2 G1 Y2, 3 G4 Y2, 1 G4
2 ,3 2 ,1 2

,2 2, 1 2, 3

w2 Y2, 1 G1, 1 G w2 Y2, 3 G1, 3 G
2 2

w3 Y2, 2 G
2

2

1, 2

G2,

w3 G3, 1 Y2, 2 G1, 2 G3
2

,2 ,3 ,3

Y2, 1 G2, 1 w3 G3, 1 Y2, 2 G2, 2 G4 Y2, 1 G2, 1 w2 G3, 1 Y2, 3 G1, 3 G3 Y2, 2 G1, 2 w3 G3, 2 Y2, 1 G1, 1 G4 Y2, 2 G4
2 ,2 2 2

Y2, 1 G1, 1 w2 G3, 1 Y2, 3 G2, 3 G3 Y2, 1 G1, 1 w2 G4, 1 Y2, 3 G1, 3 G3 Y2, 3 G4
2 4, 2 2 ,3 2 2

w1 Y2, 2 G1, 2 G2
1, 3

2

,2

w3 Y2, 1 G1, 1 G
2

2

2, 1

Y2, 2 G
,2

w1 Y2, 3 G
2

G2,

3 ,1

Y2, 1 G2, 1 w3 G4, 1 Y2, 2 G1, 2 G4

Y2, 3 G1, 3 w2 G4, 3 Y2, 1 G1, 1 G3

114


Y2, 3 G1, 3 w1 G4, 3 Y2, 2 G1, 2 G3 Y2, 2 G2, 2 w3 G4, 2 Y2, 1 G1, 1 G4 Y2, 3 G1, 3 w1 G2, 3 Y2, 2 G3
2 2 2 ,2 2

2

,2 ,1

Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G2, 1 G3 Y2, 3 G1, 3 w2 G2, 3 Y2, 1 G3
2 1, 2 2 2 ,1 ,1

2

,1

Y2, 2 G
,3

w3 G4, 2 Y2, 1 G1, 1 G3
2

Y2, 2 G1, 2 w1 G4, 2 Y2, 3 G1, 3 G3 Y2, 1 G2 Y2, 2 G2
2 ,1 2 ,2

Y2, 2 G1, 2 w1 G4, 2 Y2, 3 G2, 3 G4
2 2, 1 2 2, 2

,3

w2 Y2, 3 G3, 3 G w1 Y2, 3 G3, 3 G
2 2

2

4, 3 4, 3

Y2, 1 G Y2, 2 G
,3 ,2 ,1

w3 Y2, 2 G w3 Y2, 1 G
2 ,3 2 2

2

3, 2 3, 1

G4, G4,

2 1 ,1 ,2

Y2, 2 G2, 2 w1 G4, 2 Y2, 3 G1, 3 G4 Y2, 3 G1, 3 w1 G3, 3 Y2, 2 G2, 2 G3 Y2, 3 G2, 3 w2 G4, 3 Y2, 1 G1, 1 G4 Y2, 1 G1
2 ,1 2 2

Y2, 3 G3

w2 Y2, 1 G1, 1 G2
2

Y2, 3 G2, 3 w1 G4, 3 Y2, 2 G2, 2 G3
2

Y2, 1 G2, 1 w2 G4, 1 Y2, 3 G1, 3 G4, 3 )
2

( w2 w3
2 2 2 2 2

w2 w
2 ,2 2 ,3

2

2 3 2 3 2 2

Y2, 2 G1

2 ,2

w1 w3
2 ,3 2

2

Y2, 3 G1
2 2

2 ,3

w1 w
2 ,1

2

2 2 2 3

Y2, 1 G2

2 ,1

Y2, 2 G2 Y2, 3 G3

w1 w w1 w
2 2

2

Y2, 3 G2

w1 w2
2

Y2, 1 G3
,1 2

w2 w

2

Y2, 2 G3
2 2

,2 ,1

w1 w3

2 Y2, 1 G1, 1 w2 w3 G2
2 4, 1 4, 2 4, 3 3, 2 4, 3 4, 3

2 Y2, 1 G1, 1 w2 w3 G3
2 2 3 1 2 1

2 Y2, 1 G1, 1 w2 w3 G 2 Y2, 2 G1, 2 w1 w3 G 2 Y2, 3 G1, 3 w1 w2 G 2 Y2, 2 G2, 2 w1 w3 G 2 Y2, 3 G2, 3 w1 w2 G 2 Y2, 3 G3, 3 w1 w2 G > assign(mc):
2 2 2 2 2 2 2 2 2 2

2 Y2, 2 G 2 Y2, 3 G 2 Y2, 1 G 2 Y2, 2 G 2 Y2, 1 G Y2, 1 G rz

1, 2 1, 3 2, 1 2, 2 3, 1 2

w1 w3 G2, w1 w2 G2, w2 w3 G3, w1 w3 G4, w2 w3 G4,
2 2 2 2 2 2 2 2 2 2

2

2 Y2, 2 G1, 2 w1 w3 G3 2 Y2, 3 G1, 3 w1 w2 G3 2 Y2, 1 G2, 1 w2 w3 G4 2 Y2, 3 G2, 3 w1 w2 G3 2 Y2, 2 G3, 2 w1 w3 G4
2 ,2 2 2 2 2 2 2 2 2

2

,2 ,3 ,1 ,3 ,2 2

4, 1

w2 w3

Y2, 2 G4

w1 w3

2

2

Y2, 3 G

4, 3

w1 w2 )

2

2

rx

> pn:=5;B:=readdata(".\\Output.st",24+pn); pn := 5 B := [ [ 2.000000, .8700000, 3000.000, 44000.00, 400.0000, 14.04726, 16.39794, 27.35160, 14.04726, 16.39795, 27.35160, 14.04726, 16.39794, 27.35160, 14.04726, 16.39795, 27.35160, -.0006753510, .0006737952, -.0006730062, -.0006753488, -.0007192560, -.0006304674, -.0006753497, -.0006737934, .0006730067, -.0006753526, .0007192575, .0006304678 ] ] , , > i:=1: w:=; > Y:=Matrix(<,
115


,17+pn]|B[i,18+pn]>,,>); rx0:=B[i,1]; rz0:=B[i,2]; Ixx0:=B[i,3]; Izz0:=B[i,4]; Izx0:=B[i,5]; Mass0:=15000; w := [ 14.04726, 16.39794, 27.35160 ] -.0006753510 -.0006753488 Y := -.0006753497 -.0006753526 .0006737952 -.0007192560 -.0006737934 .0007192575 -.0006730062 -.0006304674 .0006730067 .0006304678

rx0 := 2.000000 rz0 := .8700000 Ixx0 := 3000.000 Izz0 := 44000.00 Izx0 := 400.0000 Mass0 := 15000 > G:=multiply(Y, inverse(multiply(transpose(Y),Y))); -370.1777506 -370.1788100 G := -370.1783343 -370.1776585 346.8408488 -395.6877211 -370.2442808 -370.6773229 -346.8421470 395.6870886 370.2428281 370.6766316

> P:=-2000; Lx:=4; Lz:= 1.74; P := -2000 Lx := 4 Lz := 1.74 > <<" ","rx","rz","Izz","Ixx","Izx","Mass">|<" ",rx 0,rz0,Izz0,Ixx0,Izx0,Mass0>|<" ",rx,rz,Izz,Ixx,Izx,Mass >|<" ",abs(rx0-rx),abs(rz0-rz),abs(Izz0-Izz),abs(Ixx0-Ix x),abs(Izx0-Izx),abs(Mass0-Mass)>|<" %",100*abs(rx0-rx)/rx,100*abs(rz0-rz)/rz,100*abs(Izz0-Izz)/Iz z,100*abs(Ixx0-Ixx)/Ixx,100*abs(Izx0-Izx)/(Izx+1),100*abs(Mas s0-Mass)/Mass>>; " "rx" "rz" "Izz" "Ixx" "Izx" "Mass" > " " " 2.000000 .8700000 44000.00 3000.000 400.0000 15000 " " 2.000005143 .8699994285 44031.41489 3005.936169 399.9900534 15007.77822 " " -5 .5143 10 -6 .5715 10 31.41489 5.936169 .0099466 7.77822 %" " .0002571493387 .00006568969832 .07134653764 .1974815387 .002480510406 .05182792473

116