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

2015 4(31)

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

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

Вып.

МЕХАНИКА МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 531.01+004.94

Численные методы исследования механических систем с дополнительными связями

В. Н. Иванов

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

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

Ключевые слова механические системы; голономные связи; уравнения движения; численное интегрирование; модифицированные функции Лагранжа; математическое моделирование.

Введение

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

[1-15].

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

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

"MSC.ADAMS" [16], "Универсальный меха-низм-иМ" [17], "ФРУНД" [18], "Е^Е^" [19] и другие.

Универсализм перечисленных комплексов программ имеет отрицательную сторону -это увеличенные временные затраты при проведении расчетов и изменении математической модели.

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

На ранних стадиях проектирования в качестве основной расчетной схемы при создании математической модели обычно ис-

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

Известно, что математическая модель или система уравнений движения любой связки твердых тел является системой линейных дифференциально-алгебраических уравнений относительно различных групп переменных: обобщенных или декартовых ускорений тел системы, множителей Лагранжа, реакций связей, импульсов. Матрица системы уравнений в общем случае зависит от обобщенных координат и является переменной во времени. На каждом шаге численного интегрирования необходимо приводить уравнения движения к явному виду, т.е. разрешать относительно указанных групп переменных. Это требует определенных вычислительных затрат, объем которых зависит от выбранного метода решения, плотности заполнения матрицы системы и ее структуры.

Математические модели механических систем со структурой дерева имеют рекуррентную структуру, и для них в настоящее время разработаны достаточно эффективные методы решения: метод составных тел для уравнений движения в форме Лагранжа II рода [4, 7, 15], методы «прогонки» (отдельных тел) для расширенной системы общих уравнений динамики в декартовых координатах и уравнений движения в форме Лагранжа I рода [8, 9, 11, 15], итерационные методы решения линейных систем с плотно заполненной или разреженной матрицей системы [2, 3, 5].

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

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

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

1. Постановка задачи

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

Мд = Q (д, д, г), (1)

где

д = д (t)е Rп - вектор обобщенных (абсолютных) координат, dд

д =--вектор скоростей,

dt

.. (12а

д = —- - вектор ускорений,

Л 2

М = М (д)е Япхп - матрица инерции, Q = Q (д, д, t)е Яп - вектор внешних и

внутренних активных сил и сил инерции, п - число степеней свободы или число уравнений (1).

Пусть на систему наложено т дополнительных линейно-независимых голономных, склерономных связей (т < п) :

g(д) = о, (2)

где g е Rm - дважды дифференцируемая вектор-функция.

Требуется найти частное решение системы уравнений (1) при заданных дополнительных ограничениях - связях (2) с начальными условиями: q (0) = q0 , q (0) = q0 .

2. Методы учета связей в механических системах

В настоящее время для учета дополнительных связей в механических системах используются следующие методы:

- классический метод Лагранжа [10, 20, 21, 22];

- метод стабилизации связей Баумгарта [22-26];

- методы замены связей упруго-демпфирующими элементами (УДЭ) [12, 15, 22];

- методы учета связей с помощью штрафных функций (ШФ) [25, 26];

- численные методы прямого интегрирования дифференциально-алгебраических систем уравнений (DAE) [14, 25-30, 32];

- методы разделения переменных состояния - координат q (t) на зависимые и независимые и проектирования уравнений на подпространства независимых координат [21, 25, 26, 31, 33].

Подробные обзоры по методам учета связей в механических системах можно найти в [25, 26].

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

Классический метод Лагранжа. Дополнительные связи (2) учитываются в уравнениях (1) с помощью множителей Лагранжа:

Mq + GTX = Q, (3)

где Л е Rm - вектор множителей Лагранжа, R = -GT Л - вектор обобщенных реакций связей.

Для замыкания динамических уравнений связи (2) дважды дифференцируются:

g = Gq = 0, g = Gq + (Gq = 0, (4)

G^d rtmxm г

= —— е r - матрица базиса орто-dq

гонального многообразия к дополнительным связям.

В результате получается замкнутая система уравнений (3)-(4) относительно пары векторов (q,X) .

Если матрица G имеет полный ранг m, то систему (3)-(4) легко разрешить относительно ускорений.

Из (3) выражаем ускорения:

q = M- (Q - GT Л).

Далее подставляем в (4), находим множители Лагранжа:

л = (GM~lGT )-1 (GM_1Q + Gq). (5) Подставляем их в (3), находим ускорения:

q = ^ E - GT (GM~1GT )-1 G jM_1Q -

-m ~1GT (gm ~1GT )-1 Gq.

Здесь Z = E - GT (GM~1GT )-1 GM

матрица проектирования «-мерных векторов в подпространство, касательное к связям (2).

Недостатки метода Лагранжа: требуется решать дополнительные СЛАУ, возможен уход со связей, так как дифференциальные уравнения (4) имеют общее решение g (t ) = c1t + c2, т. е. со временем (при наличии

ошибок вычислений) происходит сход со связи.

Для стабилизации связей используют два основных подхода:

1. Метод Баумгарта [22-26].

2. Численные методы прямого интегрирования дифференциально-алгебраических систем уравнений (DAE), в которых происходит проектирование численного решения ДУ (3), (4) на алгебраические связи g = 0 и g = 0 [14, 25-30, 32].

Метод Баумгарта. Баумгарт предложил заменить дифференциальные уравнения (ДУ) связей (4) другими ДУ, которые имеют то же тривиальное частное решение вида (2) g (t) = 0, но асимптотически устойчивое в

смысле Ляпунова. Основная идея стабилизации связей - демпфирование ускорений связей g = 0 с помощью обратных связей по на-

рушениям ограничений на перемещения g ^) и скорости g ^):

g + 2Dg + Cg = 0 или в развернутом виде:

Gq = -2DGq - Cg - Gq .

(7)

где D, С - положительные диагональные матрицы параметров, С - матрица коэффициентов жесткости - круговые частоты), 2В - матрица коэффициентов демпфирования колебаний невязок связей g ^) вокруг

тривиального решения (логарифмических декрементов).

Объединяя уравнения (3) и (7) получим замкнутую систему уравнений относительно пары векторов (д,Х) :

Мд + ОтХ = Q, (8)

од = -2Вод - Cg - од. (9)

Введение упругих Cg и демпфирующих членов 2ВОд в уравнения ускорений связей позволяет исключить рост нарушения ограничений в процессе их интегрирования и сделать решение g (t ) = 0 асимптотически устойчивым.

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

C = OI-2 h

= O

(

g

\

ся паразитные вариации траектории, искажаются реакции связей.

Численные методы прямого интегрирования DAE. При этом подходе методы численного интегрирования DAE применяются к расширенной системе уравнений движения механических систем, которая включает в себя динамические уравнения (3), уравнения алгебраических связей (2), кинематических связей g = 0 и связей на ускорения g = 0:

Mq + GTX = Q, (10)

g(q) = 0, (11)

g = Gq = 0, (12)

g = Gq + Gq = 0. (13)

В методах интегрирования DAE обычно используется метод Ньютона решения нелинейных уравнений (11) для проектирования численного решения системы ДУ (10) и (13) на алгебраические связи (11) и (12).

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

M(q)qu + GT (q)a = Q(q,q,t), (14) G ( q ) qu = -2DG ( q ) q - Cg ( q )-G ( q ) q, (15)

M (q - qu ) = GT (qu ) ц, G (q - qu ) = - g (qu).

(16) (17)

В = О (л/с ),

где И - шаг интегрирования, ®тах- наивысшая частота колебаний (собственных или вынужденных), учитываемая при моделировании. Коэффициенты демпфирования В выбираются из условия отсутствия колебаний функции g ^) , характеризующей нарушение связи, вокруг стационарного решения g (t ) = 0.

Недостатки метода Баумгарта: требуется решать дополнительные системы линейных уравнений, повышается жесткость системы, уменьшается шаг интегрирования, появляют-

Уравнения (14), (15) являются системой DAE относительно промежуточных переменных qu и Л. Уравнения (16), (17) - это система алгебраических уравнений относительно переменных q и ц . Они задают итерационную процедуру метода Ньютона проектирования решения q (t) на связи g (q) = 0.

Обычно ограничиваются одним шагом алгоритма. Заметим, что эффективность методов DAE возрастает, если они сочетаются с методом Баумгарта, как в системе уравнений (14)-(17).

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

Методы замены связей упруго-демпфирующими элементами (УДЭ). Связи

> ютах или

разрезаются и заменяются линейными упруго-демпфирующими элементами:

Мд = Q - GT С + 2Dg),

(18)

СП отхт

, d е к - диагональные матрицы жесткости и демпфирования в разрезанных шарнирах.

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

С точки зрения методов оптимизации этот метод есть метод внешнего квадратичного штрафа в задачах условной минимизации [32, 33].

Недостатки метода УДЭ: минимум энергии достигается не в точке g(д) = 0, а в

точке X = Сg + 2Dg , где X - вектор множителей Лагранжа в уравнениях Лагранжа (3), (4). Удовлетворительная точность достигается, только если элементы матриц С, D ^да, появляются дополнительные высокочастотные малые колебания в системе, замедляется численный расчет.

Методы учета связей с помощью штрафных функций (ШФ). При выводе уравнений движения из интегральных принципов механики ограничения (11), (12) и (13) будем учитывать не методом Лагранжа, а методом внешнего штрафа [34, 35].

Из принципа наименьшего действия Гамильтона-Остроградского:

t+h

SS = $ 5цт (Q - Мд) -аЗС + 2Dg + £)) dt =

t

t +h

= | Sqт (0 -Мд - аОт С + 2Dg + g))dt = 0,

где а - штрафной коэффициент,

СП тлту-т

, d е к - диагональные матрицы жесткости и демпфирования, получаем уравнения:

(М + аGTG ) д--

(19)

= Q -аОт С + 2DGq + Gq).

В уравнениях (19), в отличие от уравнений (18), элементы матриц С и D можно

брать небольшими. Однако этот метод имеет те же недостатки, что метод УДЭ.

Недостатки метода ШФ: минимум энергии достигается в точке X = а (Сg + Dg + £),

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

Упомянутые выше методы учета дополнительных связей с помощью разделения переменных состояния - координат д ^) на зависимые и независимые и проектирования уравнений движения на подпространства независимых координат являются традиционными в классической механике. Их преимущества и недостатки хорошо известны. Практические рекомендации по их применению можно найти в [10, 21, 25, 26, 31, 33]. Заметим только, что проекционные методы трудно применять для автоматизации этапа формирования математической модели механической системы и поэтому они не рассматриваются в этой статье.

3. Метод модифицированных функций Лагранжа

При выводе уравнений движения из принципа Гамильтона-Остроградского дополнительные связи можно учитывать методом модифицированных функций Лагранжа в форме Пауэлла [34, 35]. Тогда уравнения движения преобразуются к виду:

Мд = 0 - Gт (а (£ + е)),

где е е Rт - вектор модифицированных множителей Лагранжа, а - вектор штрафных коэффициентов. Заметим, что в точке минимума действия по Гамильтону существует связь между обычными и модифицированными множителями Лагранжа: X = ае*. Это означает, что для определения модифицированных множителей Лагранжа уравнения движения можно замкнуть, как и в классическом варианте, соотношениями (4). При этом для стабилизации связей можно так же использовать метод Баумгарта. Таким образом, в окончательном виде замкнутая система динамических уравнений имеет вид:

Мд = 0 - GT (а(£ + е)), (20)

од = -2Вод - Cg - Од.

(21)

С другой стороны, для пересчета переменных е можно использовать различные приближенные итерационные формулы. Например, формулу с линейной скоростью сходимости: е = ^ (д) + ^ (д), где К, К -коэффициенты усиления обратной связи соответственно по нарушению и ускорению нарушения связей.

Если ввести новые переменные

и = g+К | gdt■.

то уравнения движения пре-

образуются к виду:

мд = Q - от (Аи+Bg), и=^ (д)+g (д )=^ (д)+Од,

где А = aku, В = akw .

(22) (23)

Уравнения движения в форме (22), (23) можно рассматривать как уравнения системы автоматического регулирования, где управление берется в форме ПИД-регулятора [36]. Физически это означает, что к системе добавляются некоторые силовые элементы с управляемыми параметрами, которые за счет расширения фазового пространства позволяют обеспечить заданную точность моделирования.

Выбор параметров ПИД-регулятора. Настройка ПИД-регулятора в уравнениях (22), (23) заключается в подборе 3 векторных параметров: ^, С, В. Можно использовать стандартные методы из теории автоматического регулирования, например метод Зигле-ра-Никольса [36].

Однако для поставленной задачи более пригодным оказывается другой подход. Пусть на механическую систему накладывается только одна дополнительная связь (2). Если из системы уравнений (22), (23) исключить координаты механической системы, убрать внешнее силовое воздействие и линеаризовать, то получим линейное однородное ДУ, описывающее собственные колебания модифицированного множителя Лагранжа е:

е• + ОМ~1От (Ве + Ае + ^Ае) = 0. (24)

Потребуем, чтобы характеристическое уравнение для ДУ (24) имело вид:

^ + а0)(k + (а + ш))(к + (а -с)) = 0, (25)

где а0, А\ - коэффициенты демпфирования, ш - круговая частота регулятора (23).

Полагаем: со > штах , /30 = О (ра),

а = о (рс), р = о (1).

Сопоставляя (24) и (25), получаем связь между коэффициентами

ом ~1отв = а0 + 2а, ОМ ~1От А = 2£„а + а2 + с2,

ом-1очиА = а0 (а2 + с2).

Из установленных соотношений определяем следующие параметры ПИД-регулятора:

В = ( ом ~1от ) 3 рс, А = (оМ _1от )-1 (3 р 2 +1) К = с( р3 + р )/(3 р2 +1).

22 с ,

(26)

Равенства (26) задают коэффициенты усиления обратной связи, демпфирующие и упругие параметры ПИД-регулятора (23).

С практической точки зрения преимущество уравнений (22), (23) заключается в том, что система уравнений (22) имеет ту же структуру, что и исходная система (1) без дополнительных связей. Это означает, что для ее разрешения относительно ускорений можно использовать эффективные методы решения рекуррентных систем уравнений, если исходная механическая система имеет структуру дерева. Уравнения (22), (23) построены таким образом, что при возникновении отклонений Ag (д) изменяется в первую очередь ненапряженная длина дополнительного упруго-демпфирующего элемента Аи + Вg так, чтобы создаваемая этим элементом дополнительная сила обеспечивала скольжение механической

системы по дополнительной связи

g (д)«0.

4. Вычислительные эксперименты

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

Сравнивались характеристики следующих методов:

1. Метод Лагранжа, уравнения (3), (4).

2. Метод Баумгарта, уравнения (8), (9).

3. Метод замены связей упруго-демпфирующими элементами (УДЭ), уравнения (18).

4. Метод модифицированных функций Лагранжа и стабилизации связей по Баумгар-ту (МФЛ+Б), уравнения (20), (21).

5. Метод модифицированных функций Лагранжа с управляющими функциями (МФЛ+У), уравнения (22), (23).

Вычисления проводились в среде САВ Mathematica. Интегрирование выполнялось встроенной процедурой NDSolve с установленными по умолчанию параметрами. Заметим, что во всех случая из-за наличия множителей Лагранжа процедура NDSolve интегрировала предлагаемые уравнения движения как дифференциально-алгебраические системы уравнений. Поэтому в расчетах отдельно не выделялся алгоритм (14)-(17) интегрирования DAE.

Первая система представляет собой два многозвенных математических маятника

Mi...Mn и Mn+i.

M

2n :

связанных невесо-

мыми нерастяжимыми балками М2кМ4к,

п

М2кМ4к -1, к = 1,2,...,— и упруго-демфирую-

щими связями между массами М2к-1 и М4к-1 (рис. 1). В результате получается система с замкнутыми кинематическими цепями.

Рис. 1. Многозвенный маятник

Для каждого маятника можно записать: геометрические соотношения

Xj = xi_1 +1 cos qi, yi = y_i +1sin qi, Уп+1 = h +l sin qn+i и кинематические соотношения *i = Xj_1 _ qjl sin qi, yi = yj_i + qil cos q, где l = MiMi+1 - длина каждого звена маятников, h = MtM 2i - расстояние между двумя цепочками маятников, qt - обобщенные координаты - абсолютные углы отклонения звеньев от вертикали. Результаты расчетов приводятся для случая n = 6, l = 1.8 м , h = 4 м, mt = 10 кг .

Уравнения движения для каждого маятника строились в форме Лагранжа II рода:

d

( дК^

dql

д( К _ П)

dqi

= Qi, i = 1,2n,

где К, П - кинетическая и потенциальная энергии:

2n 1

К = 1 ± m (xf + y2 ) , i=1 2

n

П=Z mg ((il _ xi)+(il _ xi+n)). i=1

Дополнительные связи:

gk (q) = (X2k _ X4k )2 + (У2k _ У4k )2 _12 = 0

gk+n/2 (q) =(X2k _X4k_1 )2 +(У2к _У4к_1)2 _l2 _h2 = 0

k = 1,2,..., n¡ 2.

Обобщенные силы демпфирующих связей:

от упруго-

0 = -22 (сРк (д)+dpk (д)) дРк (д), где

к=1 д1

Рк (д) = (Х2к-1 - Х4к-1 )2 + (У2к-1 - У4к-1 )2 -12 =

В расчетах принимались следующие коэффициенты жесткости и демпфирования: с = 50000 н / м, d = 100 н ■ с / м .

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

На рис. 2 приведены графики колебаний всех обобщенных координат (результаты

по всем методам).

Рис. 2. Колебания всех звеньев

Следующие рисунки (рис. 3-7) представляют графики невязок дополнительных связей у = gk (д), k = 1,2,3, получаемых раз-

личными методами.

Рис. 4. Метод Баумгарта

Рис. 5. Метод УДЭ

Графики (рис. 3-7) показывают, что метод Лагранжа приводит к постепенному сходу со связей. Наименее точным является метод замены связей упруго-демпфирующими элементами. Видно, что он плохо справляется с компенсированием отклонений в положении звеньев системы, вызванных действием силы тяжести. Лучшие результаты показывают методы модифицированных функций Лагранжа.

В табл. 1 представлено время, затраченное на моделирование различными методами.

Таблица 1

Метод Время, с

Лагранжа 2.83

Баумгарта 3.78

УДЕ 13.63

МФЛ+Б 2.19

МФЛ+У 5.77

Из таблицы следует, что наиболее быстрым является метод МФЛ+Б, а наиболее медленным - метод УДЭ.

В качестве второй системы выбрана квадратная решетка, состоящая из 40 одинаковых твердых тел в форме тонких квадратных параллелепипедов (рис. 8). На рисунке для наглядности изображена часть решетки, содержащая 21 тело). Ширина каждого параллелепипеда равна а = 1/7 м, толщина -И = 0.01 м, масса - т = 5/4 кг. Каждое тело системы связывается с соседними телами шаровыми шарнирами, расположенными в центрах боковых граней тел. В центре каждой группы из 9 звеньев решетки отсутствует одно

тело, что обеспечивает статическую определенность всей системы. Первое тело системы закреплено в абсолютной системе координат шаровым шарниром в точке на первом теле с координатами (-а/2, -а/20, 0). В начальном положении решетка располагается в первой четверти горизонтальной плоскости OXY. Ось 02 направлена вниз, вдоль направления действия силы тяжести. Данная конструкция взята из работы [37]. Ранее она использовалась для сравнения эффективности различных программных комплексов ("Универсальный механизм-иМ" [16], "ФРУНД" [17], "Е^Е^" [18]), предназначенных для моделирования систем твердых тел.

Для обеспечения жесткости всей конструкции в шаровых шарнирах, связывающих тела системы, предусмотрены линейные упруго-демпфирующие моменты с коэффициентами жесткости с = 100 н ■ м / рад и демпфирования d = 0.01 н ■ м ■ с / рад . Время моделирования Т = 4 с.

Рис. 8. Решетка твердых тел

При отсутствии дополнительных связей каждое тело системы имеет 6 степеней свободы, 3 линейных - абсолютных координат центров тяжести тел г =(х{,у1,zj') и 3 угловых -

углов Эйлера дг = (ц>г,0,).

Уравнения движения для каждого маятника строились в форме Лагранжа II рода:

d

d

( дКЛ

дГ

дК дд1

\

д(К - П)

дг

д(К - П)

дд,

= &, , = 1,40,

= 0,, i = 1,40,

где К, П - кинетическая и потенциальная энергии:

К = 2X(тГ2 + , 2 ,=1

40

П=х ,

,=1

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

Уравнения дополнительных связей в шаровых шарнирах:

= Г +1} ■(а/2,0,0)т -Гк + 1к ■(-а/2,0,0)т ,

= г} +1} ■( 0, а/2,0)т - Гк + 1к ■( 0, -а/2,0)т ,

= г1 +11 ■ (-а /2, -а / 20,0)т ,

где I^ - матрицы преобразования координат

из систем координат, связанных с телами в абсолютную систему координат. Индексы ] и к перебираются по всем номерам соседних тел в решетке.

Обобщенные силы от упруго-демпфирующих связей:

е,=-я ср,к (д)+чРк (д) )дрдд,

] ,к д,

где р}к (д) = д} - дк.

На рис. 9-10 приведены графики колебаний всех обобщенных координат (результаты по всем методам).

Рис. 9. Линейные колебания всех звеньев

Рис. 10. Угловые колебания всех звеньев

Следующие рисунки (рис. 11-15) представляют графики невязок дополнительных

связей y = gk (q), k = 1,2,3, получаемых различными методами.

i. X Ii'

4. х 10

J. х ID'

-J. x 10'

5. xlO'

-J. x 10'

-1. xlO

-1.5 xlO"

Рис. 12. Метод Баумгарта

Рис. 13. Метод УДЭ

-1. xlO'

Рис. 14. Метод МФЛ+Б

Графики (рис. 11-15) показывают, что метод Лагранжа приводит к постепенному сходу со связей, наименее точным является метод замены связей упруго-демпфирующими элементами. Лучшие результаты показывают методы модифицированных функций Лагранжа.

В табл. 2 представлено время, потраченное на моделирование различными методами.

Таблица 2

Метод Время, с

Лагранжа 17.02

Баумгарта 19.77

УДЕ 72.83

МФЛ+Б 14.17

МФЛ+У 56.39

Рис. 15. Метод МФЛ+ У

Из табл. 2 следует, что наиболее быстрым является метод МФЛ+Б, а наиболее медленным - метод УДЭ. Выводы

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

На примерах показано, что метод имеет преимущества перед существующими методами учета дополнительных связей.

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

1. Иванов В.Н. Применение метода модифицированных функций Лагранжа для учета дополнительных связей в механических системах // Вестник Пермского университета. Математика. Механика. Информатика. 2013. Вып. 4 (23). С. 19-28.

2. Иванов В.Н., Шимановский В.А. Использование итерационных алгоритмов разрешения уравнений движения механических систем при их численном интегрировании // Вестник Пермского университета. Математика. Механика. Информатика. 2006. Вып. 4 (4). С. 28-38.

3. Иванов В.Н., Шимановский В.А. Применение итерационных методов для разрешения уравнений движения систем связанных твердых тел // Вестник Пермского универ-

ситета. Математика. Механика. Информатика. 2008. Вып. 4 (20). С. 109-116.

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

5. Иванов В.Н. Основные свойства обратного итерационного алгоритма решения систем линейных уравнений с положительно определенными матрицами // Вычислительные методы и программирование: Новые вычислительные технологии. 2012. Т. 13. С. 366-376. URL: http://num-meth.srcc. msu.ru/ (дата обращения: 30.09.2015).

6. Бячков А.Б., Иванов В.Н., Шимановский

B.А. Классификация форм уравнений динамики систем твердых тел со структурой дерева // Вестник Пермского университета. Математика. Механика. Информатика. 2009. Вып. 7(33). С. 21-25.

7. Шимановский В.А., Иванов В.Н. Формирование уравнений движения механических систем в обобщенных координатах // Проблемы механики и управления: Нелинейные динамические системы: межвуз. сб. науч. тр. / Перм. ун-т. Пермь. 2005. Вып. 37. С. 188-201.

8. Шимановский В.А., Иванов В.Н. Методы составления уравнений движения систем связанных твердых тел в декартовых координатах // Проблемы механики и управления: Нелинейные динамические системы: межвуз. сб. науч. тр. / Перм. ун-т. Пермь. 2007. Вып. 39. С. 248-262.

9. Шимановский В.А., Иванов В.Н. Уравнения движения систем связанных твердых тел в канонических переменных // Вестник Пермского университета. Математика. Механика. Информатика. 2013. Вып. 2 (21).

C. 76-82.

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

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

12. Виттенбург Й. Динамика систем твердых тел. М.: Мир. 1980.

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

14. Погорелое Д.Ю. Введение в моделирование динамики систем тел. Брянск: БГТУ, 1997.

15. Погорелое Д.Ю. Алгоритмы синтеза и численного интегрирования уравнений движения систем тел с большим числом степеней свободы // VIII Всеросс. съезд по теоретической и прикладной механике: ан-нотац. докл. Пермь, 2001. С. 490-491.

16. Adams: система виртуального моделирования машин и механизмов // ООО "Эм-Эс-Си Софтвэр РУС", 2001-2012. URL: http://www.mscsoftware.ru/products/adams (дата обращения: 30.09.2015).

17. Универсальный механизм: динамика машин и механизмов, динамика автомобилей и железнодорожных экипажей, прикладная механика, кинематика, обратная кинематика // Лаборатория вычислительной механики / Брянский государственный технический университет. Брянск, 2012. URL: http://www.umlab.ru (дата обращения: 30.09.2015).

18. ФРУНД: моделирование динамики систем твердых и упругих тел // Волгоградский государственный технический университет. Волгоград, 2005. URL: http://frund.vstu.ru /frund.htm (дата обращения: 30.09.2015).

19. EULER: программный комплекс автоматизированного динамического анализа многокомпонентных механических систем // ЗАО "АвтоМеханика". М. 1993-2011. URL: http://www.euler.ru (дата обращения: 30.09.2015).

20. Суслов Г. К. Теоретическая механика. М.: Гостехиздат. 1946.

21. Shabana A.A. Computational Dynamics. New York: Wiley. 2001.

22. Wittenburg J. Dymamics of Multibody Systems. Berlin: Springer-Verlag. 2008.

23. Nikravesh P.E. Some Methods for Dynamic Analysis of Constrained Mechanical Systems: a Survey // NATO ASI Series: Computer Aided Analysis and Optimization of Mechanical System Dynamics. 1984. Vol. F9. Р. 351368.

24. Flores P., Ambrosio J. Revolute joints with clearance in multibody systems // Computers and Structures. 2004. Vol. 82. P. 1359-1369.

25. Tseng F.-C., Ma Z.-D., Hulbert G.M. Efficient numerical solution of constrained multi-body dynamics systems // Computer Methods in Applied Mechanics and Engineering. 2003. Vol. 192. P. 439-472.

26. Jalon J. G., Bayo E. Kinematic and Dynamic Simulation of Multibody Systems. The RealTime Challenge. Springer-Verlag, New York. 1994.

27. Terze Z., Lefeber D., Muftic O. Null Space Integration Method for Constrained Multibody Systems with No Constraint Violation // Multibody System Dynamics. 2001. Vol. 6. P. 229-243.

28. Betsch P. Energy-consistent numerical integration of mechanical systems with mixed ho-lonomic and nonholonomic constraints // Computer Methods in Applied Mechanics and Engineering. 2006. V. 195. P. 7020-7035.

29. Betsch P. The discrete null space method for the energy consistent integration of constrained mechanical systems Part I: Holonomic constraints // Computer Methods in Applied Mechanics and Engineering. 2005. Vol.194. P.5159-5190.

30. ArnoldM., Fuchs A., Fuhrer C. Efficient corrector iteration for DAE time integration in

multibody dynamics // Computer Methods in Applied Mechanics and Engineering. 2006. Vol. 195. P. 6958-6973.

31. Blajer W. Geometrical Interpretation of Multibody Dynamics: Theory and Implementations // NATO ASI Series II: Mathematics, Physics and Chemistry / Virtual Nonlinear Multibody Systems. 2002. Vol. 103. Р. 17-36.

32.Mukharlyamov R.G., Beshaw A.W. Solving differential equations of motion for constrained mechanical systems // Вестник РУДН. Серия математика, информатика, физика. № 3. 2013.С. 81-91

33. Yu Q. Chen I.-M. A Direct Violation Correction Method in Numerical Simulation of Constrained Multibody Systems // Computational Mechanics. 2000. Vol. 26. P. 52-57.

34. Гилл Ф., Мюррей У., Райт М. Практическая оптимизация. М.: Мир. 1985.

35. Nocedal J., Wright S.J. Numerical Optimization. Berlin: Springer. 2006.

36. Ощепков А.Ю. Системы автоматического управления: теория, моделирование в МАТLAB: учеб. пособие. СПб.: Лань. 2013.

37. Boykov V., Gorobtsov A., Pogorelov D. Benchmarks for Testing MBS Software Efficiency // EUROMECH Colloquium 495 / Advances in simulation of multibody system dynamics. 2008. P. 19-20.

Numerical methods for the study of mechanical systems with additional constraints

V. N. Ivanov

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

The results of research of the use of modified Lagrange functions, developed in numerical optimization methods, to account for the additional constraints in the holonomic mechanical systems are presented. The problems of estimating the parameters of the modified Lagrange functions that ensure the movement of mechanical systems along the additional constraints are investigated. This methodology allows to apply of numerical simulation algorithms, developed for mechanical systems with the structure of a tree, on a system with closed kinematic chains. Comparative efficiency of the method is shown by test examples.

Key words: mechanical systems; holonomic constraints, the equations of motion, numerical integration, the modified Lagrangian function, math modeling.