Научный журнал
Фундаментальные исследования
ISSN 1812-7339
"Перечень" ВАК
ИФ РИНЦ = 1,074

МЕТОД КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ ДИНАМИКИ СИСТЕМ СВЯЗАННЫХ ТВЁРДЫХ ТЕЛ

Шимановский В.А. 1
1 ФГБОУ ВО «Пермский государственный национальный исследовательский университет»
Работа посвящена разработке методов компьютерного моделирования динамики сложных механических систем. В статье представлен новый алгоритм разрешения уравнений движения систем абсолютно твердых тел со структурой дерева относительно ускорений, ориентированный на использование ЭВМ. Алгоритм основан на применении схемы Холецкого для решения системы дифференциально-алгебраических уравнений движения. Получены рекуррентные формулы для определения всех кинематических и динамических переменных, входящих в уравнения. Вычислительная сложность решения системы уравнений с помощью данного алгоритма растёт по линейному закону в зависимости от числа тел в механической системе. Проведено сравнение предлагаемого алгоритма с классическим методом А. Ф. Верещагина. На примерах интегрирования уравнений движения систем тел с большим числом степеней свободы показано преимущество разработанного алгоритма перед классическим методом.
система связанных тел
уравнения движения
обобщенные координаты
динамика
математическое моделирование
численные методы
разложение Холецкого
1. Иванов В.Н., Домбровский И.В., Набоков Ф.В. и др. Классификация моделей систем твердых тел, используемых в численных расчетах динамического поведения машиностроительных конструкций // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. – 2012. – № 2. – С. 139–155.
2. Wittenburg J. Dynamics of Multibody Systems. Berlin: Springer-Verlag, 2008. – 223 p.
3. Shabana A.A. Dynamics of Multibody Systems. Cambridge University Press, 2005. – 385 p.
4. Шимановский В.А., Иванов В.Н. Методы составления уравнений движения систем связанных твердых тел в декартовых координатах // Проблемы механики и управления: Нелинейные динамические системы. – 2007. – Вып. 39. – С. 248–262.
5. Верещагин A.Ф. Метод моделирования на ЦВМ динамики сложных механизмов роботов-манипуляторов // Изв. АН СССР. Техническая кибернетика. – 1974. – № 6. – С. 89–94.
6. Golub G.H., Van Loan C.F. Matrix Computations // The Johns Hopkins University Press. – 2012. – 790 p.
7. Hairer E., Norsett S.P., Wanner G. Solving ordinary differential equations I. Nonstiff Problems // Springer-Verlag. – 2011. – 528 p.

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

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

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

Уравнения движения системы связанных твёрдых тел

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

Пусть N – число тел в системе (не считая тела «0», движение которого относительно инерциальной системы координат Oxyz задано функцией времени). Остальным телам системы присвоим номера таким образом, чтобы для любого тела номер предшествующего ему тела был меньше. Для полного описания структуры взаимосвязей в системе достаточно одного целочисленного массива himan01.wmf, на i-м месте которого расположен индекс тела, предшествующего i-му. Массив k позволяет определить для каждого тела системы целочисленные массивы Si номеров тел, для которых i-е тело является предшествующим.

Обозначим через Oi полюс ортогональной декартовой системы координат, связанной с телом. Положение i-го тела в пространстве относительно базового для него ki-го тела полностью определяется матрицей-столбцом ρi, элементами которой являются проекции вектора himan02.wmf на оси системы координат, связанной с телом ki и матрицей Gi направляющих косинусов между базисными векторами, связанными с телами ki и i. Абсолютную скорость каждого тела зададим 6-мерной матрицей-столбцом νi = (ωi, νi)T, где ωi, νi – векторы мгновенных угловой и линейной скоростей i-го тела. Аналогично, абсолютное ускорение каждого тела системы определим с помощью 6-мерной матрицы-столбца wi = (εi, wi)T, где εi, wi – мгновенные угловое и линейное ускорения i-го тела. Векторы νi, ωi, wi и εi будем задавать в системе координат, связанной с телом.

Для каждого соединения тел в соответствии с видом уравнений связи введём матрицу-столбец локальных обобщённых координат himan03.wmf (ni ≤ 6). Тогда элементы матриц ρi и Gi будут являться в общем случае функциями времени и выбранных координат относительного движения. В качестве обобщённых параметров системы примем совокупность himan04.wmf локальных обобщённых координат himan05.wmf.

Принимая движение i-го тела за относительное, а предшествующего ему ki-го тела за переносное, можно записать рекуррентные формулы для вычисления скоростей и ускорений тел механической системы [1]:

himan06.wmf (1)

himan07.wmf (2)

где матрицы Ci имеют вид

himan08.wmf

Ai – матрицы локального касательного базиса многообразия относительного движения в i-м сочленении; himan09.wmf – матрицы-столбцы составляющих относительных скоростей тел системы, не зависящих от обобщённых скоростей himan10.wmf; а в матрицы-столбцы himan11.wmf вошли все члены переносного, кориолисова и относительного ускорений. Здесь и далее символ «∼» используется для обозначения кососимметричной матрицы [2].

Обозначим через fi и μi главный вектор и главный момент активных сил, действующих на i-е тело. Уравнения движения выпишем в форме уравнений Ньютона – Эйлера для свободного твёрдого тела [3]:

himan12.wmf (3)

со следующими матрицами

himan13.wmf

где mi – масса i-го тела, Ji – тензор инерции i-го тела в связанной с телом системе координат, himan14.wmf – радиус-вектор центра масс i-го тела в связанной с ним системе координат, Ri – матрица-столбец реакции связей в i-м сочленении, himan15.wmf.

Уравнения (2) и (3) замыкаются соотношениями

himan16.wmf (4)

отражающими то обстоятельство, что при идеальности связей в сочленениях реакции являются ортогональными к конфигурационному многообразию относительных положений, определяемому параметрами qi. Таким образом, получается замкнутая система 12N + n скалярных дифференциальных и алгебраических уравнений (2)–(4) относительно такого же количества неизвестных величин wi, Ri и himan17.wmf.

При интегрировании уравнений движения на ЭВМ численными методами необходимо на каждом шаге один или несколько раз вычислять ускорения himan18.wmf при заданных координатах q и скоростях himan19.wmf. В нашем случае определение ускорений himan20.wmf сводится к простой алгебраической задаче решения системы линейных уравнений (2)–(4) с учётом соотношений (1).

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

himan21.wmf (5)

где

himan22.wmf

himan23.wmf

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

Разрешение уравнений движения методом прогонки

Метод прогонки тел по существу является модификацией метода Гаусса решения систем линейных алгебраических уравнений с ленточной структурой. В данном методе при прямом ходе, который выполняется начиная с последнего тела системы, из группы уравнений (3) исключаются реакции связей носимых тел. При обратном ходе по явным формулам вычисляются последовательно для каждого тела системы обобщённые ускорения himan24.wmf, декартовые ускорения wi и реакции связей Ri.

Суть метода подробно изложена в статье [4]. Здесь же приведём только алгоритм метода прогонки, обобщённый на случай систем тел со структурой дерева:

for i = N:1

himan25.wmf

himan26.wmf

end

for i = 1: N

himan27.wmf

himan28.wmf

himan29.wmf

end

Трудоёмкость вычислений по указанным рекуррентным формулам линейно растёт с ростом числа тел в системе. При этом в алгоритме требуется обращение симметричных положительно определённых матриц himan30.wmf, порядок которых равен числу степеней свободы в i-м сочленении, т.е. не превышает шести. Этим и обуславливается высокая эффективность описанного выше алгоритма разрешения уравнений движения относительно ускорений.

Впервые подобные формулы были получены А.Ф. Верещагиным в 1974 г. [5] для пространственных механизмов, содержащих ряд звеньев, соединённых в цепь с одной степенью подвижности между звеньями.

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

Идея исключения Гаусса – это преобразование исходной системы уравнений в эквивалентную треугольную систему с последующим её решением. В матричном виде это есть LU-разложение [6]. В случае симметричной матрицы системы уравнений (5) лучше использовать LTDL-разложение. Заметим, что при разложении ленточных систем треугольные множители L также являются ленточными.

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

himan71.wmf

со следующими матрицами

himan33.wmf

himan34.wmf

где himan35.wmf и himan36.wmf – факторы Холецкого симметричных положительно определённых матриц VN и himan37.wmf соответственно, матрица VN находится из матричного уравнения himan38.wmf,
а матрица himan39.wmf вычисляется по формуле himan40.wmf.

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

himan43.wmf

где

himan44.wmf

himan45.wmf

himan46.wmf

himan47.wmf

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

himan48.wmf (6)

himan49.wmf (7)

где

himan50.wmf, himan51.wmf,

himan52.wmf,

himan53.wmf.

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

himan54.wmf

himan55.wmf

himan56.wmf

где himan57.wmf.

Решая систему (7) и исключая переменные zi, yi, получим

himan58.wmf

himan59.wmf

himan60.wmf

Обобщая полученные формулы на случай системы твёрдых тел со структурой дерева, можно записать следующий алгоритм решения расширенной системы уравнений (2)–(4):

for i = N:1

himan61.wmf

himan62.wmf

Найти фактор Холецкого himan63.wmf матрицы himan64.wmf

himan65.wmf

himan66.wmf

end

for i = 1: N

himan67.wmf

himan68.wmf

himan69.wmf

end

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

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

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

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

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

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

himan1a.wmf himan1b.wmf

а) б)

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

himan2.wmf

Рис. 2. График отношения T1/T2

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

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

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

Заключение

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


Библиографическая ссылка

Шимановский В.А. МЕТОД КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ ДИНАМИКИ СИСТЕМ СВЯЗАННЫХ ТВЁРДЫХ ТЕЛ // Фундаментальные исследования. – 2017. – № 8-1. – С. 104-109;
URL: http://www.fundamental-research.ru/ru/article/view?id=41629 (дата обращения: 09.12.2019).

Предлагаем вашему вниманию журналы, издающиеся в издательстве «Академия Естествознания»
(Высокий импакт-фактор РИНЦ, тематика журналов охватывает все научные направления)

«Фундаментальные исследования» список ВАК ИФ РИНЦ = 1.074