CC BY
5_вестник пермского университета_
2022 • Математика. Механика. Информатика • Вып. 4(59)
Механика
Научная статья
УДК 531.01+004.94
DOI: 10.17072/1993-0550-2022-4-18-28
Уравнения движения в гамильтоновых переменных систем твердых тел с замкнутыми кинематическими цепями
Владимир Николаевич Иванов
Пермский государственный национальный исследовательский университет, Пермь, Россия precol@psu.ru, https://orcid.org/0000-0002-2634-4296
Аннотация. Рассматриваются методы построения математических моделей механических систем, в структуре которых содержатся кинематические циклы. Все связи голономны. На первом этапе механическая система сводится к системе со структурой дерева с помощью раздвоения определенных вершин (тел) первичного графа. Уравнения движения строятся в матричной форме относительно гамильтоновых переменных: обобщенных координат и импульсов. На втором этапе дополнительные связи, приводящие к замыканию кинематических цепей, учитываются в уравнениях движения с помощью множителей Ла-гранжа на уровне кинематических соотношений. Уравнения движения формируются по рекуррентным формулам с использованием минимального набора первичной информации о структуре, геометрических, кинематических, масс-инерционных и силовых характеристиках механической системы. Построены итерационный и конечный алгоритмы решения полученных уравнений относительно переменных интегрирования. На примере одной механической системы, в структуре которой содержатся два замкнутых цикла, продемонстрирован весь процесс подготовки первичной информации при формировании математической модели в предлагаемой форме. Результаты моделирования показывают, что алгоритмы, описываемые в статье, приводят к более устойчивому к ошибкам вычислений численному решению, чем методы, основанные на классическом подходе, в котором уравнения движения замыкаются дважды продифференцированными уравнениями дополнительных связей. Ключевые слова: система абсолютно твердых тел; уравнения движения; динамика; математическое моделирование; обобщенные координаты; импульсы Пуассона; матрично-геометрический метод Для цитирования: Иванов В. Н. Уравнения движения в гамильтоновых переменных систем твердых тел с замкнутыми кинематическими цепями // Вестник Пермского университета. Математика. Механика. Информатика. 2022. Вып. 4(59). С. 18-28. DOI: 10.17072/1993-0550-2022-4-18-28.
Статья поступила в редакцию 24.05.2022; одобрена после рецензирования 18.10.2022; принята к публикации 30.10.2022.
Mechanics
Research article
Motion Equations of Rigid Bodies Systems with Closed Kinematic Chains in Hamiltonian Variables
Vladimir N. Ivanov
Perm State University, Perm, Russia
precol@psu.ru, https://orcid.org/0000-0002-2634-4296
Abstract. Constructing methods of mathematical models for mechanical systems with closed kinematic loops are considered. All connections are holonomic. At the first stage, the mechanical system is reduced to a system with a tree structure by certain vertices (bodies) splitting in the original graph. The motion equations are constructed in matrix form with respect to Hamiltonian variables: generalized coordinates and impulses.
Эта работа О 2022 Иванов В. Н. лицензируется под СС ВУ 4.0. Чтобы просмотреть копию этой лицензии, посетите http://creativecommons.Org/licenses/by/4.0/
At the second stage, additional constraints are taken into account in the motion equations using Lagrange multipliers at the level of kinematic relations. These connections lead to closure of kinematic chains. The motion equations are formed according to recurrent formulas using a minimum set of primary information about the structure, geometric, kinematic, mass-inertial and power characteristics of a mechanical system. Iterative and finite algorithms for solving the obtained equations with respect to the integration variables are constructed. The whole process of primary information preparation is demonstrated in the formation of a mathematical model in the proposed form on the example of one mechanical system with two closed cycles. The simulation results show that the algorithms described in the article provide a numerical solution that is more resistant to computational errors than methods based on the classical approach, in which the motions equation are closed by doubly differentiated equations of additional constraints.
Keywords: multibody system; equations of motion; dynamic; mathematical modeling; generalized coordinates;
Poisson impulses; matrix-geometric method
For citation: Ivanov V. N. Motion Equations of Rigid Bodies Systems with Closed Kinematic Chains in Hamil-
tonian Variables // Bulletin of Perm University. Mathematics. Mechanics. Computer Science. 2022;4(59):18-28.
(In Russ.). DOI: 10.17072/1993-0550-2022-4-18-28.
The article was submitted 24.05.2022; approved after reviewing 18.10.2022; accepted for publication 30.10.2022.
Введение
Рассматриваются методы построения математических моделей систем абсолютно твердых тел (САТТ), граф структуры которых содержит замкнутые циклы. Выводятся уравнения движения в матричной форме относительно избыточного числа переменных состояния механической системы: обобщенных координат, квазискоростей, импульсов и множителей Ла-гранжа. Аналогичный подход к моделированию САТТ со структурой дерева применен в работах [1, 2]. Кроме того, в [1] описана взаимосвязь рассматриваемого подхода с классическими методами формирования математических моделей САТТ. Настоящее исследование распространяет полученные ранее результаты на САТТ с замкнутыми кинематическими цепями.
1. Описание структуры САТТ
В традиционном подходе структуру любой САТТ представляют в виде графа. Вершинами графа являются тела системы, а ребрами -шарнирные соединения между телами. Выделяют САТТ со структурой дерева и с замкнутыми циклами или цепями. Различия состоят в количестве путей, связывающих любые две вершины графа системы, т.е. является ли путь единственным или нет (рис. 1).
a)
б)
Пусть САТТ имеет древовидную структуру. Пронумеруем ребра и вершины так, чтобы для любых из них номер дочерних ребер и вершин был больше. Такую нумерацию назовем правильной. Тогда структуру САТТ можно задать с помощью одного массива к = км |, в котором к. - индекс вершины или ребра, предшествующих 7-му, / = 1, N, N - число тел в САТТ.
Пусть в структуре первичного графа САТТ есть циклы. На первом этапе введем вспомогательную систему, вторичный граф которой имеет структуру дерева. Вспомогательную систему построим из исходной САТТ, выполнив раздвоение одной из вершин в каждом замкнутом цикле. Соединим пунктирными линиями все пары образов и прообразов раздвоенных вершин первичного графа (рис. 2).
Рис. 1. Примеры САТТ: а) дерево, k = {0,1,1,2,3,3}; б) с замкнутыми циклами
а) б)
Рис. 2. Пример САТТ с замкнутыми циклами: а) первичный граф; б) вторичный граф, к = {0,1,1,3,3,5,2,3}
Пусть во вторичном графе со структурой дерева выполнена правильная нумерация вершин. Для упрощения записей выражений в последующих формулах введем два вспомогательных массива {р }, {и } .
Элементы массивов определим следующим образом: р - множество индексов ребер, составляющих путь между нулевой и 7-й вершинами; и - множество индексов ребер, для которых 7-я вершина является предшествующей.
Условия существования во вторичном графе раздвоенных вершин зададим с помо-
щью матрицы s =
Ь)
ненулевые элементы
которой определяют дополнительные связи между вершинами. Пусть в матрице е число строк совпадает с числом замкнутых циклов, а число столбцов равно числу вершин вторичного графа. На месте элемента е.. поставим
"+1", еслиу - номер вершины первичного графа, принадлежащей 7-му циклу, и "-1", если у -номер вершины вторичного графа, являющейся образом соответствующей раздвоенной вершины 7-го замкнутого цикла первичного графа. Каждая строка в матрице е содержит только по одной паре ненулевых элементов 1 и -1. Например, для САТТ, представленной на рис. 2, матрица е имеет вид
s =
( 0 0 0 1 0 0 0 0
0 -10 ^ 1 0 -1
2. Уравнения движения механических систем с замкнутыми цепями
Для вывода уравнений движения САТТ с замкнутыми цепями возьмем за основу методику формирования матричных уравнений движения САТТ, изложенную в статьях [1, 2].
Пусть в САТТ выполнено размыкание замкнутых циклов и построен вторичный граф системы древовидной структуры.
Введем обозначения: и. - число обобщенных координат /-го шарнира; q = col(01,...,qN) , p = col(р^...,pN) - векторы-столбцы обобщенных координат и импульсов (импульсов Пуассона);
q, = (q)^ql') , pt = (p),...,p") - вект°ры-
столбцы соответствующих координат /-го шарнира; v = col(v1,...,vN) , ¡u = col(p1,...,^N)
- векторы-столбцы шестимерных абсолютных декартовых линейных и угловых скоростей и импульсов (множителей Лагранжа) всех тел САТТ; p = p( q¡, t) - радиус-векторы начала /-й системы координат (СК), в k. -й СК;
Gj = Gt (qt, t) - ортогональные матрицы преобразования координат из k -й в /-ю СК (матрицы поворота); v, = col(%, a¡) , где % , ai -
векторы-столбцы абсолютных линейной и угловой скоростей тел САТТ в проекциях на оси /-й СК;
?1. Q = diag(Qb...,QN),
где coi , vi - кососимметрические матрицы, используемые в матричной записи векторного произведения.
Кинематические уравнения САТТ, связывающие абсолютные и относительные скорости всех тел системы, можно выписать либо в неявном виде
-Sv + Aq = -vr
(1)
либо в явном виде
v = т(Aq+vr) = Bq + (p. (2)
В уравнениях (1) используется матрица кинематической структуры £ = ( ^ ), которая состоит из следующих элементов:
S] =
E
] =', _
j = kt, ', j = 1, N.
06X6, j Ф ' V ki ,
'6 -с,
Матрицы Ci =
G, -Gtpt
0
G,
задают перенос-
ные движения, а матрицы A = diag (A1,..., AN ) и vr = col(v[,...,vrN) определяют относительные движения тел в системе (] = 1, n) :
vr, =
V®* У
A =
a% = Gx
dq]
% = G P ' ' dt
cG
cG
дЧ1 1 ' а
В уравнениях (2) используются матрицы: В = ТА - касательного локального базиса пространства возможных перемещений всех
тел системы и р = Туг - вектор-столбец явно зависящих от времени компонент вектора скоростей V.
Элементы обратной к матрицы Т = могут быть вычислены по рекуррентным формулам [1]:
T =<
ч
E6, 0,
6 х 6,
j = i, _
j е Pk , i, j = 1, N
j * P,
Будем считать, что в раздвоенных телах системы координат расположены одинаково. Тогда условия замкнутости кинематических цепей можно записать в виде дополнительных соотношений, определяющих равенство шестимерных координат и скоростей V раздвоенных вершин в абсолютной СК:
евТ^ = евТа(Вд + <р) = 0,
S7T =S\
(3)
ЛtGTvdz = 0, (4)
. а у '
где Оо = diag [о',..., О^ ) - матрица преобразования шестимерных координат из абсолютной системы координат в СК, связанные с телами;
О" =
fG0 ^
0
О0 у
О'0 = ОО
Масс-инерционные характеристики и действующие силы произвольным образом распределим между раздвоенными вершинами вторичного графа.
В работе [2] уравнения движения механических систем с замкнутыми циклами выведены из принципа Гамильтона-Остро-градского. При этом связи (3) учитывались с помощью множителей Лагранжа т] . Уравнения выписаны в гамильтоновых переменных и имеют следующий вид:
Mv-STju = 0, -Sv + Aq = -vr, ATju = p*+BTG sTr¡, sGlv = 0,
p* =(ÁT-ATá)n + Q,
(5)
(6)
mE
J<
где M = diag (Mx ,...,MN ), Mt =
F = col(Fl,...,FN),Fl=col(fl°,m°\, v = col(V..,vw) ; Q = col(Ql5...5Qn) = AtTtF;
mi, J - масса и тензор инерции; rc - радиус-вектор центра масс /-го тела, заданный в системе координат, связанной с ним; f° , mo - главные вектор и момент сил, действующих на /-е тело и заданных в /-й СК; Q - вектор-столбец обобщенных сил.
Уравнения (5), (6) выписаны относительно расширенного множества гамильтоновых переменных {q, p} и {v,/} при известных множителях Лагранжа ц.
В [1, 2] отмечается, что такая форма записи уравнений движения облегчает процесс алгоритмизации их разрешения относительно переменных интегрирования в процессе компьютерного моделирования динамики САТТ.
Кроме того, дополнительные связи, обеспечивающие замыкание контуров вторичного графа, входят в уравнения (5), (6) в виде кинематических уравнений (3), в отличие от классического метода, в котором замыкание уравнений движения происходит на уровне ускорений, т.е. продифференцированных один раз уравнений (3). Приведенный в работе пример исследования динамического поведения одной механической системы с замкнутыми контурами иллюстрирует тот факт, что при таком подходе ошибки вычислений не приводят к сходу численного решения с геометрических связей (4).
3. Методы разрешения кинематических уравнений (5)
Выпишем алгоритм разрешения уравнений (5) относительно обобщенных скоростей с/ . Будем считать, что вектор // задан. Тогда (5) является линейной системой уравнений относительно скоростей v, q и импульсов / с матрицей системы
rM -ST 0 > -S 0 A 0 AT 0
V /
которая является симметричной, квазитрех-диагональной, разреженной матрицей.
Для приведения системы (5) к двухдиаго-нальному виду можно использовать несимметричное (метод Верещагина) или симметричное (метод Холецкого) гауссово исключение [3, 7]. В данной работе для разрешения системы (5) относительно вектора q будем использовать несимметричное гауссово исключение.
Алгоритм решения системы (5)
for i = N: 1
M = Mi+Z mCTHJMjCi
D, = (ATM*At)-1, H = E-M^DA
Р, = Р, + AJ/*i
CT
j^U, J
MJAjDjPj + Hj (Mjvj + *j )
end for
i = 1: N
q =D
ll I
= M* V,
end .
Построим теперь итерационный алгоритм определения вектора ц из условия удовлетворения решения системы уравнений (5) дополнительным кинематическим связям (3). Во-первых, заметим, что, если из первого уравнения системы (5) выразить вектор v, из второго уравнения - вектор q, из третьего -вектор л и подставить их последовательно в четверное уравнение системы (5), то получим систему линейных уравнений относительно только вектора ц:
sGTav = Фц + b = 0, (7)
где Ф = sGT B (BTMB )-1 BT Ga sT - симметричная, положительно полуопределенная матрица, b = sGa (B (BTMB )-1 ( p* - BTMp) + p) .
В работе [2] для определения вектора ц вместо уравнения (7) использовалось уравнение с несимметричной матрицей системы.
В общем случае система (7) может иметь бесконечно много решений. Однако в силу исходной механической постановки задачи, в которой реализуемая траектория движения доставляет минимум действия по Гамильтону, следует, что решение системы уравнений (7) должно быть найдено в средне-квадратическом смысле, т.е. как нормальное псевдорешение, имеющее минимальную длину [4]. Обычно такое решение ищут методом сингулярного разложения матрицы Ф и построения псевдообратной к ней матрицы [4, 5], что приводит к дополнительным вычислительным затратам.
Кроме того, как видно из описания, вычисление матриц Ф и b требует множества дополнительных матричных операций.
В связи с этим такой способ вычисления вектора ц с помощью решения уравнения (7)
не является оптимальным. Построим алгоритм итерационного построения и решения системы уравнений (7).
Раздвоение вершин первичного графа САТТ будем проводить таким образом, чтобы во вторичном графе со структурой дерева все раздвоенные вершины имели бы по одному различному прообразу и масс-инерционные характеристики раздвоенных тел совпадали.
Тогда, сделав замену переменных
/л = /-TTGasTц, систему (5) можно привести к виду:
Mv-STju* =GasTri, -Sv + Aq = -vr, Ат ¿и* = p*, sGT v = 0.
(8)
(9)
Выразим из первого уравнения системы
(9) скорость V и подставим в четвертое уравнение. В результате получим, что вектор т]
*
связан с вектором / соотношением:
Т = -(еМ- 1еТ )-1 еGTaM- 1£Т/,
где М-1 = ОТМ~1Оа .
В последней формуле сделаем обратную замену / на вектор л :
Т = ~(еМ~ 1еТ )_1 еОТаМ-1 БТ (/ -ТТGaеТт) =
= Т - (еМ- 1еТ )-1 еОТаМ-1 Ба/ =
= Т-(еМ- 1еТ )"1еО>.
Выведенную формулу можно записать в итерационной форме:
Т( к) =Т( к-1) ~{еМ~а е)-1е^^ , (10) где к - номер итерации, т]^ - "новое" уточненное значение вектора т], и -
"старые" значения векторных переменных.
Легко показать, что формула (10) является формулой простых итераций [4, 5] приближенного решения системы линейных уравнений (7):
Т(к) =(Е"(еМа1 ет )1 Ф)Л(к- 1 ) " ( еМ^1 е Т )"1 ь .
Для этого достаточно выразить в формуле
(10) сомножитель еОТ^^^ из равенства (7).
Из теории известно, что последовательность {]} сходится к точному решению системы уравнений (7), если любая из норм матрицы (Е -(еМ- 1еа ) 1 ф) меньше единицы.
В предлагаемом далее алгоритме нахождения вектора т не требуется сходимость итераций, определяемых формулой (10). Однако заметим, что для многих технических систем последовательность {])}
сходится к точному значению т . Косвенным обоснованием этого факта служит то, что матрица (еМ-1еТ) 1 содержит масс-
инерционные характеристики только раздвоенных тел механической системы, а матрица Ф включает в себя обратные матрицы к матрицам, содержащим суммы масс-инерцион-ных характеристик всех вершин, входящих в замкнутые циклы первичного графа. Соответственно, норма произведения этих матриц оказывается меньше единицы. Обозначим:
V/(]) = еОТ V = Ф т] + Ь .
Тогда
/ (](к)) " V/(0) = О (^к) " V(o)) = ФТ(к).
Пусть = 0, и на каждой к-й итерации сначала вычисляются векторы ^^, с помощью представленного выше несимметричного гауссова исключения системы уравнений (5), а затем векторы вычисляются по формуле
(10). При этом матрицы М*, Д и Н вычисляются только один раз на первой итерации.
В процессе итераций сформируем мат-
V,
(V/(^(k))-Vf(0)} и X = k)} .
рицы АР = {V/
Процесс заканчиваем, когда число столбцов в матрице АР становится равным или большим числу строк в матрице Ф (длине вектора т]). В результате получаем систему линейных уравнений: АР = фх. Из этой системы следует, что псевдообратную матрицу Ф+ можно вычислить по формуле
Ф+= X АР + . (11)
При этом псевдообратную матрицу АР+ находим с помощью сингулярного разложения [4].
Зная матрицу Ф+ , вектор т] находим как нормальное псевдорешение системы уравнений (7):
T = -Ф+Vf ( 0) . (12)
Описанный метод разрешения системы
ф,)
уравнений (5) относительно неизвестных переменных f] , q и v окончательно можно оформить в виде следующего алгоритма.
Алгоритм решения уравнений (5), (6)
Функция NPF вычисления вспомогательных векторов:
[Л , Р, Ф }:=NPF for i = N: 1
Л* = sj Л + Z jeUCTrf-
p, = p*+ а!л*
Ф =Z [M ADp+ H (m; v; +
end
end NPF.
Функция QVвычисления скоростей: {q,v}=QV
for i = 1: N
\=Сг\ +Mi+Vi
л = m>, +Ф
end
end QV.
Подготовительный блок: for i = N: 1
M = M +1 ]ucthmtj d, =(atm*a,)-1, H = E-mada
end
Л = 0
i >)' P(o)' 0(o) { ?(o)' v(o) >= QV
V/(0) = sGTaV(0).
Основной цикл: for k = 1: length (л)
Л(к)=Л(k-1)-(sMa"1ST )-1 V/(k-1) {ikyP{k)^{k))=NPFAq(^v(k,}=QV
V/ (k ) = sGTav,
(k)
end
AP = {V/ (k)-V/ (0)}, X = {^k)}
Сингулярное разложение:
aFT = UiVT , AF+ = UIT1VT, Ф+ = XAF+.
Решение СЛАУ (5):
Т = -Ф+Vf ( 0 )
in , pJ}=NPF, {q,v}=QV. Вычисление правой части уравнения (6): for i = N: 1
M = M* vi + Ф
end.
Замечание. Во многих случаях можно в итерационном алгоритме вместо матрицы
(sM- 1sT) 1 использовать более просто вычисляемую матрицу
1 (sMsT ) .
4. Пример
Применим рассмотренную методику для моделирования динамики одной механической системы - шагающего механизма (ШМ), состоящего из четырех тел: корпуса механизма, содержащего элементы привода (тело 1), диска-кривошипа (тело 2), приводящего в действие ноги шагающего механизма и двух ног (тела 3 и 4). Каждая нога, корпус и диск образуют в совокупности кривошипно-шатунный механизм, с помощью которого осуществляется подъем и перенос ног.
Расчетная схема и первичный граф системы представлены на рис. 3, 4.
Взаимодействие каждой ноги с опорной поверхностью осуществляется по дуге окружности с радиусом, совпадающим с длиной ноги. Задача рассматривается в плоской постановке. Ось Ох перпендикулярна плоскости движения. Система построена таким образом, что, в случае фиксации положения ног в механизме, тот совершает устойчивые угловые колебания на опорной ноге.
Взаимодействие ног с опорной поверхностью моделируется с помощью упруго-демпфирующих сил и сил трения таким образом, чтобы обеспечивалось качение рабочей ноги по опорной поверхности без проскальзывания.
— Кривошипно-шагунные механизмы
Нога 1
Рис. 3. Шагающий механизм (ШМ)
Рис. 4. Первичный граф ШМ
ШМ включает две замкнутые цепи. Размыкание проведем раздвоением ног на два тела: на тела 3-5 и 4-6. Новая механическая система со структурой дерева и ее структурная схема представлены на рис. 5, 6. Масс-инерционные характеристики ног разделим пополам между раздвоенными вершинами. Взаимодействие с опорной поверхностью осуществляют тела 5 и 6.
Рис. 5. Схема размыкания ШМ
Рис. 6. Вторичный граф ШМ
к = {0,1,1,1,2,2}, е = -
Для вторичного графа ШМ структурные матрицы имеют вид:
'0010 -1 0 1 0 0 0 1 0 -1].
Начала СК, связанных с корпусом и диском, совпадают и находятся в центре крепления диска к корпусу. Начала СК раздвоенных ног находятся в их вершинах.
Пусть е1 = 1.5 м - длина ноги,
е2 = 0.05 м - радиус точки крепления ноги на диске, е3 = 0.3 м - координата на ноге точки соединения ноги с диском, щ = 20 , щ = 3 , щ = щ = щ = щ = 0.6, г£ =( 0,0, -0.1),
^хх = и Луу = 07, ^ = 0.5,
^ = 0.004, Згуу = 0.0025, J2z = 0.0025,
Гс = Г5С = (0,0,-0.75), г4с = г6с =(0,0,-0.75),
хх = ^Ахх = хх = хх = 0.4675 , уу = '4уу = ^уу = Луу = 04675 ,
^ = = и = ^ = 000015 - массы (кг), координаты центров масс (м) и моменты инерции (кгм2).
Опишем обобщенные координаты ШМ. Корпус (тело 1) имеет три степени свободы: линейные перемещения - у (по горизонту), г (по вертикали) и угловой поворот вокруг оси Ох - р. Диск (тело 2) вращается вокруг оси Ох, угол поворота - а . Тела 3 и 4 (прообразы раздвоенных ног) совершают возвратно-поступательные движения вдоль вертикальной оси Ог в СК корпуса (тела 1), соответствующие обобщенные координаты обозначим через г\ и гг. Кроме того, тела 3 и 4 совершают угловые перемещения относительно тела 1, соответствующие обобщенные координаты обозначим через у\ и у 2. Угловые перемещения тел 5 и 6 (образов раздвоенных ног) относительно диска (тела 2) определим углами /3\ и Д. Векторы обобщенных координат:
91 =(^ Z <Р), 92 = (a), Ь =(Г\),
94 =( Z2, ^2 ) , 95 =(Д ) , 96 = (P2 ) . Система имеет десять степеней свободы. Зададим радиус-векторы рг = Ok^ Oi:
P1 =( 0, y, z + ¿ - ¿3) , P2 =( 0,0,0), P3 =(0,0, Z + ¿3 ) , P 4 =(0,0, Z2 + ¿3 ) , P =( 0, - e3 sin ft, - e2 + e3 cos P), р = (0, - e3 sin P2, e2 + e3 cos P2).
G =
Матрицы направляющих косинусов:
(10 0 ^ (10 0 ^ 0 cosp sinp 0 - sinp cosp
G =
0 cos a sin a 0 - sin a cosa
(
G3,4 =
(
G5,6 =
1 0 0
0 cos sin Yx 2
0 -sin ;к12 cos
1 0 0
0 cos Pl2 sin ft 2
0 - sin P 2 cos P 2
Л
Л
Зададим кинематические матрицы относительного движения A :
(
A =
A3,4 =
0 0 0 1 í 01
cosp sinp 0 0
- sinp cosp 0 0 0 1 = 0 1
0 0 0 0
0 0 0 V l0 V
í 0 0 1 í 0 1
sin^1,2 0 ¿3
cos Y\ 2 0 A - = 0
0 , 1 ? 5,6 1
0 0 0
l 0 0 V l 0 V
А =
0
0
0 1
-<psinq> <pcoscp 0 -<pcos(p -tpsinq) 0
0 0 0
0 0 0
Л.4 ~
( 0 0^1
yh2 cos/L2 0
-yh2únyl2 0
о о
0 0
0 0
6x1 •
Считаем, что ШМ находится в однородном поле силы тяжести. Ускорение свободного падения g = 9.8 м / с2.
Поступательное движение ШМ происходит в результате вращения тела 2 с угловой
скоростью а0 = рад / с и периодом Т = 2.6 с под действием вращающего момента та = -6000(а - а0 ■ ^ - 0.4(а -сг0) и'л' •
Дальнейшее формирование матричных блоков уравнений (5), (6) может быть выполнено по рекуррентным алгоритмам, описанным в предыдущих разделах работы.
Численное моделирование выполнялось в САВ Mathemat7ca. Ставилась задача определения скорости вращения диска, обеспечивающей устойчивую ходьбу ШМ. Результаты исследования представлены на рис. 7 а, б, в.
На рис. 8, для примера, приведены графики модулей отклонений абсолютных координат точек контакта между телом 2 (диском) и телом 3 (ноги 1).
Рис. 8. Отклонения в дополнительных связях
Графики на рис. 8 показывают, что представленный алгоритм решения задачи, основанный на использовании кинематических связей (3), обеспечивает хороший уровень точности в удовлетворении дополнительных геометрических связей (4) на длительном промежутке интегрирования даже без дополнительных процедур по их стабилизации.
Всплески на графиках на рис. 8 происходят в моменты ударного касания ног об опорную поверхность после их переноса.
Для сравнения данная задача была решена с использованием классического подхода к формированию уравнений движения механических систем с дополнительными связями. Уравнения динамики ШМ строились в форме уравнений Лагранжа первого рода в обобщенных координатах [6].
Для определения множителей Лагранжа
*
Т динамические уравнения замыкались один раз продифференцированными уравнениями кинематических связей (3). Получаемая в результате система дифференциально-алгебраических уравнений (ДАЕ) относительно обобщенных ускорений и множителей Лаграгжа интегрировалась встроенной процедурой NDSolve. Как показал вычислительный эксперимент, интегрирование такой замкнутой системы уравнений движения приводит в процессе построения численного решения к значительному уходу с геометрических связей (4). Результаты моделирования представлены на рис. 9.
Рис. 7. Обобщенные координаты
Рис. 9. Уход со связи при использовании уравнений Лагранжа первого рода
Эти результаты показывают, что при замыкании уравнений движения на уровне ускорений требуется дополнительная стабилизация связей, например, методом Баумгарта [8, 9]. При этом система ДАЕ становится более жесткой, что приводит к дроблению шага численного решения и замедлению счета. В отличие от классического подхода, интегрирование уравнений (5), (6) избавляет от указанных проблем.
Заключение
В работе получена новая матричная форма уравнений движения систем связанных твердых тел с замкнутыми кинематическими цепями. Уравнения содержат импульсы Пуассона, обобщенные координаты, квазискорости и множители Лагранжа. Выписаны рекуррентные формулы, предназначенные для компьютерного формирования уравнений из простейших основных блоков, описывающих структуру, масс-инерционные, геометрические и кинематические характеристики отдельных звеньев механической системы.
Группа кинематических уравнений, входящих в расширенную систему уравнений (5), (6), имеет квазитрехдиагональную блочную структуру относительно переменных {¿¡. V. //}. что позволяет для их разрешения относительно обобщенных скоростей использовать методы "прогонки", а для определения множителей Лагранжа Т] , появляющихся в уравнениях в результате учета дополнительных связей (3), - итерационные методы. В отличие от работы [2] множители Т] вычислялись с использованием системы линейных уравнений с симметричной положительно полуопределенной матрицей системы.
Предлагаемые в статье алгоритмы позволяют сократить временные затраты на моделирование механических систем, содержащих длинные кинематические цепи, в сравнении с алгоритмами моделирования САТТ, основанными на использовании уравнений в форме Лагранжа второго или первого рода.
Приведены все этапы формирования исходной информации и результаты компьютерного моделирования одной демонстрационной механической системы описанным в работе методом. Результаты показывают, что предлагаемый метод позволяет получать удовлетворительные результаты численных расчетов стандартными методами интегрирова-
ния дифференциальных уравнений, встроенными в существующие пакеты прикладных программ, без дополнительной стабилизации уравнений связей.
Список источников
1. Иванов В.Н. Матричные уравнения движения систем твердых тел в гамильтоновых переменных. Системы со структурой дерева // Вестник Пермского университета. Математика. Механика. Информатика. 2019. Вып. 3(46). С. 38-46.
2. Иванов В.Н. Матричные уравнения движения систем твердых тел в гамильтоновых переменных. Системы с замкнутыми цепями // Вестник Пермского университета. Математика. Механика. Информатика. 2019. Вып. 4(47). С. 13-20.
3. Иванов В.Н. Алгоритмы решения уравнений движения в импульсах Пуассона систем твердых тел со структурой дерева // Вестник Пермского университета. Математика. Механика. Информатика. 2017. Вып. 4(39). С. 25-31.
4. Голуб Дж., Ван Лоун Ч. Матричные вычисления. М.: Мир, 1999. 548c.
5. Шевцов Г.С. Численные методы линейной алгебры: учеб. пособие / Г.С. Шевцов, О.Г. Крюкова, Б.И. Мызникова. М.: Финансы и статистика: ИНФРА-М, 2008. 480 с.
6. Суслов Г.К. Теоретическая механика. М.: Гостехиздат, 1946. 656с.
7. Featherstone R. Rigid Body Dynamics Algorithms. New York: Springer, 2008. 272 p.
8. Wittenburg J. Dynamics of multibody systems. Berlin: Springer-Verlag, 2008. 223 p.
9. Fan-Chung Tseng, Zheng-Dong Ma, Gregory M. Hulbert. Efficient numerical solution of constrained multibody dynamics systems // Computer Methods in Applied. Mechanics and Engineering. 192 (2003). P. 439-472.
References
1. Ivanov V.N. Matrichnye uravneniya dvizheni-ya sistem tverdyh tel v gamil'tonovyh peremennyh. Sistemy so strukturoj dereva. Vestnik Permskogo universiteta. Matematika. Mekhanika. Informatika. 2019;3(46):38-46. (In Russ.).
2. Ivanov V.N. Matrichnye uravneniya dvizheni-ya sistem tverdyh tel v gamil'tonovyh peremennyh. Sistemy s zamknutymi cepyami. Vestnik Permskogo universiteta. Matematika. Mekhanika. Informatika. 2019;4(47):13-20. (In Russ.).
3. Ivanov V.N. Algoritmy resheniya uravnenij dvizheniya v impul'sah Puassona sistem tverdyh tel so strukturoj dereva. Vestnik Permskogo universiteta. Matematika. Mekhanika. Informat-ika. 2017;4(39):25—31. (In Russ.).
4. Golub Dzh., Van Loun Ch. Matrichnye vychisleniya. M.: Mir, 1999. 548 s.
5. Shevcov G.S., O.G. Kryukova, B.I. Myznikova. Chislennye metody linejnoj algebry: ucheb. Posobie. M.: Finansy i statistika: INFRA-M, 2008. 480 s.
6. Suslov G.K. Teoreticheskaya mekhanika. M.: Gostekhizdat, 1946. 656 s.
7. Featherstone R. Rigid Body Dynamics Algorithms. New York: Springer, 2008. 272 p.
8. Wittenburg J. Dynamics of multibody systems. Berlin: Springer-Verlag, 2008. 223 p.
9. Fan-Chung Tseng, Zheng-Dong Ma, Gregory M. Hulbert. Efficient numerical solution of constrained multibody dynamics systems. Computer Methods in Applied. Mechanics and Engineering. 2003;(192):439-472.
Информация об авторе:
В. Н. Иванов - кандидат физико-математических наук, доцент, доцент кафедры высшей математики механико-математического факультета Пермского государственного национального исследовательского университета (614068, г. Пермь, ул. Букирева, 15), AuthorlD 13318.
Information about the author:
V. N. Ivanov - Candidate of Physical and Mathematical Sciences, Associate Professor, Associate Professor of the Department of Higher Mathematics, Faculty of Mechanics and Mathematics, Perm State University (Bukireva Street, 15, Perm, 614068), AuthorlD 13318.
