В опытно-конструкторских работах при разработке новых изделий машиностроения широко используется компьютерное моделирование. На ранних стадиях такого проектирования исследуемую механическую систему обычно идеализируют системой связанных абсолютно твердых тел. Требование точности компьютерного моделирования приводит к необходимости увеличивать число тел, на которые разбивается механическая система. C ростом размерности математической модели увеличивается и трудоемкость моделирования. Поэтому разработка методов, позволяющих ускорить процесс математического моделирования, является актуальной задачей.
Известно, что система уравнений движения связки твердых тел является системой линейных дифференциально-алгебраических уравнений (ДАУ) относительно различных групп переменных: обобщенных или декартовых ускорений тел системы, множителей Лагранжа, реакций связей, импульсов. Элементы матрицы системы уравнений в общем случае зависят от обобщенных координат и времени. На каждом шаге численного интегрирования необходимо разрешать эти уравнения относительно указанных групп переменных, что требует определенных вычислительных затрат, объем которых зависит от выбранного метода решения, количества ненулевых элементов матрицы системы и ее структуры в целом.
Математические модели механических систем со структурой дерева имеют рекуррентную структуру. Для них в настоящее время разработаны достаточно эффективные методы анализа: метод составных тел для уравнений движения в форме уравнений Лагранжа II рода [5, 15], методы «прогонки» (отдельных тел) для расширенной системы общих уравнений динамики в декартовых координатах и уравнений движения в форме уравнений Лагранжа I рода [1, 5], настроенные на решение механических задач итерационные методы решения систем линейных алгебраических уравнений (СЛАУ) с плотно заполненной или разреженной матрицей системы [2, 3] и др.
Однако структуры большинства моделей технических систем включают замкнутые циклы. Наличие дополнительных связей, порождающих замкнутые циклы, нарушает рекуррентную структуру уравнений движения и приводит к необходимости организации дополнительных матричных вычислений для исключения зависимых переменных в процессе интегрирования. Это влечет существенное увеличение времени, затрачиваемого на проведение вычислительных экспериментов.
Описываемая в данной работе новая методика учета дополнительных голономных связей, использующая технику модифицированных функций Лагранжа, позволяет включить в область применения эффективных алгоритмов численного моделирования, разработанных для систем со структурой дерева, системы с замкнутыми кинематическими цепями.
Постановка задачи
Требуется исследовать динамическое поведение механической системы (системы твердых тел), дифференциальные уравнения (ДУ) движения которой в обобщенных или абсолютных координатах имеют следующий вид:
(1)
где 
– вектор обобщенных (абсолютных) координат; 
– вектор скоростей; 
– вектор ускорений; 
– матрица инерции, 
– вектор внешних и внутренних активных сил и сил инерции; n – число степеней свободы или число уравнений (1).
Пусть на систему наложено m дополнительных линейно-независимых голономных склерономных связей (m < n):
g(q) = 0, (2)
где g ∈ Rm – дважды дифференцируемая вектор-функция.
Требуется найти частное решение системы уравнений (1) при заданных дополнительных ограничениях (2) с начальными условиями: q(0) = q0, 
В настоящее время для учета дополнительных связей при компьютерном моделировании динамики механических систем используются следующие методы (подробные обзоры по этому вопросу можно найти в работах [10, 14]): классический метод Лагранжа [15]; метод стабилизации связей Баумгарта [12, 15]; методы замены связей упруго-демпфирующими элементами (УДЭ) [4]; методы учета связей с помощью штрафных функций (ШФ) [11]; численные методы прямого интегрирования дифференциально-алгебраических систем уравнений (ДАУ) [6]; методы разделения координат q(t) на зависимые и независимые и проектирования уравнений на подпространства независимых координат [8] и др.
Среди основных недостатков указанных методов отметим следующие: необходимость решения дополнительных СЛАУ, возможность ухода со связей, появление нежелательных дополнительных высокочастотных малых колебаний в системе, искажение реакций связей и уменьшение скорости численных расчетов в результате дробления шага интегрирования.
Метод модифицированных функций Лагранжа
При выводе уравнений движения из принципа Гамильтона – Остроградского дополнительные связи можно учитывать методом модифицированных функций Лагранжа (МФЛ) в форме Пауэлла [13]. Тогда уравнения движения преобразуются к виду
(3)
где e ∈ Rm – вектор модифицированных множителей Лагранжа; 
– матрица базиса многообразия ортогонального к множеству векторов, удовлетворяющих уравнениям связей (2); α – штрафной коэффициент.
Заметим, что в точке минимума действия по Гамильтону существует связь между обычными λ и модифицированными множителями Лагранжа e вида λ = αe. Это означает, что для определения последних уравнения движения можно замкнуть дважды продифференцированными по времени уравнениями связей (2), как и в классическом варианте. При этом для стабилизации связей можно использовать метод Баумгарта. Таким образом, в окончательном виде система динамических уравнений включает в себя уравнения (3) и замыкающие уравнения:
(4)
где C, D – положительные диагональные матрицы коэффициентов жесткости и демпфирования колебаний невязок связей g(t) вокруг тривиального решения.
С другой стороны, для пересчета компонент вектора e можно использовать различные приближенные итерационные соотношения, например формулу с линейной скоростью сходимости
,
где ku, kw – коэффициенты усиления обратной связи соответственно по нарушению и ускорению нарушения связей.
Если ввести новые переменные u:
то уравнения движения (3) преобразуются к виду
(5)
(6)
где A = αku; B = αkw.
Уравнения в форме (5), (6) можно рассматривать как уравнения системы автоматического регулирования, где управление выбирается в форме ПИД-регулятора [7]. Физически это означает, что к системе добавляются некоторые силовые элементы с управляемыми параметрами, которые за счет расширения фазового пространства позволяют обеспечить заданную точность моделирования.
Выбор параметров ПИД-регулятора
Настройка управляющих слагаемых в уравнениях (5), (6) заключается в подборе трех векторных параметров ku, С, D. Для этого можно использовать какой-либо из стандартных методов теории автоматического регулирования, например метод Зиглера – Никольса [7].
Однако для решения поставленной задачи более пригодным оказывается другой подход. Пусть на механическую систему накладывается только одна дополнительная связь (2). Если из системы уравнений (5), (6) исключить компоненты вектора состояний механической системы, удалить внешнее силовое воздействие, а затем линеаризовать эту систему, то получим линейное однородное ДУ, описывающее собственные колебания модифицированного множителя Лагранжа е:
(7)
Потребуем, чтобы характеристическое уравнение для ДУ (7) имело вид
(8)
где β0, β1 – коэффициенты демпфирования; ν – круговая частота регулятора. Предположим, что ν > νmax, β0 = O(p ν), b1 = O(p n), p = O(1), т.е. время переходного процесса в системе управления должно быть меньше, чем период наивысшей частоты колебаний механической системы. Сопоставляя характеристическое уравнение для ДУ (7) и уравнение (8), получаем следующие параметры ПИД-регулятора:
(9)
Равенства (9) задают коэффициенты усиления обратной связи, демпфирующие и упругие параметры ПИД-регулятора (6).
С практической точки зрения преимущество уравнений (5), (6) заключается в том, что система уравнений (5) имеет ту же структуру, что и исходная система (1) без дополнительных связей. Это означает, что для ее разрешения относительно ускорений можно использовать те же методы, что и для исходной системы уравнений (1).
Уравнения (5), (6) построены таким образом, что при возникновении отклонений Δg(q) изменяется в первую очередь ненапряженная длина дополнительного упруго-демпфирующего элемента 
так, чтобы создаваемая этим элементом дополнительная сила обеспечивала скольжение механической системы по дополнительной связи g(q) ≈ 0.
Вычислительные эксперименты
Эффективность различных подходов к учету дополнительных связей исследовалась на двух модельных механических системах: первая – двойной многозвенный математический маятник (рис. 1), вторая – квадратная решетка, составленная из связанных плоских пластин (рис. 2).
Рис. 1. Многозвенный маятник
Рис. 2. Решетка твердых тел
Сравнивались характеристики следующих методов: a – метод Лагранжа; b – метод Баумгарта; c – метод замены связей УДЭ; d – метод МФЛ и стабилизации связей по Баумгарту (МФЛ + Б), уравнения (3), (4); e – метод МФЛ с управляющими функциями (МФЛ + У), уравнения (5), (6).
Вычисления проводились в среде САВ Mathematica. Интегрирование уравнений движения выполнялось с помощью встроенной процедуры NDSolve с установленными по умолчанию параметрами. Заметим, что во всех случаях из-за наличия множителей Лагранжа эта процедура интегрировала предлагаемые уравнения как дифференциально-алгебраические. Поэтому в расчетах алгоритм интегрирования ДАУ отдельно не выделялся.
Первая система представляет собой два многозвенных математических маятника с четным числом тел (рис. 1), связанных невесомыми нерастяжимыми балками и упруго-демпфирующими связями между массами. В результате получается система с замкнутыми кинематическими цепями.
Уравнения движения для каждого маятника строились в форме уравнений Лагранжа II рода. Дополнительные связи были представлены в следующем виде:
k = 1, 2, ..., n/2,
где l – длина любого звена; h – расстояние между цепочками; xi = xi(q); yi = yi(q) – абсолютные декартовые координаты тел; q = (q1, ..., q2n) – вектор абсолютных углов отклонения звеньев от вертикали (обобщенных координат); 2n – число тел (точечных масс) в системе.
Обобщенные силы от упруго-демпфирующих связей в уравнениях были определены формулами вида
i = 1, 2, ..., 2n,
где
Расчеты проводились для следующих значений параметров модели: n = 6, l = 1,8 м; h = 4 м; mi = 10 кг; c = 50 000 Н/м (коэффициент жесткости), d = 100 Н•с/м (коэффициент демпфирования). Время моделирования составило 20 с.
Конструктивные параметры методов выбирались таким образом, чтобы в графических представлениях результатов расчетов переменных состояния механической системы, полученных разными методами, не наблюдалось значительных визуальных отличий.
На рис. 3 приведены графики колебаний нечетных звеньев маятника.
Следующий рисунок в логарифмическом масштабе (рис. 4) представляет графики невязок связей 
, полученных различными методами.
Графики на рис. 4 показывают, что:
а) метод Лагранжа приводит к постепенному сходу со связей;
б) наименее точным является метод УДЭ. Он хуже других методов справляется с компенсированием отклонений в положении звеньев системы, вызванных действием силы тяжести;
в) лучшие результаты показывают методы МФЛ и Баумгарта.
Рис. 3. Колебания нечетных звеньев маятников
Рис. 4. Графики невязок связей. Методы: a – Лагранжа; b – Баумгарта; c – УДЭ; d – МФЛ + Б; e – МФЛ + У
Таблица 1
|
Метод |
Лагранжа |
Баумгарта |
УДЭ |
МФЛ + Б |
МФЛ + У |
|
Время, с |
2,83 |
3,78 |
13,63 |
2,19 |
5,77 |
В табл. 1 представлены сведения о затратах компьютерного времени на моделирование системы 1 различными методами. Из анализа этой таблицы следует, что наиболее быстрым является метод МФЛ+Б, а самым медленным – метод УДЭ.
В качестве второй системы была выбрана квадратная решетка, состоящая из 40 одинаковых твердых тел в форме тонких квадратных пластин (рис. 2). На этом рисунке для наглядности изображена только часть решетки, содержащая 21 тело. Ширина каждой пластины равна а = 1/7 м, толщина – h = 0,01 м, масса – m = 1,25 кг. Каждое тело системы связано с соседними телами шаровыми шарнирами, расположенными в центрах боковых граней тел. В центре каждой группы из 9 звеньев решетки отсутствует одно тело, что обеспечивает статическую определенность всей системы. Первое тело системы закреплено шаровым шарниром в точке этого тела с координатами (–а/2, –а/20, 0) в абсолютной системе координат. В начальном положении решетка располагается в первой четверти плоскости OXY. Ось OZ направлена вниз вдоль направления действия силы тяжести.
Данная конструкция взята из работы [9]. Ранее она использовалась для сравнения эффективности программных комплексов «Универсальный механизм – UM», «ФРУНД», «EULER», предназначенных для моделирования систем твердых тел.
Для обеспечения жесткости всей конструкции в шаровых шарнирах, связывающих тела системы, заданы линейные упруго-демпфирующие моменты с коэффициентами жесткости c = 100 Н•м/рад и демпфирования d = 0,01 Н•м•с/рад. Время моделирования составило 4 с.
Все шарниры в системе рассматривались как дополнительные связи. При отсутствии дополнительных связей каждое тело системы имело 6 степеней свободы, включая 3 линейные – абсолютные координаты центров тяжести тел ri = (xi, yi, zi) 
и 3 угловые – углы Эйлера qi = (Ψi, θi, φi).
Уравнения движения системы без связей строились в форме уравнений Лагранжа II рода, а уравнения дополнительных связей в шаровых шарнирах для рассматриваемой системы выглядели так:
где Ij – матрицы преобразования систем координат, связанных с телами в абсолютную систему координат. Индексы j и k перебираются по всем номерам соседних тел в решетке. Обобщенные силы от упруго-демпфирующих связей имеют вид
Рис. 5. Линейные колебания первых пяти звеньев системы
Рис. 6. Графики невязок связей. Методы: a – Лагранжа; b – Баумгарта; c – УДЭ; d – МФЛ + Б; e – МФЛ + У
Таблица 2
|
Метод |
Лагранжа |
Баумгарта |
УДЭ |
МФЛ + Б |
МФЛ + У |
|
Время, с |
17,02 |
19,77 |
72,83 |
14,17 |
56,39 |
На рис. 5 приведены графики абсолютных линейных координат центров тяжести первых пяти тел системы.
На следующем рисунке (рис. 6) в логарифмическом масштабе представлены графики невязок дополнительных связей 
, описывающих первые четыре шарнира.
Графики на рис. 6 показывают, что:
а) метод Лагранжа приводит к постепенному сходу со связей;
б) наименее точным является метод замены связей упруго-демпфирующими элементами;
в) лучшие результаты показывают методы Баумгарта и МФЛ.
В табл. 2 представлены сведения о затратах компьютерного времени на моделирование системы 2 различными методами. Из анализа этой таблицы следует, что наиболее быстрым является метод МФЛ+Б, а самым медленным – метод УДЭ.
В результате проведенных численных экспериментов было установлено, что предлагаемый новый метод учета дополнительных связей в механических системах, основанный на модифицированных функциях Лагранжа и алгоритмах адаптивного управления, имеет определенные преимущества перед существующими схемами учета дополнительных связей.