CC BY
65
УДК 681.3.06, 531.36
А.В. БАНЩИКОВ, Л.А. БУРЛАКОВА, В.Д. ИРТЕГОВ, Т.Н. ТИТОРЕНКО ЗАДАЧИ МЕХАНИКИ И КОМПЬЮТЕРНАЯ АЛГЕБРА_____________________
Abstract: The paper suggests a description of the functional content of software complexes stability, LinModel, etc., which are intended for modeling and qualitative investigation of both systems of interconnected bodies and electromechanical systems. The complexes have been elaborated on the basis of the computer algebra system MATHEMATICA. Examples of application of this software are given. In particular, investigation of the invariant manifold of steady motions in the Kirchhoff’s problem, which is related about the motion of a body in ideal fluid in a special case, has been fulfilled.
Key words: functional content, invariant manifold, computer algebra of the body.
Анотація: Стаття містить опис функціонального наповнення комплексів програм стійкість, LinModel та ін., призначених для моделювання та якісного дослідження систем взаємозалежних тіл і електромеханічних систем. Комплекси створені на базі системи комп'ютерної алгебри MATHEMATCA. Наведені приклади використання цього програмного забезпечення. Зокрема, виконане дослідження інваріантного різноманіття стаціонарних рухів у задачі Кірхгофа про рух тіла в ідеальній рідині в спеціальному випадку.
Ключові слова: функціональне наповнення, інваріантне різноманіття, комп’ютерна алгебра, рух тіла.
Аннотация: Статья содержит описание функционального наполнения комплексов программ
устойчивость, LinModel и др., предназначенных для моделирования и качественного исследования систем взаимосвязанных тел и электромеханических систем. Комплексы созданы на базе системы компьютерной алгебры MATHEMATICA. Приведены примеры использования этого программного обеспечения. В частности, выполнено исследование инвариантного многообразия стационарных движений в задаче Кирхгофа о движении тела в идеальной жидкости в специальном случае.
Ключевые слова: функциональное наполнение, инвариантное многообразие, компьютерная алгебра, движение тела.
1. Вступление
Системы компьютерной алгебры (СКА) [1] уже давно широко используются для решения сложных задач в различных отраслях знаний, в том числе в механике. Регулярно пополняемый обзор [2] дает достаточно полное представление о развитии этого направления в нашей стране. Безусловно, символьные вычисления требуют достаточно больших ресурсов вычислительных машин и для задач большой размерности приводят к очень громоздким выражениям. Однако многолетний опыт показывает, что использование СКА оправдано на этапах получения математических моделей (систем дифференциальных уравнений движения) и в случаях, когда необходимо проводить качественный анализ систем, в частности, параметрический.
В лаборатории "Устойчивости движения" Иркутского вычислительного центра (ныне Институт динамики систем и теории управления СО РАН) с 1976 года ведутся исследования в области механики и устойчивости движения с использованием САВ. Разработаны алгоритмы и созданы эффективные проблемно-ориентированные системы аналитических вычислений для ЭВМ типа «Эльбрус» на языке АЛГОЛ-ГДР и для ЕС ЭВМ на языке PL-1, которые были использованы в качестве инструментальных средств для пакетов прикладных программ "Динамика" и "Механик" [3, 4]. С помощью этих пакетов решены некоторые прикладные задачи.
В основе алгоритмов, которые реализованы в этих пакетах, лежат классические методы аналитической механики и теории устойчивости движения. Указанное программное обеспечение позволяет по геометрическому описанию системы твердых и деформируемых тел составлять их конечномерные модели в виде уравнений Лагранжа второго рода и Эйлера-Лагранжа, а также получать линейные уравнения в окрестности установившегося движения. Эти пакеты отличаются от других программных разработок тем, что они содержат крупные блоки, позволяющие проводить
82 © Банщиков А.В., Бурлакова Л.А., Иртегов В.Д., Титоренко Т.Н., 2008
ISSN 1028-9763. Математичні машини і системи, 2008, № 4
в символьном виде качественные исследования полученных дифференциальных уравнений: строить первые интегралы, линейные и квадратичные по обобщенным скоростям (квазискоростям); для линейных (циклических) и квадратичных по скоростям интегралов выделять стационарные движения и инвариантные многообразия стационарных движений и исследовать их устойчивость.
В девяностых годах прошлого века возникла необходимость в переходе на персональные компьютеры. В настоящее время широко используются постоянно совершенствующиеся универсальные системы компьютерной алгебры, например, Maple, Mathematica и др., не только для решения конкретных задач, но и для создания проблемно-ориентированных комплексов программ в различных областях науки, в том числе и в механике. Элементы символьных преобразований становятся также неотъемлемой частью многих вычислительных систем, например, MatLab и др.
Для IBM PC авторами разработаны на базе универсальной системы Mathematica программный комплекс LinModel [5], информационно-исследовательская система "Устойчивость" [6] и комплекс программ по моделированию и качественному анализу электромеханических систем и электрических цепей [7]. Все эти комплексы содержат модули, обеспечивающие построение дифференциальных уравнений движения, и модули качественного исследования, которые можно использовать автономно. Различие - в специфике объектов исследования (система твердых тел, система твердых и деформируемых тел, электромеханическая система) и в функциональных возможностях блоков качественного анализа.
Дадим краткое описание функционального наполнения упомянутых разработок.
2. Описание моделей
2.1. Модель механической системы
Пакеты позволяют автоматизировать получение математических моделей (дифференциальных уравнений) таких механических систем, которые допускают представление в виде системы взаимосвязанных твердых и деформируемых тел S. (i = 1,____,N). Следуя [8], рассматриваем
деформируемое тело по следующей приближенной модели. При отсутствии деформации определяется твердый "каркас" - отсчетная конфигурация. При деформации точки тела совершают малые движения в окрестности тех положений, которые они занимали бы в абсолютно твердом теле. С каркасом в некоторой его точке O жестко свяжем систему координат S . Положение точки
M'(x1, x2, x3) каркаса определяется в S радиусом - вектором р = OM' , а смещение точки M от M при деформации определяется вектором u (при отсутствии деформации u = 0 точка M совпадает с M ); u зависит от положения точки в недеформированном состоянии и от времени:
где qo.it) - обобщенные координаты деформации, иа(хх, х2, х3) - собственные формы колебаний или их функции, п! - некоторое конечное число. Тела в системе могут быть связаны
так, что существуют общие точки 0г е и 01 е (у = 1,_____,Ы, г = 2,...,Ы) , или могут совершать
относительные поступательные или винтовые движения друг относительно друга (рис. 1).
(2.1.1)
Рис. 1. Модель механической системы
Тело 5 является носителем, его положение определяется относительно инерциальной
системы координат Е0. С твердым телом 5 свяжем систему координат Е1, имеющую начало в
точке 01 е 5, с телом Б. (твердым или "каркасом" деформируемого) в точке О. е Б. свяжем
систему координат X, поворот X относительно X определяется матрицей а1, элементы которой являются функциями обобщенных координат поворота. В точке О е Б (О е S¡) введем
систему координат X, неподвижную относительно тела 5. в малой окрестности точки 01. При
отсутствии деформации тела Б. Е1 и X совпадают. Обобщенные координаты для тела могут
выбираться так, что положение X задается по отношению к Е1, при этом определяется матрица
2.1.1. Кинетическая энергия системы
Кинетическая энергия такой сложной системы может быть вычислена как сумма кинетических энергий Т каждого тела 51:
ю' - угловая скорость поворота X1 относительно X за счет деформации; эта угловая скорость
учитывается, если движение X определяется по отношению Е1:
(2.1.2)
Здесь М] — масса тела, V0 — абсолютная скорость точки О.:
(2.1.3)
ю. — абсолютная угловая скорость тела:
ю = ю + ю1 + к ю
(2.1.4)
(2.1.5)
юі - угловая скорость относительного движения X' по отношению к Xі (кі = 0) или Xі (кі = 1) ; у! - относительная скорость точки О., если тело Бі не имеет общих точек с :
■О = Роі+кіі0=,ио;?0"; (2.1.6)
г, — радиус-вектор центра масс в теле:
Сі і (* '
■Сі = РС+к, (20=, ім,иоі0т ум, = Рсі+к , 20=, ио; «о; чл .7)
0° — тензор инерции тела относительно точки О, :
&О; = 0°+2к, 2І, л° +к,ХІЖь € ?; (2.1.8)
00і — тензор инерции твердого (отвердевшего) тела; матрицы Л ° и определяются формами и00 (і = 1,2,3) из (2.1.1) [4]:
В ; = Г [ р х и а'■ ] йш + У Ь [ иь х и а. ] йш ; 4. = Г (и а'■ • ) йш ;
¿М. р ¿М. р ¿М.
к. = 0 , если тело твердое; к{ = 1, если тело деформируемое.
2.1.2. Силовая функция системы
Программное обеспечение позволяет вычислять силовые функции некоторых полей. Силовая функция для тела 5. (/ = 1,...,М), находящегося в ньютоновском поле тяготения, с точностью до второго порядка малости относительно величин, равных отношению характерных размеров тела к расстоянию точки 01 до неподвижного притягивающего центра О' (в котором помещена
неподвижная система координат X0), вычисляется по формуле
=пм—^[ _ ) + — ^[ о|У;+(Г0. г4) ]+пМ.[ ^+
г г3 1 ' ^ ' г3 2 ' ' г5 2
01 01 °1 °1
+(rO ■ rO>(r“ -^)]+п(®0; +©O; +®0з)-b)TQO’b ■ (2.1.9)
Здесь п- постоянная тяготения; r° =
о - расстояние от точки O1 до притягивающего центра; Ь = auЬ, Ь = a01g, Ь — матрица 3ХІ косинусов углов, образуемых вектором r° с осями X'; g — матрица 3 Х1 косинусов углов, образуемых вектором r° с осями X0 ; —
элементы тензора инерции 0° (1.1.8); X' предполагаются главными осями деформации.
Для тела S' в поле постоянной тяжести U' = M . (g • r4 ) + const, здесь g — вектор
ускорения постоянной силы тяжести; r° = r,° + rj + r' . Если тело S деформируемое, то силовая
Ч °1 Oi Чі '
функция упругих сил представляется в виде квадратичной формы обобщенных координат
2.2. Описание электрической цепи
Электромеханические аналогии [7, 9] позволяют использовать методы аналитической механики для описания и изучения электрических цепей.
Рассмотрим линейную электрическую цепь, в которой произвольным образом соединены элементы сопротивления (Р), катушки индуктивности (Ь), конденсаторы (С), источники тока (I) и напряжения (Е). Независимыми переменными, описывающими состояние цепи, могут быть либо токи в участках цепи, либо узловые напряжения, либо токи и узловые напряжения. Выбор того или иного множества переменных влияет на структуру уравнений.
Если в качестве переменных состояния выбрать множество токов, то оно определяется по методу контурных токов. Общая идея метода контурных токов состоит в разбиении электрической цепи на независимые геометрические контуры и задании в каждом контуре контурного тока. С задачей построения независимых контуров тесно связана задача нахождения множества фундаментальных циклов графа. Каждому фундаментальному циклу в графе, соответствующему заданной электрической цепи, можно сопоставить уравнение, определяющее закон Кирхгофа для напряжений. Ни одно из этих уравнений не будет зависеть от остальных. Процедура построения фундаментальных циклов графа цепи проводится по известному в теории графов алгоритму. Найдем значение тока в каждой ветви как сумму или разность соответствующих контурных токов, обозначая через ^ ток в 1-ом контуре и считая его положительным в ребре {пк,пш} , если
пк < пш , и отрицательным в противном случае. Согласно [9] кинетическая энергия Т, силовая функция и , функция Релея Я линейной электрической цепи имеют вид
где Ь,, — индуктивность, С', — электроемкость, К, — омическое сопротивление, ?' —
заряд.
Для построения функции состояния нелинейной электрической цепи используется алгоритм [7, 10].
3. Уравнения движения
После вычисления кинетической энергии и силовой функции находим характеристическую функцию системы и строим уравнения движения.
3.1. Уравнения Лагранжа
Пусть характеристической функцией является лагранжиан Ь = Т + и , который представлен в форме:
деформации: и*е = —1/2^?а') Оір- Силовая функция всей системы тел: и = ^ ^и ‘'
і = Т+т + т+и = 2 У *=, У*= л„ (я, () я + У *= а (я, і) + лдя, І),
(3.1.1)
где к — число степеней свободы системы, q¡ — обобщенные координаты, Т0(д,¿) - часть кинетической энергии, не зависящая от обобщенных скоростей;
А(Я,О = T0(q,¿) + и(q) ; Л = Л^. . Тогда уравнения движения строятся как уравнения Лагранжа второго рода:
У * л д+ - У * У *
=1 п 1г 2^^ г=1 •*=1
ЭЛ ЭЛ дЛ
Л
Эя. Эяг Эя
,=1
і У
V Эя- Эя
і У
ХП* Эл,.,- . ЭЛ. эа .
+У —-<7+—L-----------0 = О , (ї = 1, —, * )•
ді у ді дді ‘
0і — обобщенные силы непотенциальной природы.
(3.1.2)
3.2. Уравнения Эйлера-Лагранжа
Пусть конфигурация системы задается независимыми параметрами я , як и квазискоростями -
линейными формами по обобщенным скоростям я. (в общем случае неоднородными) [8]:
п, = Ум (я) Я , ( 5 = 1>---,* +1)’ (321)
здесь п*+1 = Я*+1 = 1 (а*+1.(я) = <5*+и — символ Кронекера, Я*+1 = і ) . Предполагается, что матрица 5 = ||йг (я)|| невырожденная. Линейные формы дифференциалов обобщенных
координат йж8 = У*=+15(я)^, (5 = я,• ••, * +1) определяют дифференциалы квазикоординат. Кинетическая энергия системы Т в квазискоростях в общем случае имеет вид
т = 2 У *+1 У *=11 ЛК?) п, п + У *=1 л» п, + л*(я) = т + Т+Т
(3.2.2)
Уравнения Эйлера-Лагранжа, когда кинетическая энергия не содержит время I, имеют вид
У * л* п . + У *+1 У * У * у Л . п . п + У *+1 У * ул* п +
=1 , і т=1 і=1 г=1 т, п і т ^т=^ ^ у=1 т, г т
+У * У * (Элт — 1ЭЛт) п п + У * (ЭЛ^ — ЭЛ-) п —ЭЛ — — = О
¿-1т=Шу=1Ч^_ О З — ' у т ¿—¡у=1Ч^_ "Л у ’
Эр 2 Эр
"Эр Эр
Эр Эр
(3.2.3)
я = У*+1 Ь п , (, = 1,*),
1 , ^^^у = 1 ,У у“^ “ 5/5
где Q* =z := brs Qr — обобщенная сила непотенциальной природы по квазикоординате
Z: ,—,: да да
Z (—-----------—) b b , ( m, s = 1,_,:; и = 1, ,: +1) - трехиндексные
r=1 i=1 dq dq. 555 55 / i m
символы Больцмана-Гамеля ( b. — элементы матрицы b = ail).
Автоматизирован вывод уравнений Эйлера-Лагранжа и при наличии неголономных связей, а также избыточных обобщенных координат [4].
Для построения уравнений Гамильтона [8] выполняется переход к гамильтоновым переменным q., p. (p. = dL(q,q)/dqi); для построения уравнений Рауса - переход к переменным
Рауса [9]. Для электрических цепей характерно: det||д2L(q,q)/Э<7;Э<7.|| = 0,(/, j = 1,_,к) .
Поэтому не всегда удается выполнить стандартное преобразование обобщенных скоростей к гамильтоновым переменным. В этом случае осуществляется частичное преобразование [11].
4. Качественные исследования. Основные формулы и алгоритмы
4.1. Исследование по первому приближению
Полученные дифференциальные уравнения после разложения в ряд в окрестности стационарного движения могут быть представлены в виде
M x + (D + 2G) x + (K + P) x = Q( x, x), (4.1.1)
где x, x, H — матрица-столбец отклонений обобщенных координат, скоростей, ускорений от их значений в невозмущенном движении, Q(x, x) — нелинейные по x, x, слагаемые выше первого порядка малости, Q(0,0) = 0 , M = MT , D = DT - матрица диссипативных сил, K = KT -
матрица потенциальных сил, G = —GT - матрица гироскопических сил, P = —PT - матрица неконсервативных сил. Если принять Q(x, x) = 0 , то получим уравнения первого приближения. Формируем характеристическое уравнение для уравнений первого приближения:
det(M m2 + (D + 2G) m + (K + P)) =
(4.1.2)
hQm + TrD m2:—1 + h2m —1 + _ + К:—1m + det(K + P) = 0.
Используя классические теоремы Ляпунова, по уравнению (4.1.2) можно провести исследование устойчивости по первому приближению. Для этого реализованы критерий Рауса-Гурвица отсутствия корней с положительными вещественными частями и критерий, в соответствии с которым все корни полинома (4.1.2) отрицательные вещественные. Анализ полученных неравенств выполняется Mathematica-программой “Reduce”. По уравнению (4.1.1) исследуется задача о влиянии структуры сил на устойчивость тривиального решения системы (4.1.1) [12, 13]. Если выполнить в уравнении (4.1.2) подстановку m = n + X , где X— заданное смещение, то можно решать задачи оценки заданной степени устойчивости для системы (4.1.1).
4.2. Построение первых интегралов дифференциальных уравнений движения
Алгоритмы построения первых интегралов для описанных систем базируются на инвариантных соотношениях (дифференциальных следствиях) уравнений движения. Для уравнений движения (3.2.3) некоторые из дифференциальных следствий приведены в работе [14]. Одно из них:
d_
dt
f
a t+,
=1 r, к+1
ЭТ
\
dW
-f (T+и)-Уk=iS' (K,+i-Wr) -
■I к
dW
a
- + W
(4.2.1)
г
> ч -/
При определенных условиях из (4.2.1) следует первый интеграл, квадратичный относительно квазискоростей. Например, если лагранжиан (3.1.1) и кинематические соотношения (3.2.1) не содержат явно время и непотенциальные силы О* отсутствуют, то имеет место обобщенный интеграл Якоби:
Т,(П, q) - T(q) - и (q) -1
í
ai
ЭТ ЭТ
\
+ -
= const.
(4.2.2)
Г э a
* эа
+I—(ь а + ь ,+1)
"ч / j \ sm m s, к+1 /
dq¡ s=i dq¡
ЭТ ЭТ
+-------= 0, и Эи/Эд1 = 0 для q. (i = 1,...,l;
Эаr Эq¡
l < k), то система (3.2.3) имеет циклические интегралы ^ as. дТ/dWs = const.
Пусть для несущего тела системы обобщенные координаты q1 , q2 , q3 описывают вращательное движение, q4 , q5 , q6 - перемещения точки O1, координаты q. (j = 7,...,k) -относительные перемещения несомых тел. Выберем квазискорости для несущего тела равными проекциям абсолютной угловой скорости w и линейной скорости v точки OJ на связанные с телом оси. Пусть квазискорости для несомых тел равны проекциям относительных скоростей или обобщенным скоростям. При таких условиях для уравнений (3.2.3) имеют место дифференциальные следствия:
d_
dt
d_
dt
ЭТ'1 ЭТ
Эу ) Эу
дТ Y ЭТ
Эю ) Эу
=2 Г ЭТ ь; d
I Эу ) dt
V
Л=( эт Tf+Г ЭТ Тм
ЭТ) ЭТ
Эю) Эю
Г„ЭТ^ ЭТ 0^Т ЭТ + 1 v— I — = 2M —; I Эу ) Эю Эю
(4.2.3)
)
Эю
Эу
где M = col (P, P2, P), F = col (p P, p) квазикоординатам ж ,
матрицы обобщенных сил по
Г 0 -V3 V2 ']
v = V3 0 -V1
V-V2 V1 0 )
ЭТ ;/ЭТ ЭТ ЭТЧ ЭТ ;/ЭТ ЭТ ЭТ
, — = col (—,------,—), — = col (--------,----,----
Эу dv, dv dv Эю dw dw dw
).
r=1
Из соотношений (4.2.3) при соответствующих условиях можно получить первые интегралы дифференциальных уравнений движения. Например, если (дТ/dv)TF = 0 , то из (4.2.3) следует (дТ/dv)T дТ/dv = const.
Если система описана уравнениями Лагранжа, то можно пересчитать (4.2.1), (4.2.3) в переменных q, q, учитывая (3.2.1), и получить дифференциальные следствия для уравнений
(3.1.2).
Другие возможности существования первых интегралов для систем (3.1.2) и (3.2.3) приведены в [14]. При проверке условий существования интегралов на компьютере повышенные требования предъявляются к возможностям СКА в тригонометрических преобразованиях.
Для уравнений вида (4.1.1) алгоритмы построения дифференциальных следствий и инвариантных соотношений описаны в работах [12, 13]. Для линейных систем (4.1.1) задача построения первых интегралов квадратичного вида сводится к решению систем алгебраических уравнений:
B + BT - (D - 2G)M-1N - NM-1 (D + 2G) = 0;
(K - P)M-1N - NM-1 (K + P) + BTM-1 (D + 2G) - (D - 2G)M-lB = 0;
(K - P) M-1B + BTM-1 (K + P) = 0;
2L = (K - P) M-1N + NM-1 (K + P) + BTM-1 (D + 2G) + (D - 2G) M-lB
(4.2.4)
относительно матриц квадратичной формы:
V = xTN x + xTLx + xTBTx + xTBx = const, L = LT, N = NT.
Система (4.2.4) является системой однородных линейных уравнений и может иметь или единственное тривиальное решение или множество решений. Если существует нетривиальное решение N*,B*, то решениями являются N* = N* ±N*T, Bj* = B* ±B*T и
N**= CN*C1T + CjN*CT, B**= CTB*C1 + C1TB*C, где C, Cj - матрицы, перестановочные с
M-(D + 2G), M~\K + P). Если в (4.1.1) Q(x,x) Ф 0, то интеграл ищем в виде
V = xTNx + xTLx + xTBTx + xTBx + F (x, x), и тогда к (4.2.4) нужно добавить уравнения
— + 2NM-Q - (D - 2G)M-1 — = 0; QTM- — = 0; 2BTM lQ - (K - P)M-1 dF = 0. dx dx dx dx
4.3. Исследование устойчивости вторым методом Ляпунова
Инвариантные соотношения, первые интегралы, получаемые из них, и другие инварианты можно использовать для нахождения стационарных движений (инвариантных многообразий стационарных движений - ИМСД) и построения функций Ляпунова при исследовании устойчивости последних [12]. Для этого используются алгоритмы теоремы Рауса-Ляпунова [15] и ее обобщений [16]. Как известно, эти алгоритмы сводятся к некоторой задаче на экстремум.
Пусть для системы обыкновенных дифференциальных уравнений
x = X(x), x е Rk (4.3.1)
известно некоторое число первых интегралов V0( x) = c0 , V1( x) = Cj , ..., Vm ( x) = cm .
Выбирая эти интегралы за базис, определяем алгебру первых интегралов системы над R , например, элементами этой алгебры будут интегралы вида
KJ = ¿ЛД (x)> J = {io,•••,h}î (0,1,.,m), l = c°nst. (4.3.2)
j=i
Для полноты анализа множества стационарности системы нужно рассматривать и нелинейные комбинации базовых интегралов. Стационарными назовем те решения и инвариантные многообразия системы (4.3.1), которые доставляют стационарное значение одному из первых интегралов алгебры первых интегралов задачи. Будем находить стационарные движения (ИМСД) из уравнений
dKj/dx = f (xj,...x ,1,...,\) = °, (i = 1,. , k), V1(x) = ch, ■■■, Vn (x) = cn. (4.3.3)
В регулярном случае, когда работает теорема Рауса-Ляпунова, достаточные условия устойчивости этих движений можно найти как условия знакоопределенности второй вариации d2KJ при нулевых значениях первых вариаций остальных первых интегралов. Эта задача имеет
хороший алгоритм решения для случая, когда все интегралы квадратичны. При вырожденности системы (4.3.3) (это характерно для электромеханических систем) или при наличии первых интегралов другой структуры трудности исследования возрастают, т.к. в этом случае система
(4.3.3) становится нелинейной. Для анализа таких систем современные СКА предоставляют возможность использовать метод базисов Грёбнера и другие методы [17]. Например, СКА Maple включает в себя пакет программ “PolynomialIdeals'' для работы с полиномиальными идеалами. С помощью программ этого пакета можно проанализировать множество решений полиномиальной системы уравнений, получить ответ на вопрос, конечное или бесконечное число решений имеет исследуемая система, найти размерность множества ее решений, если число решений бесконечно, и т. п .
Если в системе нет первых интегралов, для построения функций Ляпунова при исследовании устойчивости вторым методом Ляпунова можно использовать дифференциальные следствия уравнений движения (4.2.1), (4.2.3) и др. [14]. В литературе чаще всего используются дифференциальные следствия, определяющие закон изменения энергии. Но применение других дифференциальных следствий иногда приводит к смягчению условий устойчивости.
Для линейных уравнений (4.1.1) один из алгоритмов построения функции Ляпунова описан в [6] и состоит в следующем: умножаем слева (4.1.1) на матрицы вида xTL1, xTL2, xTL3 и
симметризуем полученные выражения, выделяя полную производную по времени. Матрицы L. равны матрицам исходного уравнения или их комбинациям [6].
5. Примеры
С помощью описанных пакетов в символьном виде исследованы различные конкретные механические системы. В частности, получены нелинейные уравнения движения ряда
робототехнических систем (робот-манипулятор шестизвенник, шагающая платформа на четырех ногах с двумя степенями свободы каждая); рассмотрены вопросы динамики и устойчивости стационарных движений искусственных спутников (гироскопическая рама в ньютоновском поле сил [18], спутник с гравитационной системой стабилизации на круговой орбите при различных законах управления гравитационным стабилизатором [12], спутник с солнечным парусом, спутник с гиродинами [19] и др.); исследованы ИМСД для уравнений, описывающих движение твердого тела с неподвижной точкой [20]; тела в жидкости [21] и др.
5.1. Пример 1
В качестве конкретного объекта исследования выбрана механическая система из 18 абсолютно твёрдых тел с 32 степенями свободы, находящаяся в ньютоновском центральном поле сил. Рассмотрим для этой системы тел ограниченную задачу, когда точкаО1 носителя движется по
круговой орбите радиуса Я (на рис. 1 | гО01 = Я) с постоянной скоростью. Введем в
точкеО1 орбитальную систему координатX1 так, что осьО1 X направлена по радиусу орбиты от центра, ось О1X по касательной к орбите в направлении вектора линейной скорости, ось О1 х^ по нормали к плоскости орбиты. Считаем орбитальную систему телом номер 1 (на рис. 1 обозначено как ^). Положение несущего тела $2 (спутника) определяется в X1. Центр масс спутника
находится в точке О1. Пренебрежем изменением положения центра масс системы при движении по орбите.
Для осуществления ориентации спутника в орбитальной системе установлены три гиродина (рамка и ротор каждого из них обозначены, соответственно, $4, $5; £6, $7; $8, $9). К спутнику
присоединены также антенна £3 (две степени свободы), две платформы с двигателями £10, £П (две
степени свободы у каждой) и тело $12 (три степени свободы), несущие тела £13-$18 , образующие
две «гантели» из трех тел каждая (приближенная модель солнечных панелей).
В соответствии с конкретной конфигурацией механической системы её геометрическое и кинематическое описание в программном комплексе ИпМоЬе! происходит в диалоговом режиме, где в определенных окнах программы для каждого тела пользователь задаёт входные данные. В качестве примера представим входные данные первых двух тел системы.
Тело 1 - орбитальная система координатX1. Масса и тензор инерции тела равны нулю. ПоложениеX1 описывается в инерциальной системе координатX0 с центром в притягивающей точке О'. Номера осей вращения {3, 0, 0} и углы поворотов {с, 0, 0} , т.е. положение орбитальной системы координат относительно инерциальной системы координат определяется одним поворотом на угол с вокруг третьей осиОх3° . Нули среди номеров осей вращения указывают на
то, что повороты вокруг этих осей отсутствуют. Радиус-вектор точки О1 в системе координат X0:
r<0 = (R cosc, RsinX, 0) ■ Радиус-вектор центра масс тела вS1: =(0,0,0) . Вектор
абсолютной линейной скорости точки O1: v0 = (-Rw0 sinX, Rw0 cosX, 0) , где w0 = X = const
- угловая скорость орбитальной системы координат.
Тело 2 (спутник) движется относительно тела 1. Масса тела - М2. Номера осей
вращения {1,3,2} и соответствующие углы поворотов {у, в, р}, т.е. положение спутника (S2)
относительно орбитальной системы координат S1 определяется последовательностью трех поворотов: вокруг первой оси на угол у; вокруг третьей оси на угол в; вокруг второй оси на угол
р . Радиус-вектор точки O2 соединения тел S1 и S2 : r¿2 =(0,0,0), т.е. точки O1 и O2 совпадают. Радиус-вектор центра масс: r1 = (0, 0, 0) , т.е. центр масс спутника находится в точке
O . Вектор относительной линейной скорости точки O2: vO = (0, 0, 0) . Тензор инерции тела в
точке O2 диагональный, т.е. Jx, Jy, Jz - главные моменты инерции спутника.
С помощью LinModel в символьном виде на PC для рассматриваемой механической системы решены следующие задачи: 1) для каждого тела построены матрица направляющих
косинусов aJ и вектор относительной угловой скорости
(J < i; i = 1,18; J = {0,1,2,4,6,8,12,13,16}) по заданной пользователем последовательности поворотов; по формулам (2.1.3) и (2.1.4) подсчитаны абсолютные линейные скорости точек Q. и
абсолютные угловые скорости всех тел; 2) получены кинетическая энергия системы тел (2.1.2) как квадратичная форма по обобщенным скоростям и силовая функция (2.1.9) приближенного ньютоновского поля тяготения; 3) построены 32 нелинейных уравнения движения (3.1.2) в форме Лагранжа 2-го рода; 4) получены условия существования заданного невозмущённого движения; 5) построены уравнения возмущенного движения в первом приближении (3.1.1); 6) из коэффициентов уравнений возмущенного движения сформированы M, G, K - квадратные матрицы 32 порядка характеристического уравнения системы (4.1.2). Для решения этих задач минимальные требования к оперативной памяти РС - 250 Мбайт. Время вычислений зависит от возможностей процессора. Например, на PC с двухъядерным процессором Core 2 Duo 2,4GHz было затрачено 28 минут процессорного времени. Построенная кинетическая энергия в текстовом формате занимает 18 Мбайт дисковой памяти.
Для подсистемы тел S1 -S9 (спутник с гиродинами) с помощью программного комплекса
LinModel исследована устойчивость по первому приближению двух найденных стационарных движений, одно из которых оказалось неустойчивым. Для другого стационарного движения найдены условия вековой и гироскопической устойчивости. Области гироскопической стабилизации имеют графическое 2D-представление■
5.2. Задача Кирхгофа о движении твердого тела в идеальной жидкости
Рассмотрим уравнения Кирхгофа для случая [22], когда они записываются в виде
5 =аг 52-ос\ Г + (( - э2)(Ьг2 - 53), Г = г2(огх + (г2 + 253) - гъ52,
= (о2 + Ь2)Г1 Г3 - (ог +ЬГ2)5! + (ог, -5^, т2 = Г,-Г1(ОГ1 + ( + 25,), (5.2.1)
53 = ((Г1 ОГ2 )5з , Г3 = Г 52 - Г2 51.
и обладают первыми интегралами:
2К0 = (512 + 522 + 2532) + 2(ОГ1 + (г2)53 - (о2 + Ь) Г = 2к, У1 = 51 Г1 + 52 Г2 + 5, Г, = С1,
V = Г12 + Г22 + Г32 = С2, 2И, = (г1 51 + Г2 52)((02 + Ь2)(Г 51 + Г2 52) + 2(051 + (52 )53) + (5.2.2)
+532 (512 + 522 + (оГ1 +РГ2 + 53)2) = 2с,.
Здесь 5 = (51,52,53), г = (г1,Г2,Г3) - векторы “импульсивного момента'' и “импульсивной силы'',
о, (- некоторые постоянные.
Для уравнений (5.2.1) находим инвариантные многообразия и исследуем их на устойчивость по алгоритму, описанному в п. 4.3, используя линейную комбинацию интегралов
(5.2.2):
К Л-\V\-lyj2.-1^/2.
Запишем условия (4.3.3) стационарности К по всем переменным:
151 -Дг1 - (о2 + (2)1 Г1(Г1 51 + г2 52) -Я3 (ог2 +(Г1) 52 53 -Я3(2оГ1 + 53)5153 = 0, Л0 52 -Дг2 - (о2 + (2)1 Г2(Г151 + Г2 52)-Я3 (оГ2 + (Г1) 51 53 -Л3(2(Г1 + 53 )52 53 = 0,
(5.2.3)
Л(оГ1 +(Г2 + 253)-Л1Г3 -Л3 (о51 + МХГ1 51 + Г2 52)-Л3(оГ1 +(Г2)253 - Л ^ + 52 ) 53 -
- 3Л(оГ + ( г2 )52 - 2^533 = 0, о Л0 53 - Л 51 - Л2Г1 - (о2 + (2)Л3 51 (Г151 + Г2 52 ) --Л3 5153(о51 +(52)-оЛ3 (оГ1 +(Г2 + 53)532 = 0, (Л0 53 -Л52 -Л2Г2 - (о2 + (2)Л3 52 (Г151 +
+Г2 52 ) - Л 52 53 (о51 + ( 52 ) - рЛъ (оГ1 + ( Г2 + 53 )53 = 0, Л 53 + ((о2 + (2)Л0 + Л2 ) Г3 = 0.
Это система нелинейных алгебраических уравнений с параметрами о,(,Л0,Л,Л2,Л3. Нас
интересуют решения системы (5.2.3), содержащие одну и более переменных свободными. С помощью программы “НИЬегЮ^епзюп" из пакета программ “Ро!упот1а!!Ьеа!8'' удалось установить, что уравнения (5.2.3) имеют своими решениями множества размерности 1, т.е. решения, содержащие одну свободную переменную. Для выделения этих решений построим базис Грёбнера уравнений (5.2.3) при исключающем мономиальном упорядочении:
((°2 + (2) Л0 + Л2 )Г3 + Л Г3 = 0 53 (оЛ153 + ((°2 + (2) Л0 + Л2 )51 ) / (51, 52 , 53) = 0
((о2 +(2)Л 53 +а((а2 +(2)Л0 +Л2) 51 +(((°2 +(2)Л0 +Л2 ) 52) ./2 (51, 52 , 53) = 0,
53 /з(51, 52 , 53) = 0 53 /4 (51, 52 , 53) = 0 /5 (51, 52 , 53 ) = 0, /б (51, 52 , 53 ) = 0,
52 , 53 ) = а 53/8 (51, 52 , 53 ) = а /9 (Г1, 51, 52 , 53 ) = 0 /10 (Г2 , 51, 52 , 53 ) = 0
здесь / ( = 1,____10)- полиномы степени 2-6 от переменных Г1,Г2,Г3 , 51,52 , 53 .
(5.2.4)
Программа “НИЬегЮ^епзюп" позволяет установить, какие из подсистем системы (5.2.4) имеют решения - множества размерности 1. Ниже приведены некоторые из найденных решений:
(о2 + (2)Л;Г2 52 + Л = о51 + (
(о2+(2)Л) 51 , 3 о2 + (
s Г = + VЛ + (а2 + ^ХЛ - Л^2 ) „ =м^1+ a +Р2)(10 -13Sl2)
, Г2 = + , О о. ПТ , S2 = ±
л]а2 + fb 2 (a +ь2)л1Лз 2 л] (a + fr)1
S3 = 0 при 1 = -(a2+b2)1®-
a(а2 +Ь2)1з si ±>Д +(a +b2)(10 -13si2) „ = Л _ Л 2 , e2
(5.2.6)
(a2 + Ь2Гл/Л
Используя определение инвариантного многообразия и средства СКА Mathematica (или Maple), можно проверить, что (5.2.5) и (5.2.6) определяют семейства ИМСД уравнений (5.2.1).
Исключив из уравнений (5.2.1) переменные Г1,r3,s3 с помощью (5.2.5), получим, что векторное поле на элементах семейства ИМСД (5.2.5) описывается уравнениями
= (a, +bS2)S, + (b(a2 + b2)13S1 Г2 -a((a +b2)13S2 Г2 +^1))((a2 +b2)13S2 Г2 +11)
Г2 2 , 02
a2 + b (a2+b2)2^32 s,2
S = a(1 s2 + (^2 + b )Л3Г2 (s2 + S2 )) s = b(11S2 + (a + b )Л3Г2 (si2 + s2)) (5 2 7)
1 (a2 + b2)1 s1 ’ 2 (a2 + b2)13s1 '
Уравнения (5.2.7) имеют первые интегралы:
2
К = (bs1 aS2 )2 К = (as1 + bS2 )2 + ((a2 +Ь2)Л3Г2S2 +^1)2 + Г2
0 2(a2 +b2)’ 1 (a2 +b2)2 (a2 +b2)1 s? 2'
Для уравнений (5.2.7) рассмотрим задачу выделения ИМСД 2-го уровня (на ИМСД (5.2.5)). Как и в предыдущем случае, составим линейную комбинацию из первых интегралов задачи:
к = m К - m К (m = const) (5.2.8)
и запишем условия (4.3.3) стационарности K по всем переменным:
2 ((а2 + (()(2т - 2а2т К + 2(а2 + (2) ЛтГ5^ + 4(а2 + (2 ЛЛтГ52 - а(Л2 ((а2 + (2 )т0 + +т1)513 52 + 2Л2т = 0, а(Л3 ((а2 + (2 )т0 + 2т к3 - л ((а2 + (2 )а т0 - 2(2т1)512 52 +
+ 2(а2 + (2) Ла^252 + 2(а2 + (2 = 0,
(Л (а + ( )*251 + Л (а + ( )*2 52 + Л52 )т = 0. (5.2.9)
Используя технику базисов Грёбнера, разрешим уравнения (5.2.9) относительно переменных г2,52 и параметров т0,т1. Одно из найденных решений следующее:
Г2 =abЛ/((a2 +(2)Л51), 52 =-о51/(, т =-2(4Л2м/((о2 +Ь2)Л254). (5.2.10)
Проверка показывает, что (5.2.10) определяют семейство ИМСД дифференциальных уравнений (5.2.7). Векторное поле на элементах этого семейства ИМСД описывается уравнением
s0 = 0. (5.2.11)
Следовательно, последнее уравнение (5.2.10) представляет тривиальный первый интеграл дифференциального уравнения (5.2.11). Исследуем на устойчивость решения, принадлежащие элементам семейства ИМСД (5.2.10) и соответствующие параметру s0. Введем отклонения
переменных от их значений в невозмущенном движении: z0 = r2-ap\j((a2 +b2)1si),
z2 = s2 + asjp, z3 = s0 -s0. Вычислим вторую вариацию^2K интеграла K (5.2.8) по этим
отклонениям. Исключив z3 из S2K с помощью уравнения dV0 = (a2 + p2) z3 -a pz2 = 0 , получим квадратичную форму
(a2 +fb)m , + 2(a2-p2Um, /((a2 +р2)УЧ2-a2(a2 + в/Ют 2 (5212)
- /Э2 Z0 + 2 2 Z0 Z2 2 2n5 04 «2 z 2, (5.2.12)
p s0 (a2 +p2)1 (a +p)si 2
знакоопределенную, если выполнен критерий Сильвестра, что приводит к неравенствам
2 +b2)m/b2 < 0, ((a2 +p2)312 s04 -/(/-00a2)1)l ((a2 +p2)42 ^) > 0. (5.2.13)
Неравенства (5.2.13) совместны при выполнении условий
m< 0 л si. ф 0 л ((((V0üa >pAp> 0) v (b< 0 a V0üa + p > 0)) a l ф 0)
/ г- п I- п ,p4p2 - 10a22 „ pjp2 - 10a2 2 0 ллЛ (5.2.14)
v((y[\0a<bvy¡\0a + b< 0)л(—*--------------^ <2 v^2-----------3/г + 2 < 0))).
s0 (a2 +p2) s0 (a2 +p2)
Решения, принадлежащие элементам семейства ИМСД (5.2.10), для которых выполняются условия (5.2.14), устойчивы в смысле Ляпунова. Подобное исследование на устойчивость было проведено и для семейства ИМСД (5.2.6), векторное поле на элементах которого также описывается уравнением вида (5.2.11). Для решений, принадлежащих элементам семейства ИМСД
(5.2.6) и соответствующих параметру s00, получены, кроме того, условия неустойчивости по первому приближению.
6. Заключение
Разработанное программное обеспечение позволяет автоматизировать и, следовательно,
существенно ускорить процессы моделирования и качественного анализа сложных систем,
избежать ошибок на всех этапах исследований. Многолетний опыт работы с программными
комплексами, инструментальными средствам которых являются системы аналитических
вычислений, подтверждает их эффективность и перспективность. Они могут быть использованы как в научно-исследовательской деятельности, так и образовательном процессе.
Работа выполнена при частичной поддержке INTAS-SB RAS, грант №06-1000013-9019.
1. Дэвенпорт Дж. и др. Компьютерная алгебра / Дж. Дэвенпорт, И. Сирэ, Э. Турнье. - Москва: Мир, 1991. -350 с.
2. Грошева М.В. и др. История использования аналитических вычислений в задачах механики / М.В. Грошева, Г.Б. Ефимов, В.А. Самсонов. - Москва: Институт прикладной математики им. М.В. Келдыша, 2005. - 87 с.
3. Бурлакова Л.А. и др. Использование символьных выкладок на ЭВМ в некоторых задачах механики / Л.А. Бурлакова, В.Д. Иртегов, М.В. Почтаренко // Аналитические вычисления на ЭВМ и их применение в теоретической физике. - Дубна: ОИЯИ, 1980. - Д-11 -80-13. - С. 137 - 142.
4. Пакет символьных вычислений для исследования динамики систем тел / А.В. Банщиков, Л.А. Бурлакова, Г.Н. Иванова и др. // ППП. Программное обеспечение математического моделирования. Серия «Алгоритмы и алгоритмические языки». - Москва: Наука, 1992. - С. 63 - 71.
5. Программный комплекс LinModel для анализа динамики механических систем большой размерности / А.В. Банщиков, Л.А. Бурлакова, В.Д. Иртегов и др. // Свидетельство о государственной регистрации программы для ЭВМ № 2008610622. ФГУ ФИПС, 1 февраля 2008 г.
6. Банщиков А.В., Бурлакова Л.А. Информационно-исследовательская система «Устойчивость» // Известия РАН. Теория и системы управления. - 1996. - № 2. - С.13 - 20.
7. Иртегов В.Д., Титоренко Т.Н. Об использовании электромеханических аналогий // Системы поддержки принятия решений для исследований и управления энергетикой. - Новосибирск: Наука, 1997. - С. 100 - 108.
8. Лурье А.И. Аналитическая механика. - Москва: ФИЗМАТГИЗ, 1961. - 824 с.
9. Гантмахер Ф.Р. Лекции по аналитической механике. - Москва: ФИЗМАТГИЗ, 1960. - 292 с.
10. Brayton R.K., Moser J.K. A theory of nonlinear networks -1 // Quarterly of Applied Mathematics. - 1964. - Vol. 22-
1. - P. 1 - 33.
11. Титоренко Т.Н. Уравнения Гамильтона в случае вырожденности динамических систем // Труды IV междунар. конф. женщин-математиков. Спецвып. «Известия ВУЗ. Радиофизика». - 1997. - Т. 4, Вып. 2. -С. 102 - 108.
12. Банщиков А.В. и др. Problems of stability of dynamic systems and computer algebra / Банщиков А.В., Бурлакова Л.А., Иртегов В.Д. // Записки научных семинаров ПОМИ им. Стеклова. - 1999. - Т. 258. - С. 262 -279.
13. Banshchikov A.V., Bourlakova L.A. Computer algebra and problems of motion stability // Mathematics and Computer in Simulation. - 2001. - № 57. - P. 161 - 174.
14. Алгоритмы качественного исследования сложных систем / А.В. Банщиков, Л.А. Бурлакова, В.Д. Иртегов и др. // Кибернетика и системный анализ. - 1992. - № 1. - С. 129 - 139.
15. Румянцев В.В. Сравнение трех методов построения функций Ляпунова // Прикладная математика и механика. - 1995. - Т. 59, Вып.6. - С. 916 - 921.
16. Иртегов В.Д. О теореме Рауса-Ляпунова // Труды IX междунар. Четаевской конф. “Аналитическая механика, устойчивость и управление движением “. - Иркутск, 2007. - Т. 2. - С. 86 - 100.
17. Кокс Д. и др. Идеалы, многообразия и алгоритмы / Д. Кокс, Дж. Литтл, О'Ши. - Москва: Мир, 2000. - 760 с.
18. Банщиков А.В. Исследование устойчивости гирорамы методами компьютерной алгебры // Математическое моделирование. - 1996. - № 4. - С. 67 - 78.
19. Banshchikov A., Burlakova L. On Stability of a Satellite with Gyrodines // Proc. of the 7th Workshop on Computer Algebra in Scientific Computing. - München: Technische Universität, 2004. - P. 61 - 69.
20. Иртегов В.Д., Титоренко Т.Н. О моделировании и исследовании некоторых задач с помощью компьютерной алгебры // Программирование. - 1997. - № 1. - С. 68 - 74.
21. Irtegov V., Titorenko T. On Some Results of Investigation of Kirchhoff Equations in Case of a Rigid Body Motion in Fluid // Springer, LNCS. - 2005. - Vol. 3718. - P. 259 - 271.
22. Соколов В.В. Новый интегрируемый случай для уравнений Кирхгофа // Теоретическая и математическая физика. - Т.129, Вып. 1. - 2001. - С. 31 - 37.
Стаття надійшла до редакції 23.09.2008
