АЛГОРИТМЫ РЕШЕНИЯ УРАВНЕНИЙ ДВИЖЕНИЯ В ИМПУЛЬСАХ ПУАССОНА СИСТЕМ ТВЕРДЫХ ТЕЛ СО СТРУКТУРОЙ ДЕРЕВА Текст научной статьи по специальности «Математика»

2017

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

Математика. Механика. Информатика

Вып. 4(39)

УДК 531.01+004.94

Алгоритмы решения уравнений движения в импульсах Пуассона систем твердых тел со структурой дерева

В. Н. Иванов

Пермский государственный национальный исследовательский университет Россия, 614990, г. Пермь, ул. Букирева, 15 precol@psu.ru; 8(342) 2-396-560

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

Ключевые слова: система абсолютно твердых тел; уравнения движения; динамика; математическое моделирование; обобщенные координаты; импульсы Пуассона; матрично-геометрический метод; разложение Холецкого

DOI: 10.17072/1993-0550-2017-4-25-31

Введение

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

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

© Иванов В. Н., 2017

ляющих положение и распределение скоростей тел системы в пространстве, используются обобщенные координаты, квазискорости и переменные, имеющие размерность обобщенных импульсов (импульсов Пуассона). Статья является продолжением работы [1]. Цель настоящего исследования: сравнение различных алгоритмов разрешения предлагаемых уравнений движения относительно обобщенных скоростей и импульсов по их вычислительной эффективности.

1. Описание механической системы

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

Пусть N - число тел в системе (не считая тела "0", движение которого во времени

относительно инерциальной системы координат (СК) Охуг задано). Тогда количество шарниров всегда будет равно N если учитывать шарнир между системой и телом "0".

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

Положение и ориентацию 7-го тела системы относительно инерциальной СК будем

определять с помощью радиуса-вектора ООг фиксированной в нем точки О и системы ор-

тонормированных векторов в]

(') ^) ^)

е', ко-

торые вместе с точкой О определяют СК ОлУ.г,, неизменно связанную с 7-м телом.

Введем следующие обозначения: р= Ок О - матрица-столбец координат точки О в к -й СК; г = ОО - матрица-столбец координат точки О в инерциальной СК; О' - матрица направляющих косинусов между базисными векторами е\1), е\1), е( 1) и е),

е^), е3(г) (матрица преобразования координат

из ]-й СК в 7-ю).

Предположим, что 7-й шарнир моделируется системой / < 6 идеальных голоном-

ных линейно-независимых связей. Тогда множество относительных положений 7-го тела в к -й СК образует щ -мерное конфигурационное многообразие [2], где щ = 6 — / . Уравнения этого многообразия можно записать в параметрической форме, введя матрицу-столбец д =(д],...,дЩ' криволинейных (обобщенных) координат его текущей точки.

Тогда матрицы рi и О = О' являются функциями обобщенных координат:

Рг =Р (Ч, , *X 0г = 0 (д, , *) .

Задание матрицы-столбца д = (д,..., д^ У

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

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

Г' = гк1 + О0 Р, О'0 = ООк .

Принимая движение 7-го тела за относительное, а предшествующего ему к -го тела

за переносное, в соответствии с правилом сложения скоростей можно записать рекуррентные формулы для вычисления проекций линейной V и угловой ф1 скоростей тел механической системы на оси 7-й СК [3]:

(1)

где vi

С1 -С1Р1

0

(аш

А —

а

( 1

а

1

—с>

а^ — С/ —~7, аУп — С/

( 0 —Р1

( 0

—

дС>

а11 — ТГ^ а10 — ,] — 1,П1.

—*1

Здесь и далее символ "~" используется для обозначения кососимметричной матрицы (это соответствует матричной записи векторного произведения [3]).

Введем блочную 6^6^матрицу с квадратными подматрицами порядка 6 по следующей формуле:

Е6, 7 =', _

— С', 1 = к, г, 1 = 1, N. (2)

06 x6, 1 кг,

Заметим, что для любой кинематической структуры эта матрица содержит в каждой строке только два ненулевых блока: £6 и

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

2

ш

уравнения кинематики системы тел (1) можно записать следующим образом:

Sv = Aq + V*,

(3)

где V = (У!,^!, ...,ум,шм)т, А = diag(A1, ...,АМ).

Рекуррентные формулы (3) можно записать в виде явных выражений:

V = Т(Ац + V*), (4)

где обратная к матрица Т = $ 1 является блочной 6^6^матрицей, подматрицы которой могут быть вычислены по рекуррентным формулам:

Т =

V

Ев,

>

0

6 X 6'

] = I, _

] е р, I, ] = 1, N. (5)

] * Р,

2. Уравнения движения в импульсах Пуассона

Введем следующие обозначения:

М1 =

т^Е —т

П =

Ь =

К

т)

0

6)1

М = diag(M1, ...,Мм), П = diag(П1, Р = (Р1, ...,РМ)Т, Q = АГТтР, где щ - масса /-го тела; Ji - тензор инерции

/-го тела; г[ - матрица-столбец проекций радиуса-вектора центра масс /-го тела на оси связанной с ним СК; /¿0, т° - проекции главного вектора и главного момента активных сил, действующих на /-е тело в /-й СК; Q -матрица-столбец обобщенных сил.

Тогда уравнения движения систем твердых тел со структурой дерева можно записать относительно расширенного множества переменных в виде:

Mv — = 0, —5v + Ац = —V*, Атц = р, р = (Ат — АТП)ц + д.

(6)

(7)

Особенность уравнений (6), (7) состоит в том, что они разрешены относительно производных обобщенных импульсов р (импульсов Пуассона). Первые три уравнения (6) образуют линейную систему с симметричной, блочной трехдиагональной разреженной матрицей коэффициентов относительно скоростей v,q и переменных ^, которые являются множителями Лагранжа [4]. Полный вывод

уравнений движения в форме (6), (7) можно найти в работе [1].

Структура трех кинематических уравнений системы (6) аналогична структуре расширенной системы уравнений динамики системы связанных твердых тел [5, 6]:

+ Ац = —ш*, АТИ = 0,

(8)

где w - матрица-столбец абсолютных линейных и угловых декартовых ускорений всех тел системы в связанных с ними СК. Правые части уравнений Р* и w*, помимо внешних сил, содержат гироскопические, центробежные и явно зависящие от времени ускорения.

Отличие уравнений (6) и (8) состоит в том, что первые уравнения служат для получения значений обобщенных скоростей ц, а не ускорений ц.

Совпадение структуры уравнений (6) и (8) означает, что все разработанные ранее алгоритмы решения системы уравнений (8), такие как метод "прогонки" (отдельных тел) А.Ф. Верещагина [7, 8], методы проекций уравнений в подпространства обобщенных координат или подпространства линейно-независимых компонент реакций связей (схемы явного формирования уравнений движения в форме Лагранжа 2 или 1 рода) [3, 9, 10, 11], применимы и для системы уравнений (6).

Преимущество формы уравнений (6), (7) состоит в том, что при наличии дополнительных геометрических или кинематических связей в механической системе, уравнения (6), (7) могут быть замкнуты только уравнениями кинематических связей, а не продифференцированными дважды уравнениями связей, как в случае использования уравнений (8). При этом методы стабилизации связей (типа метода Баумгарта [10]), применяемые для уравнений (6), (7), обеспечивают экспоненциальный закон компенсации отклонений в случае нарушения связей и не приводят к появлению дополнительных "паразитных" малых высокочастотных колебаний в механической системе. Это означает, что при применении известных численных методов интегрирования дифференциально-алгебраических уравнений к уравнениям движения вида (6), (7) возможен больший шаг интегрирования, следовательно, меньшее время компьютерного моделирования.

В работе [6] получены алгоритмы метода прогонки и Холецкого разрешения системы уравнений (8) относительно старших про-

<

с

ЩЪ

изводных. В настоящей статье аналогичные алгоритмы построены для уравнений (6), (7) в импульсах Пуассона.

3. Решение уравнений движения в импульсах Пуассона методом прогонки

Метод прогонки по существу является модификацией метода Гаусса решения систем линейных алгебраических уравнений с ленточной структурой. В данном методе при прямом ходе, который выполняется начиная с последнего тела системы, из группы уравнений (6) исключаются множители Лагранжа /. При обратном ходе по явным формулам вычисляются последовательно для каждого тела системы обобщенные скорости q, декартовые скорости v и множители Лагранжа /.

Вывод рекуррентных формул метода прогонки подробно изложен в статье [1]. Здесь же приведем только результирующий алгоритм метода.

Алгоритм метода прогонки for i = N: 1

Di = (AjMlAi)-1 Ht = E- MlAiDiA] F^Fi + Zj^CjF/ Ml = Mi+ZjeUiCfHjM;Cj Vl = ZyeUi C[[MjAjDjPj + H/M/v/ + <pj)] end

for i = 1:N

qi = Di [pi — Af (<pi + Mf(CiVki + vi))] Vi = CcJk.+ Aiqi + vi /i = Mlvi + <р\ Pi = (Af — Af fli)/i + AfFi

end

Трудоемкость решения системы уравнений (6), (7) с помощью данного алгоритма растёт по линейному закону в зависимости от числа тел в механической системе. При реализации этого алгоритма требуется обращение только симметричных положительно определенных матриц AfMlAi, порядок которых равен числу степеней свободы в i-м шарнире, причем эти матрицы симметричны и положительно определены, а их порядок всегда мал (не превышает шести). Именно этим и обусловлена эффективность этого метода.

4. Разрешение уравнений движения с помощью блочного ЬТБЬ -разложения

Запишем систему уравнений (6) в виде единого матричного уравнения. Для системы тел, имеющей структуру простой цепочки, оно будет иметь следующий вид:

0 \ / ИЦ /ф Я Я - 0 \1 Щ>\_1 ф

, (9)

0 0

D

N'

.и*,/

где

о A]

0 0 0

Di = ( A,- 0 —E ), Hi = (0 0 Ci

E M Pi

0 0 0

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

Идея исключения Гаусса - это преобразование исходной системы уравнений в эквивалентную треугольную систему с последующим ее решением. В матричном виде это есть ¿^-разложение [12]. В случае симметричной матрицы системы уравнений (9) лучше использовать

-разложение. Заметим, что при разложении ленточных систем треугольные множители Ь также являются ленточными.

Первый шаг алгоритма блочного ЬТОЬ-разложения, примененного к блочной трех-диагональной матрице системы (9), можно записать в виде

- 0 0

м = м0 =

0 — DV-1

0 — HN Dn со следующими матрицами:

Hi l = ^MA

L =

м1 =

E

0 0

D1

0 0

juw

00

E 0

pn Ln/ 00

D

N-1

ln = i lm-ian —imi-1 0 — l-mTn

-MN,

0

0

0

0

0

0

0

'0 0

V

N

PN = \0 0 L'J-iC,

DN-1 =

00 0

An-i 0

M-0

AN-1 0

—E

-N

0 E

N-1J

MN—

E 0

0

/=10

E 0

где Lm

LM-i и L—w - факторы Холецкого

\0 0 Е

, ¿м-

симметричных положительно определенных матриц Мм, М^1 и — соответ-

ственно, матрица находится из матричного

уравнения — А^МмСм, а матрица М*^_1

вычисляется по формуле Мдт_1

— М„_ 1 +

Нетрудно видеть, что матрица имеет ту же структуру, что и матрица и следовательно, процесс разложения матрицы М на множители можно продолжить:

М — М0 — ![М1!1 — цц. М2Ц2Ц1 — ••• — цц. - .„ 1211 — 1ГИЦ

где

L =

=

L1 0 P2 L2

00 7 0

0 I

00

0 0

L2/ 0 0

'Щ 0

li = ( lmt-lai —lm-1

-l—ml

00

lmu

V

Р1 — (0 0 1~М-гС{ ),

0 0 0

М? — М{ + - — ¿м 1м,

и? — А[М?А( — ^ — ^[М?^

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

— Ж,

— Ж,

(10) (11)

где X = (Xi.....Xn)t, Xi = (х0 ylt z),

F = (01,., 0N)T,W = (W1,., WN)T.

Решая систему (10), начиная с последнего уравнения, получим:

= — L

Mi

X

* i+1

(Ci+1( p*

j —T * = —LMíФi,

+ M*

i+1vi+1 + ViT+1Xi+1))

L* i т —T *

.i M- l^i+LMГlVi,

Xi = Ltf(pi — AT(<p; + M*v;)).

Решая систему (11), и исключая переменные z i, yi, получим:

4i = L—]1i(Xi — ViVi—1), vi = CivkT + Aiqi + V* Mi = M*vi + pi*.

Обобщая полученные формулы на случай системы твердых тел со структурой дерева, запишем следующий алгоритм решения расширенной системы уравнений (6), (7).

Алгоритм метода Холецкого

for i = N: 1

F[ = Fi+Zje—iCjF]* M* = Mi + Zjes^MjCj — V/V-) P* = ^jBSi(cJ (Mj vj + p**) + Vjxj) Найти фактор Холецкого L—. матрицы У*: L—iL—i = У* = AM A, 1 Решить L—.Xi = pi— AT(p* + Mfv*) Решить L—.Vi = ATiM**Ci end

fori = 1:N

Решить L—.qi = x, — ViVki vi = Civki + Aiqi + vi* Mi = M*vi + pi* Pi = (AJ—Ajni)Mi+AjF[ end

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

5. Сравнение эффективности методов

Приступая к сравнению рассмотренных выше методов разрешения системы диффе-

0

0

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

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

Метод прогонки

с 80

I

1_Г 70

го 3 60

50

и

40

о

го 30

I

о; 20

S е 10

Q. CQ 0

5 10 15 20 25 30 35 40 45 50 Число тел

Затраты на разрешение уравнений движения III Затраты на вычисление кинематических соотношений ■ Затраты численного метода интегрирования

Рис. 1

Анализ этих зависимостей подтверждает теоретические расчеты. Видно, что затраты времени растут линейно от числа тел в системе. Причем затраты времени на разрешение уравнений движения относительно ускорений занимают львиную долю всех временных затрат.

Для сравнительной оценки эффективности методов построим график отношения T / T2 (рис. 3), где T и T2 - временные затраты на разрешение уравнений движения методом прогонки и методом L DL -разложения соответственно.

Нетрудно видеть, что метод, основанный на симметричном LT DL -разложении, эффективнее метода прогонки, не учитывающего особенности уравнений движения системы тел. Из рис. 3 видно, что алгоритм

Сравнение описанных выше алгоритмов разрешения уравнений движения проводилось с помощью компьютерного моделирования динамики механических систем, представляющих собой цепочки твердых тел, соединенных двухстепенными кардановыми шарнирами. При сравнении варьировалось число тел в цепочке. Дифференциальные уравнения решались методом Штёрмера [13] с порядком аппроксимации равным шести. Эффективность методов измерялась временем выполнения одного шага этого численного метода.

На рис. 1 и 2 представлены зависимости временных затрат на выполнение одного шага интегрирования от числа тел в цепочке для рассматриваемых методов разрешения уравнений движения.

Метод ¿тО£-разложения

и 80

5 10 15 20 25 30 35 40 45 50 Число тел

К Затраты на разрешение уравнений движения

III Затраты на вычисление кинематических соотношений

■ Затраты численного метода интегрирования

Рис. 2

L DL -разложения в 1,3 - 1,4 раза сокращает временные затраты на разрешение уравнений движения по отношению к методу прогонки.

1.5 1.4 1.3 I-T 1.2 1.1 1

5 10 15 20 25 30 35 40 45 50 Число тел • Экспериментальные данные -Линия тренда

Рис. 3

Заключение

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

Предложены два рекуррентных алгоритма разрешения этих дифференциальных уравнений относительно старших производных при их численном интегрировании: метод прогонки и метод, основанный на Ь ПЬ -разложении Холецкого. Проведено сравнение предложенных алгоритмов по их вычислительной эффективности. На примерах моделирования механических систем с различным числом степеней свободы показано, что алгоритм ЬПЬ -разложения ускоряет расчеты по сравнению с методом прогонки на 30-40 %.

Список литературы

1. Иванов В.Н., Полосков И.Е., Шимановский В.А. Математические модели систем связанных твердых тел в импульсах Пуассона // Фундаментальные исследования. 2016. № 10-3. С.493-499.

2. Величенко В.В. Матрично-геометрические методы в механике с приложениями к задачам робототехники. М.: Наука, 1988. 280 с.

3. Лурье А.И. Аналитическая механика. М.: Гос. изд-во физ-мат. лит., 1961. 824 с.

4. Иванов В.Н., Суслонов В.М. Уравнения движения механических систем со структу-

рой дерева // Проблемы механики управляемого движения. Нелинейные динамические системы. 1984. С. 154-159.

5. Иванов В.Н., Домбровский И.В., Набоков Ф.В. и др. Классификация моделей систем твердых тел, используемых в численных расчетах динамического поведения машиностроительных конструкций // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2012. № 2. С. 139-155.

6. Шимановский В.А. Метод компьютерного моделирования динамики систем связанных твердых тел // Фундаментальные исследования. 2017. № 8. С. 104-109.

7. Верещагин А.Ф. Метод моделирования на ЦВМ динамики сложных механизмов роботов-манипуляторов // Известия АН СССР. Техническая кибернетика. 1974. № 6. С. 89-94.

8. Верещагин А.Ф. Принцип наименьшего принуждения Гаусса для моделирования на ЭВМ динамики роботов-манипуляторов // Докл. АН СССР. 1975. Т. 220, № 1. С. 51-53.

9. Лилов Л.К. Моделирование систем связанных тел. М.: Наука, 1993. 272 с.

10. Wittenburg J. Dynamics of multibody systems. Berlin: Springer-Verlag, 2008. 223 p.

11. Shabana A.A. Computational dynamics. New York: Wiley, 2009. 542 p.

12. Golub G.H., Van Loan C.F. Matrix Computations. The Johns Hopkins University Press, 2012. 790 p.

13. Hairer E., Norsett S.P., Wanner G. Solving ordinary differential equations I. Nonstiff Problems. Springer-Verlag, 2011. 528 p.

Algorithms for solving the equations of motion with Poisson impulses of multibody systems with tree structure

V. N. Ivanov

Perm State University; 15, Bukireva st., Perm, 614990, Russia precol@psu.ru; 8(342) 239-65-60

The article presents the new matrix form the equations of motion with generalized coordinates and Poisson impulses of multibody systems with tree structure. Two algorithms solve these equations for senior derivatives focused on the use of computers are offered. The first algorithm is a method of Chaser for solving the systems of linear algebraic equations. The second algorithm uses the Cholesky factorization schema. Recurrent formulas are obtained for all kinematic and dynamic variables that are included in the equation. The time complexity of solution of equations using data algorithms grows on linear depending with the number of bodies in a mechanical system that testifies to their effectiveness. Proposed algorithms are compared. For examples of integrating the equations of motion of systems of bodies with a large number of degrees of freedom shows the advantage of algorithm based Cholesky decomposition schema.

Keywords: multibody system; equations of motion; dynamic; numerical methods; generalized coordinates; Poisson impulses; matrix-geometric decomposition method; Cholesky decomposition.