Приведем некоторые соотношения, отражающие свойства матриц поворота, которые будут использоваться далее:

преобразование матричных представлений векторов при переходе от СКj к где x – любой вектор в трехмерном пространстве, верхний индекс в скобках здесь и далее означает систему координат, в которой представлен вектор;

матрица поворота при переходе от СКi к СКj равна транспонированной матрице обратного перехода, т.е. A ij = A T ;

при переходе от СКi к СКk в цепочке тел i,j,k матрица поворота может быть вычислена как результат матричного произведения A ik = A ij A jk.

Таким образом, положение i-го тела относительно СК0 задается радиус- вектором начала отсчета СКi r0i (рис. 1.1) и матрицей поворота A0i. Аналогично определяется положение i-го тела относительно любой СКj.

Для описания кинематики механической системы вводится понятие шарнир, который в общем случае трактуется как способ описания относительного движения некоторого тела j системы относительно тела i. Далее будем говорить, что шарнир соединяет или связывает пару тел i и j, одним из которых может быть базовое тело 0. Каждому шарниру соответствует набор локальных шарнирных координат qij. Размер вектора qij зависит от типа шарнира и равен минимальному числу координат, необходимых для описания конкретного типа связи. Модели механических систем могут быть построены с использованием шарниров двух типов: нормального и контактного.

Шарнир нормального типа, соединяющий тела i, j, может быть представлен как последовательность элементарных преобразований (сдвигов или поворотов), посредством которых осуществляется переход от СКi к СКj при любом относительном положении пары тел, допускаемом данным шарниром.

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

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

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

Перенумеруем тела и шарниры механической системы в произвольном порядке. Объединим все локальные шарнирные координаты в один векторстолбец q = [qij ], расположив их в порядке нумерации шарниров. Шарнирные координаты q механической системы со структурой дерева являются обобщенными координатами, то есть они независимы и однозначно определяют кинематические характеристики системы. В случае систем с замкнутыми кинематическими цепями набор координат q избыточен, что приводит к необходимости добавления уравнений связей.

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

1.1.2.2. Описание относительной кинематики пары тел, связанных шарниром Соотношения локальной кинематики шарнира описывают движения тела j, относительно тела i.

Положение СКj относительно СКi задается радиус-вектором начала отсчета rij и матрицей поворота Aij:

Конкретный вид выражений (1.1) определяется типом шарнира.

Соотношение для относительной линейной скорости получим дифференцированием первого из соотношений (1.1) по времени:

где vij – скорость начала отсчета СКj относительно СКi, суммирование распространяется по всем координатам в шарнире.

Обозначив diji,)m = riji ) qij,m, v (i ) = riji ) t и расположив dij,m по столбцам матрицы Dij, выражению (1.2) можно придать следующий вид:

Заметим, что Dij является матрицей Якоби.

Выражение для относительной угловой скорости получим, рассмотрев соотношения для скорости произвольной точки K тела j, положение которой задается вектором К, постоянным в системе координат j. Тогда в СКi положение точки K задается следующим соотношением:

Дифференцируя выражение (1.4) по времени, получим скорость точки K в СКi:

Кососимметрическая матрица A ij A ji, появившаяся в соотношении (1.5) в результате этой операции, определяет угловую скорость тела j относительно тела i.

Введем обозначение iij = A ij A ji. Вектор угловой скорости ij строится по кососимметрической матрице ij в соответствие с правилом:

где нижние индексы опущены для удобства записи.

Заметим, что iji ) (i ) определяет векторное произведение ij k. Это правило распространяется на любые векторы в трехмерном пространстве.

Запишем выражение угловой скорости в виде:

ния biji,)m, (i ) соответственно и переходя к векторам, соотношение для относиij тельной угловой скорости в СКi перепишем в следующем виде:

где по столбцам матрицы Bij расположены векторы-столбцы b ij,m.

Относительное ускорение полюса получим дифференцированием соотношения (1.3) по времени:

Введя обозначение a(i ) = Diji )q ij + v(i ), перепишем выражение для aij в виде:

Аналогично получим соотношение для относительного углового ускорения:

Заметим, что если выражения (1.1) не содержат время t в явном виде (склерономная связь), векторы-столбцы v, отсутствуют.

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

1.1.2.3. Уравнения кинематики системы со структурой дерева Выражения, определяющие кинематику каждого i-го тела механической системы относительно СК0, можно представить в следующем обобщенном виде:

Они строятся рекуррентно на основе анализа графа системы с использованием соотношений (1.1,1.3, 1.6-1.8). Здесь ri, vi, ai – радиус-вектор, скорость и ускорение начала отсчета СКi соответственно; i, i – угловые скорость и ускорение тела.

В графе системы, имеющей структуру дерева, имеется единственный простой путь от корневой вершины, соответствующей базовому телу 0, к любой iой вершине, соответствующей i-ому телу. Как было отмечено выше, в этом случае q является набором обобщенных координат. Рассмотрим кинематику тела с номером k+1. Предположим, что тела в пути от корня графа к телу k+1 пронумерованы последовательно 0..k. Тогда справедливы соотношения:

где rk(,kk)+1 (q k, k +1, t ), A k,k +1 (q k,k +1, t ) известны из анализа локальной кинематики шарнира.

Если выполняется символьный синтез соотношений кинематики, выражения для линейных и угловых скоростей и ускорений получаются дифференцированием соотношений (1.9) по времени. При численном выводе кинематические характеристики выражаются в явной форме в зависимости от q. Вид этих выражений приводится в [28].

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

где вектор-столбец qb объединяет локальные шарнирные координаты, соответствующие шарнирам дерева системы, q a объединяет координаты, соответствующие разрезанным шарнирам. Координаты, составляющие вектор-столбец qb, называются основными координатами системы, а составляющие векторстолбец q a – вспомогательными. Далее вывод уравнений кинематики и динамики для систем с замкнутыми цепями строится в терминах qb, поскольку положение всех тел системы однозначно определяется значением основных координат. При этом алгоритм синтеза кинематических соотношений для механических систем со структурой дерева и систем с замкнутыми цепями после разрезания соответствующих шарниров одинаков. Однако координаты qb являются избыточными (число степеней свободы меньше размерности векторастолбца). Они удовлетворяют нелинейным уравнениям связей:

где G = g qT – матрица Якоби. Если в системе отсутствуют реономные связи, вектор-столбец g отсутствует.

Вывод уравнений связей в форме (1.10) часто бывает затруднительным.

Удобнее вывести уравнения, в которые входят как основные, так и вспомогательные координаты:

значения индексов в обозначении векторов-столбцов координат соответствуют первым буквам англоязычных терминов basic (основной) и auxiliary (вспомогательный) Вычисление вспомогательных координат во многих случаях полезно, а иногда просто необходимо, например, когда от их значения зависят активные силы, действующие на систему.

Заметим, что уравнения (1.10,1.11) представляют собой голономные связи механической системы. Для неголономных связей алгебраическая часть уравнений отсутствует, а дифференциальным уравнениям связей в большинстве случаев можно придать форму (1.11.2, 1.11.3). Неголономные связи соответствуют, как правило, качению тел без проскальзывания и редко используются для моделирования технических задач. В настоящем обзоре они не рассматриваются.

Исключив из уравнений (1.11) вспомогательные координаты, получим уравнения связей в стандартной форме (1.10). Для этого разделим уравнение (1.11.1) на две части:

так, чтобы матрица Якоби h 2 q a была квадратной и невырожденной. Тогда, используя второе из полученных уравнений, вспомогательные координаты могут быть выражены через основные:

Подставив это соотношение в первое уравнение, получим:

Аналогично исключаются производные по времени от вспомогательных координат из уравнений (1.11.2, 1.11.3).

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

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

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

где нижний индекс означает номер тела, mi - масса тела, Ji – матрица тензора инерции, f ie, l ie, f ir, l ir – главные векторы и главные моменты относительно центра масс активных сил и сил реакций связей соответственно. Уравнения (1.12.1, 1.12.2) представлены в разных системах координат с целью уменьшения числа операций при выводе.

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

и означает равенство нулю суммы работ сил реакций на возможных перемещениях системы. Здесь ri( 0), i(i ) – векторы перемещения центра масс и поворота тела i при возможном перемещении системы. Выразим их через вариации координаты системы:

Тогда, подставляя в уравнение (1.13) значения f ir, l ir, выраженные из уравнений (1.12.1, 1.12.2), с учетом соотношений (1.14) получим уравнения движения механической системы в следующем виде:

Дальнейшие преобразования уравнения (1.15) зависят от структуры механической системы.

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

Поскольку вариации обобщенных координат независимы, справедливо следующее уравнение:

Используя соотношения (1.7, 1.8), введем следующие обозначения:

где M – матрица масс системы, k, f – векторы-столбцы обобщенных сил инерции и обобщенных активных сил. Тогда уравнению (1.16) можно придать более компактный вид:

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

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

Во-вторых, процедура вычисления M, k, f включает независимые расчеты соответствующих величин M i, k i, fi для каждого тела системы и последующее суммирование этих величин.

Уравнения движения систем с замкнутыми кинематическими цепями.

В соответствие с подходом, описанным при анализе кинематики, перейдем к системе со структурой дерева, разрезав соответствующие шарниры. Тогда обобщенные силы от сил реакций в разрезанных шарнирах войдут в столбец обобщенных сил f. Уравнения движения систем с замкнутыми кинематическими цепями строятся в терминах основных шарнирных координат qb, которые являются избыточными и вместе с производными по времени первого и второго порядка удовлетворяют уравнениям связей (1.10). Тогда для вариации координат справедливо соотношение С учетом этого можно доказать, что уравнение (1.15) удовлетворяется, если выполняется следующее соотношение:

где – вектор-столбец множителей Лагранжа, размерность которого соответствует числу уравнений связей. Множители Лагранжа определяют силы реакций в разрезанных шарнирах. Явный вид этой зависимости определяется типом шарнира.

Уравнения (1.19) совместно с (1.10) образуют замкнутую систему уравнений движения механической системы. С учетом обозначений (1.17) ей можно придать следующий вид:

Неизвестными в этих уравнениях являются координаты q b и множители Лагранжа.

Как упоминалось выше, уравнения (1.20) могут быть синтезированы автоматически с применением формализмов в численной или символьной форме. Численная форма уравнений предполагает их вывод на каждом шаге интегрирования. Уравнения в символьной форме синтезируются один раз для всего цикла исследований. Синтез повторяется только после изменения модели.

На рис. 1.2 представлены схемы исследований с применением упомянутых формализмов.

Рис.1.2. Схема применения численного и символьного формализмов для исследования модели посредством интегрирования уравнений движения В заключение отметим весьма высокую эффективность моделирования динамики механических систем на основе их представления системой абсолютно твердых тел. Метод хорошо зарекомендовал себя в различных областях науки и техники [12, 22, 78].

§ 1.2. Обзор подходов к формированию уравнений движения упругих тел Решение многих задач моделирования динамики механических систем с необходимой точностью возможно только при учете упругих свойств отдельных частей конструкций. Например, значения низших собственных частот кузова железнодорожного экипажа вполне сопоставимы с частотами колебаний на подвеске. Поэтому учет форм колебаний кузова, соответствующих низшим частотам, может существенно уточнить результаты анализа вибраций при движении. С другой стороны, значения низших частот железнодорожного колеса составляют 180-200 Гц, то есть весьма далеки от низших частот экипажа. Применение единого подхода для всей конструкции в такой ситуации нерационально, желательно совмещение различных подходов при построении модели.

Пример подобной задачи из аэрокосмической области приводится в работе Швертассека [82], где с учетом упругости моделируется вращение антенны, присоединенной к станции. В работе Амброзио и Перейры [40] комбинированный подход применяется для исследования динамики автомобиля.

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

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

Подход, называемый методом конечных сегментов, или методом твердотельных элементов, представлен в п.1.2.1.

Следующим шагом по пути развития методов моделирования динамики упругих тел является попытка совмещения подходов СТТ и МКЭ в рамках линейной теории динамики упругих тел (ЛТДУТ). МКЭ появился в начале 1950-х годов и первоначально разрабатывался как метод анализа статики конструкций.

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

где М,D,C постоянные, симметричные матрицы масс, внутреннего демпфирования и жесткости соответственно, f – вектор приложенных сил. Размер системы уравнений N соответствует числу узловых степеней свободы.

В настоящее время МКЭ оформился как классическая дисциплина и мощный, широко применяемый инструмент исследования динамики, статики и напряженно-деформированного состояния конструкций. Множество разработанных типов конечных элементов позволяет моделировать тела практически любой конфигурации. Описанию основ и развитию метода посвящены многочисленные публикации, среди которых можно отметить работы Зенкевича [16], Галагера [6], Постнова [33, 35]. Метод реализован в виде программного обеспечения. Самыми известными пакетами являются ANSYS, NASTRAN, ABAQUS.

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

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

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

Далее рассмотрим упомянутые методы более подробно. Кроме того, остановимся на методе последовательных приближений, который широко применялся до появления методов ПСК и АК, и методе вектора большого поворота, не получившем широкого распространения.

1.2.1. Метод твердотельных элементов В соответствии с методом твердотельных элементов упругие компоненты конструкции представляются набором твердых тел, соединенных упругодиссипативными элементами. Основными проблемами являются выбор числа твердых тел, их расположения, а также подбор характеристик упругодиссипативных элементов. Значения параметров могут выбираться, например, с использованием МКЭ. Другой проблемой является интегрирование жестких уравнений движения систем, включающих сотни тел.

В этой области можно отметить работы Крушевского и др. 1975 г. [63], Хьюстона и др. 1981 г. [62], Рау и Шилена 1986 г. [76], Леонтьева 1992 г. [18], Паскаль и Гагариной [70], Погорелова и др. [73, 74, 30, 9].

В большинстве работ рассматриваются балочные схемы, в том числе пространственные. Например, в работе [9] балочные элементы успешно применяются для определения положений равновесия и расчета собственных форм кабеля жесткого диска компьютера. В работе [50] подобная схема используется для исследования динамики ворсинок вращающейся щетки.

В работах [74, 30] предприняты попытки моделирования пластин на основе твердотельных элементов. В работе [74] показана сходимость для задач статики, однако в задаче расчета собственных частот и форм колебаний пластины наблюдается плохая сходимость. Эта закономерность подтверждается работой [30], где моделируется конвейер с подвесной лентой. Модель содержит 200 абсолютно твердых тел, размер системы уравнений движения равен 1000.

1.2.2. Линейная теория динамики упругих тел Предполагается, что упругие перемещения за счет деформаций не имеют значительного влияния на движение тела как абсолютно твердого [85]. Моделирование проводится в два этапа. На первом этапе вся конструкция представляется твердотельной моделью, с использованием которой вычисляются значения сил инерции и сил реакций связей. Эти значения служат исходными данными для второго этапа – решения уравнения (1.21) для упругих компонентов, где f представляется суммой сил, рассчитанных на первом этапе f i + f r + f e. Использованы следующие обозначения: f i – инерционные силы, в том числе силы Кориолиса и центробежные силы, f r – силы реакций, f e – активные силы. Движение деформируемых тел представляется суммой (суперпозицией) больших твердотельных движений и малых упругих перемещений.

Результаты исследований с применением описанного метода приводятся, например, в работе [17].

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

1.2.3. Метод последовательных приближений.

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

Варианты численных процедур, реализующих описанный подход, представлены в работах [55, 60].

Классическим примером использования метода последовательных приближений является часто цитируемая задача деформирования консольного стержня под воздействием момента, приложенного на свободном конце [44, 85].

Этот метод плохо моделирует твердотельную динамику при вращении тел. К тому же он не приводит к нулевым напряжениям при произвольном движении тел как абсолютно твердых. Для устранения этих недостатков используется ряд приемов, которые, как правило, сводятся к введению дополнительных систем координат [85].

1.2.4. Векторы больших поворотов Предлагается подход, при котором вращательные координаты в узлах позволяют описывать большие повороты сечений конечного элемента. Интерполяционные функции форм строятся обычным образом. То есть используется процедура классического МКЭ без линеаризации по углам поворота [89, 90].

Применение подхода к некоторым типам конечных элементов, например, к балкам Эйлера-Бернули, приводят к избыточности координат и вырождениям [83]. Вследствие этого, а также по причине появления более прогрессивного метода формулировки МКЭ в абсолютных координатах представляемый подход не получил широкого распространения.

1.2.5. Метод присоединенной системы координат Реализация метода является естественной попыткой распространить кинематическое соотношение (1.4) на упругие тела. Пространственное движение упругого тела представляется посредством связанной с ним системы координат 1, относительно которой строятся соотношения для упругих перемещений. Положение произвольной k-ой точки деформируемого тела в глобальной СК определяется следующим соотношением (рис 1.3):

Рис 1.3. Метод присоединенной системы координат Вектор d представляет перемещения за счет упругих деформаций тела. Соотношение (1.22), вообще говоря, не налагает никаких ограничений на их величину. Деформации и соответствующие перемещения могут быть описаны любым корректным способом. Однако описываемый подход получил широкое распространение при условии малости упругих перемещений, то есть в СК конструкция может быть представлена линейным МКЭ.

Используя любой известный способ построения уравнений движения, например, уравнения Лагранжа II рода или общее уравнение динамики, уравнения движения упругого тела можно записать в виде:

где M, С, D – матрицы масс, упругости и внутреннего демпфирования упругого тела, q – обобщенные координаты, k – вектор-столбец обобщенных сил инерции, fa – вектор-столбец обобщенных сил от приложенных нагрузок, fc – вектор-столбец обобщенных сил реакций. Число обобщенных координат равно 6+N; 6 координат задают положение СК1 в пространстве, N узловых координат используются для описания малых упругих перемещений в СК1. Уравнения, построенные методом ПСК, содержат нелинейные матрицу масс и вектор сил инерции как результат инерционной связи между большим твердотельным движением и малыми упругими перемещениями. Матрица жесткости, с другой стороны, не отличается от матрицы, построенной линейным МКЭ, благодаря факту, что упругие координаты определяются в СК тела.

Число узловых координат конечно-элементной модели для многих конструкций может равняться десяткам и сотнями тысяч. Этот факт, а также наличие высокочастотных составляющих в решении делают непосредственный анализ системы уравнений (1.23) путем интегрирования неэффективным. Ее размер можно уменьшить на несколько порядков, используя модальный анализ [86].

Общий вид уравнений, записанных в терминах модальных координат, совпадает с (1.23), однако вместо N узловых координат упругие свойства представляют H новых координат, соответствующих допустимым формам упругого тела,