CC BY
219
_________УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XXXI 2000
№ 3—4
УДК 519.63:512.643 539.3
МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ В ЗАДАЧАХ ДИНАМИКИ СВОБОДНЫХ КОНСТРУКЦИЙ
Л. П. Лущин, А. В. Шаранюк
При анализе динамического поведения сложных составных конструкций используется метод конечных элементов (МКЭ). Среди прочих преимуществ этот метод позволяет с единых позиций рассматривать напряженно-деформи-рованное состояние конструкции, состоящей из разнородных фрагментов, деформации которых описываются разными системами дифференциальных уравнений. Особенно часто такая ситуация возникает при расчете тонкостенных конструкций авиационной и космической техники.
Динамическое поведение подобных конструкций связано с нестационарными воздействиями двигателей, органов управления, взаимодействием с окружающей средой. Подобные воздействия приводят к изменению характера движения конструкции. Возникают динамические нестационарные нагрузки, приводящие к появлению напряжений и деформаций конструкционных элементов.
Анализ напряженно-деформированного состояния движущихся и маневрирующих конструкций часто бывает удобно проводить в подвижной системе координат, связанной с конструкцией. В отличие от системы координат, связанной с Землей, такая система координат не является инерциальной, вследствие чего в динамических уравнениях, описывающих нестационарные деформации конструкции, появляются дополнительные составляющие.
1. Движение свободной конструкции. Будем рассматривать две системы координат. Инерциальную систему , в качестве которой можно выбрать, например, систему
координат, связанную с Землей, и подвижную систему координат 0x^x2Х3, связанную с движущейся деформируемой конструкцией. Единичные векторы осей этой системы координат обозначим через в\, в2, Є3. Направляющие косинусы, задающие ориентацию осей подвижной системы координат относительно неподвижной, могут быть заданы в виде матрицы [С], состоящей из компонентов единичных векторов в инерциальной системе координат
[ С] = [ е1е2е3 ] •
С помощью матрицы [С] компоненты любого вектора а, заданного
в инерциальной системе координат, могут быть преобразованы в компоненты вектора а' в подвижной системе координат и наоборот
а' = [ С ]Та , а = [С ]а'. (1.1)
В дальнейшем векторы и матрицы, компоненты которых заданы в подвижной системе координат, будем помечать штрихом.
Каждому вектору а можно поставить в соответствие кососимметрическую матрицу [а], имеющую следующую структуру:
" 0 - а3 а2
[ а ] = а3 0 - а1
а2 а1 0
Непосредственной проверкой можно убедиться, что
Х1
axb = [ a]b,
т. е. векторное произведение двух векторов можно заменить на произведение матрицы на вектор.
В дальнейшем приняты следующие обозначения. Матрицы обозначаются квадратными скобками. Векторы — жирными буквами. Жирные символы, заключенные в квадратные скобки, обозначают кососимметрические матрицы, образованные компонентами вектора, заключенного в квадратные скобки.
Предположим, что связанная система координат ОХ1Х2Х3 вращается с угловой скоростью ш. В силу того, что Є = шхвї (і = 1, 2, 3) с учетом структуры матрицы [С], получим
(точкой обозначается дифференцирование по времени).
Вектором р задается положение произвольной точки недеформированной конструкции относительно начала связанной с конструкцией системы координат. Пространственное положение произвольной точки недеформированной конструкции в неподвижной системе координат задается с помощью вектора Я, состоящего из суммы векторов До, характеризующего положение начала
связанной системы координат, и вектора р. В процессе движения упругая конструкция претерпевает деформации, которые характеризуются вектором упругих перемещений и (см. рисунок),
Поступательное и вращательное движения упругой конструкции будем исследовать в инерциальной системе координат 2^3, а упругое поведение — в подвижной системе
0Х[Х2Х3. Положение упругой конструкции относительно инерциальной системы координат характеризуется положением начала связанной системы координат, которое определяется вектором Rq , а ориентация задается матрицей косинусов [С].
Определение угловой ориентации связанной системы координат. Матрица косинусов [С] содержит девять параметров, подчиненных шести тождествам, т. е. имеется только три независимых параметра, характеризующих угловую ориентацию связанной с конструкцией системы координат. Существует множество способов выбора такой тройки: углы Эйлера, углы Крылова, направляющие косинусы, параметры Родрига — Гамильтона и т. д. Использование углов Эйлера или Крылова приводит к необходимости интегрирования системы трех дифференциальных уравнений. Общим недостатком введения трех углов в качестве независимых параметров ориентации связанной системы координат является то, что кинематические уравнения, разрешенные относительно производных углов, содержат нули в знаменателях. Направляющие косинусы Cij (i, j = 1, 2, 3), как видно из соотношения (1.2),
определяются интегрированием кинематических уравнений Пуассона.
Использование кватернионов для задания угловой ориентации связанной системы координат позволяет уменьшить количество кинематических уравнений до четырех, которые не вырождаются при любой ориентации связанной системы координат.
[ C ]=Ы[ c ]
(1.2)
R(t) = R o(t) +p (t)+u (t).
(1.3)
Применение кватернионов основано на теореме Эйлера о том, что произвольное вращение твердого тела вокруг неподвижной точки можно осуществить одним поворотом вокруг соответствующим образом выбранной оси, проходящей через эту точку. Матрицу направляющих косинусов можно выразить через компоненты кватерниона [1]
+ А2 — А^ — А2 2 (А]А2 — АоА3 ) 2 (А]А3 — АоА2 )
2 (А]А2 + АоА3 ) Ар — А2 + А"2 — А2 2 (А2А3 — АоА]) .
2 (А] А3 — Ао А2 ) 2 (А2А3 + Ао А]) А^ — А2 — А"2 + А2
Кинематические соотношения для компонентов кватерниона имеют вид
2Ао — —Ш]А] — Ш2А2 — Ш3А3 ,
2А] — Ш]Ао + Ш3А2 —Ш2А3 ,
2А2 — Ш2Ао + Ш]Аз — шзА] ,
2Аз — шзАо + Ш2А] — Ш]А2 .
Уравнения движения свободной упругой конструкции. Упругие смещения элементов конструкции удобно исследовать в связанной системе координат. Используя соотношение (1.3), получим
щ) — яо(0+[ С](р'+«'(?)). (1.6)
Таким образом, для описания положения точек деформируемой конструкции требуется задание векторов Ло, р' , и и компонентов нормированного кватерниона, определяющего матрицу косинусов [С].
В подвижной системе координат в процессе движения положение любой точки деформируемой конструкции характеризуется вектором р' + и', что позволяет вычислить ее
скорость и' и ускорение и'. Для определения скорости произвольной точки деформируемого тела продифференцируем по времени соотношение (1.6)
Л — Ло + [С](р' + и') + [С] и. (1.7)
Предположим, что связанная с конструкцией система координат вращается с угловой скоростью ю. Тогда, с учетом соотношения (1.2), получим
Л — Ло + [ю][ С] (р' + и' ) + [С] и. (1.8)
Для определения ускорений произвольной точки упругого тела продифференцируем по
времени соотношение (1.8)
Л — Ло + [Ш ][С](р ' + и') + [ш ][с](р ' + и') + [ш ][С]и' + [с] и' + [С] а' (1.9)
и с учетом (1.2) получим
Л — Ло + [С]и ' + [Ш ][С](р ' + и' ) + [ш ] [ш ][С](р ' + и') + 2[ш ][С]и'.
Воспользуемся свойствами кососимметрических матриц [а]Ь = — [Ь]а и перепишем приведенное соотношение следующим образом:
Л — Ло + и + [р + и]Ш + [ш][ш ](р + и) + 2 [ш ] и (1.10)
для компонентов векторов в неподвижной инерциальной системе коор -
динат.
Для построения системы уравнений, описывающих движение свободной упругой конструкции, воспользуемся вариационным принципом возможных перемещений. В соответствии
[ С ] —
с ним работа внешних сосредоточенных, объемных и поверхностных сил на возможных перемещениях 5Л равна работе 5 и напряжений а на возможных деформациях 5е упругого тела
5А — 5и. (1.11)
Возможное перемещение произвольной точки конструкции получим, проварьировав соотношение (1.6)
5Л — 5Ло + [5С](р' + и') + [С]5и'. (1.12)
Умножив правую и левую части полученного равенства на [С]Т , получим выражение для 5Л'
5Л' — 5Ло + [С]Т [5С] (р' + и') + 5и '. (1.13)
Выделим три независимых параметра у; (/ = 1, 2, 3), определяющих матрицу [С], в качестве
которых могут быть три угла Эйлера, три угла Крылова, три компонента кватерниона и т. д. Тогда
[С]Т [5С] = ]Г [С]Т ® 5у; — [50 ' ], (1.14)
дУ/
г^пт д[с ] „
так как оператор [С ] ------5у; — кососимметрическая матрица, в соответствие которой можно
дУ;
поставить вектор 5© '. Свойство кососимметричности матрицы [С]Т д[С] 5у; проверяется
дУ;
следующим образом. Вычислим производную по у; от произведения [С ]Т [С ]
З[С]т[С] д[С]Т [с]+[с]Т д\С] =0
дУ; дУ; дУ; ’
которая равна нулю в силу ортогональности матрицы косинусов. Тогда из последнего равенства следует
®1СЬ[С ]—
дУ;
[С ]
Т д[С]
дУ;
= — [С]
Т д[С]
дУ;
что является признаком кососимметричности матрицы. В силу того, что линейная комбинация кососимметрических матриц есть также кососимметрическая матрица, соотношением (1.14) можно определить вектор малого поворота 5©', и тогда соотношение (1.13) примет вид
5Л' — 5Л + [ 50' ] (р' + и') + 5и'
или в неподвижной системе координат
5Л — 5Ло + [5© ](р + и) + 5и . (1.15)
Воспользовавшись свойством кососимметрических матриц
[5©](р + и) = — [р + и]5©,
вычислим работу, совершаемую объемными Ро и поверхностными Рр силами для тела, занимающего область О с границей Г, на возможных перемещениях
( \
5А — 5ЛТ
+
+ 5©Т
|[р + и] Pоd О + |[р + и] Р^ Г
Г
/
+
+ ^5иТР<о^О+ |"5иТР^Г .
(1.16)
О Г
Запишем выражение для работы, совершаемой внутренними силами демпфирования Род и
Рга, на возможных перемещениях
5Л2 = -|5иТPОddО-|5иТpгddГ •
(1.17)
О
Г
Согласно принципу Даламбера, в систему сил, действующих на тело во все моменты времени, необходимо добавить силы инерции, т. е. объемные силы
Рт = -Р Д >
и вычислить их вклад в работу на возможных перемещениях (р — плотность материала конструкции). Для этого можно воспользоваться соотношением (1.16), в котором Ро заменяется
на Р и исключается Р
Г
5А3 = -5ДТ | ^т - 50Т |[р+ и] Шт -15иТ ^т,
О
О
о
где dm — р(X], Х2, Х3)dО — элемент массы конструкции.
Для вычисления работы массовых инерционных сил на возможных перемещениях 5А3 воспользуемся выражением (1.10), что позволяет записать последнее соотношение следующим образом:
тДо + | Мт - |[р + и^тю + [ю][ш] |(р + и) dm + 2[ш] | Udm>
5Л3 =-5^Т \тКо + I Шт- | |р + и^тю + |ш||©П (Р + и)dm + 2
О О О О
- 5©Т \ | [р + и] dmRо + |[р + и] Мт - |[р + и] [р + и] dmШ +
О
Ю
О
| [р + и][ю][ш](р + и)dm + 21[р + и][ш] Ыт;
+ ^р +
О
О
15иТdmRо + 15иТийт -15иТ [р + и]dmш +
О
О
О
+ 15иТ [ш] [ш] (р + и) dm + 215иТ [ю] Шт I.
(1.18)
О О
Таким образом, вариационный принцип возможных перемещений приводит к соотношению
5Л = 5^1 + 5^2 + 5А3 = Г 5єТad О
(1.19)
О
— <
где правая часть равенства — работа внутренних напряжений конструкции на возможных деформациях.
Соотношение принципа возможных перемещений (1.19) содержит возможное смещение начала связанной системы координат 5Лд, возможный поворот связанной системы координат 50,
возможные деформации конструкции 5в и возможные упругие перемещения конструкции 5и. Возможные деформации конструкции при использовании линейной теории упругости в силу соотношений Коши могут быть выражены через возможные перемещения.
Соотношение (1.19) должно выполняться для любых кинематических допустимых смещений начала связанной системы координат 5Лд, для любых поворотов 50 и любых возможных упругих перемещений. Для выполнения этого необходимо приравнять нулю соответствующие множители при 5Лд, 50 и 5и.
Подставляя в (1.19) соотношения (1.16)—(1.18) и приравнивая нулю множитель при 5Лд, получим
| Р<^й О +1 Руй Г - тЩ - | Ыт + | [р + и] йтса
О Г О О
+ и) йт - 2 [о] / Ийт = 0.
О
(1.20)
О
Аналогично, в соотношении (1.19) приравняем нулю множитель при 50 и получим | [р + и] Р^ О +1 [р + и] Р^ Г - | [р + и] йтк§ - | [р + и] ийт +
О
О
О
+
| [р + и] [р + и] йто - |[р + и] [о][о] (р + и) йт -
О О
- 2 |[р + и] [о] Шт = 0.
(1.21)
О
С учетом уравнений (1.20), (1.21) соотношение принципа возможных перемещений (1.19) принимает вид
15єТа<і О = 15иТ Ро^О. +| 5иТ і~гй Г.
(1.22)
О
О
где приняты обозначения
РО = РО - РШ + Рт >
РГ - РГ - РГс1 •
Для описания напряженно-деформированного состояния конструкции воспользуемся уравнениями линейной теории упругости
°у, ] + ЇІ =
= -( иІ,} + иЦ ), = ЕІ]кігкі >
(1.23)
где агу — тензор напряжений, — тензор деформаций, Ецкі — тензор упругих констант материала, ^ — объемные силы. Здесь предполагается суммирование по повторяющемуся
Г
Г
индексу. Запятая перед индексом обозначает дифференцирование по соответствующей переменной. Система (1.23), дополненная уравнениями совместности деформаций, представляет собой систему линейных дифференциальных уравнений относительно неизвестных функций иI, 8у и Оу (/,] = 1, 2, 3). Напряжения во внутренних точках тела должны непрерывно переходить в напряжения на поверхности, т. е. должны выполняться краевые условия
где пI — компоненты единичного вектора нормали к поверхности, Р™ — нормальные поверхностные напряжения. Проварьировав второе уравнение системы (1.23), получим
8в у =
- [(5иг-) у + (5иу). ],
2
что позволяет ввести линейный дифференциальный оператор В
5є = В5и.
Подставив (1.24) в (1.22) и применив теорему Гаусса, получим
15иТ (Бс + РО) й О-| 5иТ (БГс + РГ) й Г = 0.
О Г
В силу произвольности поля перемещений 5и из (1.25) получим
Бст + Ро = 0,
Бг^ - Рг = 0 .
(1.24)
(1.25)
(1.26)
Эти уравнения устанавливают равновесие между внутренними напряжениями в упругом теле, внутренними объемными силами и внешними поверхностными силами. Структура линейных
операторов В и Вг получается из системы уравнений (1.23). С учетом того, что /1 — Ро/ , подставляя третье уравнение (1.23) (закон Гука) в первое и исключая вгу , получим
- Еи
\}кі (ик + иі к)+роі = 0
с краевыми условиями
- Еукі( ик ,і + иі ,к ) пу = рГі
или в операторном виде
Аи = Ро ,
Аги = Р~г .
2
Таким образом, последнее уравнение, полученное из соотношения принципа возможных перемещений (1.19), имеет вид
Р'ь - Р'ш -р(*6 + [®'] Р ' + [®'][®'] Р ')- Р ([» '] и' + [®'][®']« ')-
- 2р[ю' ] и-р и' = Ли'. (1.27)
Система уравнений (1.20), (1.21) и (1.27), описывающая динамику свободной упругой конструкции, является «гибридной» системой нелинейных уравнений. Решение подобных систем уравнений вызывает определенные трудности, и поэтому возникает необходимость в разработке эффективных численных методов, учитывающих специфику анализируемых механических объектов и процессов.
Метод конечных элементов. При проведении расчетов динамического поведения больших составных конструкций наиболее эффективным вычислительным методом является метод конечных элементов.
Следуя процедуре метода конечных элементов [2], [3], разобьем конструкцию на
непересекающиеся области — элементы, связанные в узлах. Поле перемещений аппроксимируется гладкими в пределах элемента функциями, называемыми функциями формы. В узловых точках дискретной модели записываются условия непрерывности, причем, если используются согласованные элементы, то непрерывность функций, аппроксимирующих поле перемещений, выполняется и вдоль границ. Таким образом, поле перемещений в пределах конструкции аппроксимируется непрерывной кусочно-гладкой функцией. Пусть конструкция разбита на п конечных элементов. Поле перемещений в пределах каждого элемента приближенно может быть представлено следующим образом (верхним индексом «е» будем обозначать принадлежность конкретному элементу)
ие = [С] [у'е ]q'е , (1.28)
где [у ' е ] — матрица функций формы элемента, q'е — вектор обобщенных перемещений в узлах, высота которого равна количеству степеней свободы конечного элемента. По известным перемещениям узлов элемента можно определить деформации, выражения для которых в матричной форме имеют вид
se = [С] [в 'е е , (1.29)
где [в ' e ] — матрица производных функций формы.
Возможные перемещения и возможные деформации в пределах элемента могут быть представлены выражениями
5ие = [С] [уe]дц'e
>
5ве = [С] [в 'е]^'е ,,
(1.30)
т. е. выражены через возможные узловые смещения.
Вектор напряжений, возникающих в элементе при деформировании конструкции, согласно закону Гука связан с деформациями с помощью матрицы упругих постоянных [ж' е]
О 'е = [Же] 8 'е .
Следовательно, упругая энергия деформации элемента может быть представлена следующим образом:
8ие =(5q'e)Т |[в'е]Т [ж'е] [в'е]йЬе4е (1.31)
Ье
или
Ъие — ^'е)Т[к'е]q'е , (1.32)
где [к ' е ] — матрица жесткости элемента.
Для получения выражения энергии упругой деформации конструкции необходимо провести суммирование выражений (1.32) для всех элементов, на которые разбита конструкция. Как правило, в методе конечных элементов матрица элементов строится в локальной системе координат, связанной с элементом. После чего умножением на матрицу поворота осуществляется переход в подвижную систему координат.
При построении расчетной модели конструкции разбивка на конечные элементы осуществляется таким образом, чтобы сосредоточенные усилия были приложены в узловых точках. Допустим, что в упругом теле задан конечный набор к узлов, в которых приложены силы Ру (= 1, 2, ..., к). Перепишем уравнения (1.20) и (1.21), подставив в них соотношения (1.28)— (1.30), заменив упругие деформации и на их приближенное представление через узловые перемещения q и добавив узловые сосредоточенные силы Ру :
к
2 Р} +| Рой о +1РГЛ Г - тЯ0 - [С] | [N'] йт q ' +
у=1 о г о
+| [р] йтш - [С] | [ш '][С]^' ] йт q' - [ю][ю] | рйт -
о о о
[ш][ш][С]|[N'] йтq'-2[ш][С]|[N'] йтq' = 0; (1.33)
о о
-2[Ру ](Р+ [С ] 4) у-|[Ро](р+ [С][ N' ] q') й о} =1 о
/РГ ](р + [С] [ N ] q') й Г + [Я01|(Р + [С] [N'] q') йт -
Г Г
- | р [С] [N'] йт^г ' -1 [и] ийт+| [р + и][р + и] йти +
о о о
+ [и] |[р + и] [р + и] йти — 2^[р + и] [и] ийт — 0. (1.34)
о о
Для получения конечномерного аналога системы уравнений (1.27) воспользуемся соотношением принципа возможных перемещений (1.19). Уравнение (1.19) должно выполняться для любого кинематически допустимого поля возможных упругих перемещений 5и. Потребуем выполнения этого соотношения для конечного набора функций, задаваемых матрицами функций формы [ N ]. Используя соотношения (1.28)—(1.32), получим систему обыкновенных дифференциальных уравнений
к
+ Гг л г'чТ гу л о I Гглт^То^лг ГглтТ
2 Р} + | [N']Т Ройо +1 [N']Т РгйГ -1 [N']Т Р'оайо -
}—1 о Г о
-1 [N']Т РГйй Г -1 [N']Т (Я0 + [ш '] р ' + [ш '][ш '] р ') йт -
Г о
к
- | [N']Т ([ш '] ^' ] + [ш '][ш'Ш' ])йтщ' - 21N']Т [ш ']^'] йтц -
о о
- I[N']Т[N'] йЩ' = І[В]Т [&][Б]й^'.
О О
(1.35)
Система уравнений (1.35) записана для всего ансамбля конечных элементов, операция суммирования по элементам расчетной модели опущена.
Таким образом, анализ динамики упругой конструкции сводится к решению системы нелинейных уравнений (1.33), (1.34), (1.35) совместно с кинематическими уравнениями, что позволяет определить Ло(0> ю(О, Ч (0 и А;(ґ) (і = 0, 1, 2, 3). Заметим, что область применимости полученной системы уравнений ограничена рамками малых деформаций упругих конструкций, но допускает быстрые пространственные развороты.
2. Анализ системы уравнений динамики свободной упругой конструкции. Уравнения движения свободного твердого тела. Рассмотрим твердое тело, занимающее область О в инерциальной системе координат 0*^у2^3. Свяжем с телом систему координат 0x^2Х3, ориентация
которой по отношению к неподвижной системе координат задается матрицей направляющих косинусов [С]. Движение твердого тела характеризуется поступательным движением и вращением связанной с телом системы координат. Вектор ю задает угловую скорость вращения тела. Как и в предыдущей части, вектором Ло будем характеризовать положение начала связанной системы координат относительно неподвижной. Введем вектор с, связывающий начало подвижной системы координат с центром масс тела. Следовательно, вектор Я, соединяющий начало неподвижной системы координат с центром масс тела
Л (0 = Лэ(*)+ с (О-
Введем вектор количества движения Q:
Q = тЛ
и, воспользовавшись теоремой об изменении количества движения, получим
dQ й (тЯ)
йґ йґ
где Р — главный вектор сил, приложенных к телу. Принимая во внимание, что с = [ю] с, получим
тЛо -т[с]ю + т[ю][ю]с = Р. (2.1)
Введем вектор кинетического момента тела К относительно начала неподвижной системы координат. Вектор з, имеющий начало в центре масс тела, пробегает все точки твердого тела, и тогда, по определению,
К = |[Л + ж] (Л + ж)йт = [Л] Лт + [Л]|жйт + |[ж] йтЛ + |[ж] жйт. (2.2)
О ООО
Так как вектор з берет начало в центре масс тела, то, по определению центра масс, выполняется следующее соотношение
| [ ж] йт = 0. (2.3)
О
Следовательно, второе и третье слагаемые правой части равенства (2.1) равны нулю, и выражение для кинетического момента принимает следующий вид:
К = [ Л] Лт + | [ ж] жйт.
О
(2.4)
Применим теорему об изменении кинетического момента
^ = В* + *1 ] Р , (2-5)
М
г
где р — внешние силы, приложенные к телу в точках, положение которых задается вектором . Выполняя дифференцирование по времени в (2.5), получим
(
[ Я] Яш + [ Я]
Яш - I р 1+| [ ж] ж Мш + I [ я ж Мш = £[*/ ]Р - (2-6)
V г / а а г
В соотношении (2.6) второе слагаемое в левой части равно нулю на основании теоремы о
движении центра масс системы. Первое и третье слагаемые равны нулю в силу тождества [а]а = 0,
справедливого для любого вектора а. Таким образом получим
|[ж] жйш = 2[(, ] Р . (2.7)
а г
Выполняя дифференцирование по времени в левой части (2.7), получим
I и ж йш = | [*][«]* йш + I [ж][ю][ю] ж йш =
а а а
= -|[ж][ж] Мшш - [га] |[ж][ж] Мшга. (2.8)
а а
Принимая во внимание, что по определению
[ ^ ] = -| [*][*] Мш,
а
где [/] — матрица инерции тела относительно центра масс, из (2.7) по-
лучим
1 +
[/] га + [га][/] га = £[*,■ ] р . (2.9)
Матрицу инерции тела относительно начала связанной системы координат [У*] можно выразить следующим образом:
[/ * ]= -1 [с + ж][с + ж] йш = - ш[с][с] - [с] | [ж] йш -
О О
| [ж] йш [с] -1 [ж][ж] йш = [ / ] - ш [с][с]. (2.10)
О О
Подставляя (2.10) в (2.9), с учетом (2.1) получим
ш [с] Яо +[/*] га = 2 [4 ]Рі + [с] р - [га] [/*] га. (2.11)
Система уравнений (2.1), (2.11) описывает движение свободного твердого тела под действием системы внешних сил. Сравним системы уравнений (1.20) и (2.1).
Первые два слагаемых уравнения (1.20) представляют собой внешние силовые воздействия на конструкцию. Им соответствует правая часть уравнения (2.1). Последнее слагаемое в системе
(1.20) — кориолисова сила, возникающая вследствие упругих деформаций тела. Первому слагаемому левой части системы (2.1) соответствуют третье и четвертое слагаемые системы (1.20), учитывающие упругие деформации. Аналогично можно поставить в соответствие оставшиеся
члены систем уравнений (1.20) и (2.1). Таким образом, система уравнений (1.20) отличается от системы уравнений (2.1) лишь дополнительным силовым воздействием — силой Кориолиса, а оставшиеся члены системы уравнений, описывающей движение твердого тела, аналогичны
есть матрица инерции деформированной упругой конструкции, можно провести аналогичные вышеприведенным рассуждения и установить идентичность систем уравнений. Различия состоят в том, что к силовым воздействиям добавляется момент кориолисовых сил (последний член системы уравнений (1.21), и в том, что все величины, входящие в систему (1.21), содержат поправки на дополнительную упругую деформацию. Система уравнений (1.33), полученная из
(1.20) переходом к конечномерному представлению, описывает пространственное перемещение начала связанной системы координат. Система (1.34), полученная из (1.21), описывает вращение связанной системы координат. Без учета кориолисовых сил, возникающих вследствие упругих деформаций, эти системы в каждый момент времени совпадают с уравнениями движения твердого тела, получившего
в данный момент упругие деформации и «затвердевшего».
Построение матриц конечных элементов свободной упругой конструкции. Система уравнений (1.35) описывает упругие деформации конструкции, величина которых зависит от внешних силовых воздействий и движения связанной системы координат (векторов В0(Ґ) и ю (0).
Систему уравнений (1.35) можно записать в виде
является вектором внешних силовых воздействий: сосредоточенных сил (отнесенных к узлам расчетной модели МКЭ), объемных сил, поверхностных сил.
При движении упругой конструкции в ее элементах возникают деформации, порождающие внутренние силовые факторы, к которым необходимо добавить силы демпфирования Род (/) и
Рш (/). Природа внутреннего демпфирования в упругом теле трудно объяснима и является
предметом интенсивных исследований в настоящее время. Эти силы порождаются как трением между элементами конструкции, так и потерей энергии внутри деформирующегося конструкционного материала. Ввиду того, что описание реального механизма потери энергии внутри деформируемого конструкционного материала вызывает определенные трудности, при проведении практических расчетов широко используются разнообразные упрощенные модели демпфирования, ориентированные, как правило, на организацию простой вычислительной процедуры.
Удовлетворительные результаты получаются при использовании эвристической модели демпфирования, величина которого пропорциональна скорости демпфирования упругого тела. Введя коэффициенты, характеризующие величину распределенных объемных сил демпфирования и поверхностных сил демпфирования , структуру вектора сил демпфирования (0 можно представить следующим образом:
соответствующим
членам
системы
(1.20)
с учетом поправок, возникающих вследствие упругих деформаций конструкции. Сравним системы уравнений (1.21) и (2.11). Принимая во внимание, что
О
(2.12)
і=1
О
О
/ Л
РЬ (0 = ГіО^'] Т №] й Т[№] йГ д' = [Ь] д',
ІО О /
где [О' ] — матрица демпфирования. Значения коэффициентов цо и определяются из эксперимента.
Вектор Г0(О является вектором сил, которые возникают вследствие ускоренного поступательного и вращательного движений связанной системы координат 0x^x2 Х3
Г0 (г) = | [N']Т (Я0 + [ш '] р ' + [ш '][ш'] р ' ) ат .
О
Даже при отсутствии внешних силовых воздействий существуют силовые факторы, деформирующие упругое тело вследствие его движения (вра-
щения).
Рассмотрим структуру вектора Г} (г)
Г' (г) = | [ N ' ]Т + ([ш ' ]^'] + [ш ' ][ш ' ]^ ' ])йтц' + 21 [N ' ] Т[ш ' ][N' ] йтц' +
О О
+ | [ N' ]Т [ N' ] йтц'.
О
Вектор сил инерции Г \ (г) можно записать в виде
Г} (0 = [О '+в ] ч'+[О' ] 4'+[М ' ] 4'. (2.13)
Силовые воздействия, пропорциональные вектору узловых смещений #(0, вычисляются с помощью матриц [С'] и [<2'], которые можно назвать матрицами псевдожесткости. Эти матрицы возникают вследствие вращения связанной системы координат.
Матрица [С' ] — кососимметрическая матрица
[О ' ] = | [ N' ] Т[ш ' ]^ ' ] а.т,
О
так как
[О ']Т =| [N']Т [ш '][N'] а.т = —| [N']Т [ш ']^ '] а.т = — [О '].
О О
Выражение для матрицы псевдожесткости [<2' ] имеет вид
[в'] = | [N'] Т[ш '][ш 'Ш'] ат.
О
Легко проверяется, что [<2' ] — симметрическая матрица (используется свойство [ш']Т =— [ш ' ]). Матрица [О' ], имеющая следующую структуру
[О ' ] = 21 [ N ]Т [ш ' ]^' ] а.т,
О
которая характеризует воздействие кориолисовых сил на упругую конструкцию. Эта матрица в представлении вектора Г (г) входит в виде произведения на вектор скорости упругих перемещений </' (г), следовательно, ее можно интерпретировать как матрицу псевдодемпфирования. Следует заметить, что в отличие от демпфирования наличие этого члена в выражении для Г} (г) не приводит к диссипации энергии системы. Матрица [О'] — кососимметрическая. Это свойство проверяется аналогично свойству кососимметричности матрицы [С' ].
Матрица [М' ], входящая в (2.13), является матрицей масс
[М'] = | [ N']Т [ N ' ] ат.
О
Вектор Г'(г), входящий в правую часть (2.12), является вектором
упругих сил
Г'(г) = |[в']Т [Ж '][£ ']йОд' = [К']4(г).
О
Матрицы [М' ] и [К ' ] — симметрические матрицы.
Второе и третье слагаемые вектора Г0(г), первые два слагаемых вектора Г} (г) являются векторами сил, возникающих вследствие вращения связанной системы координат.
Таким образом, используя введенные матрицы, соотношение (2.12) можно записать в виде
[М' ] 4' + [Б ' + О' ] 4 ' + [К ' + О ' + в ] 4 ' = Г ' (г) — (г). (2.14)
Система уравнений (2.14) для замкнутости должна быть дополнена начальными условиями
4 1/=0 = "о, 4 |/=о = Ь0. (2.15)
Система обыкновенных дифференциальных уравнений (2.14), (2.15) содержит столько уравнений, сколько степеней свободы используется в конечноэлементной модели упругой конструкции. В отличие от традиционных задач динамики упругих конструкций, которые исследуются с помощью метода конечных элементов, система уравнений (2.14), (2.15) содержит несимметрические матрицы, что в значительной степени усложняет построение вычислительного алгоритма.
Таким образом, для анализа динамического поведения свободных упругих конструкций требуется интегрирование систем нелинейных уравнений (1.33), (1.34),
(2.14), (2.15) для определения -^О(г), ю (г), (г) (5 =
= 0, ... , 3), д '(0. Очевидно, что интегрирование столь сложной системы нелинейных
дифференциальных уравнений требует разработки специализированных вычислительных
процедур.
В зависимости от типа исследуемой конструкции и характера ее движения возможны некоторые упрощения приведенной выше системы уравнений. Так, например, в случае медленных движений упругой конструкции при слабых внешних воздействиях уравнения (1.33), (1.34) можно заменить на (2.1), (2.11), т. е. исследовать пространственное движение и вращение связанной системы координат в предположении о малом влиянии упругих деформаций на характер этого движения. Однако даже такое сильное
упрощение уравнений тем не менее сохраняет необходимость интегрирования системы (2.14),
(2.15).
Рассмотрим некоторые частные случаи движения упругой свободной конструкции, допускающие упрощение системы (2.14), (2.15).
1) При установившемся движении конструкции (ш' = 0) матрица [С' ] становится нулевой
[М' ] 4'+[Б ' + О' ] 4' + [К ' + в ] 4 = Г ' (г)—Г0 (г).
Получаем систему уравнений демпфируемой гироскопической системы. Матрицы [М' ] и [К ' + Q' ]
— симметрические, [О' + О' ] — несимметрическая. Если пренебречь демпфированием ([О' ] = 0), получим систему уравнений движения недемпфируемой гироскопической системы ([О ] — кососимметрическая матрица).
2) При отсутствии вращения связанной системы координат получим систему уравнений
[М' ] 4' + [Б ] 4' + [К ' ] 4' = Г' (г)—Г0 (г)
с симметрическими матрицами [М'], [О' ] и [К']. Исследование подобных систем уравнений значительно упрощается и подробно описано в [3], [10].
3) Наиболее простой вид система (2.14), (2.15) приобретает в случае, кода можно пренебречь вращением и ускорением подвижной системы координат, а также демпфированием
[М' ]4' + [К' ] 4 ' = Г ' (г)
Подобный случай движения возможен лишь при определенной структуре вектора узловых сил Р '(О (система внешних сил должна быть самоуравновешенной, момент внешних сил должен быть равен нулю).
3. Выбор связанной системы координат. Для описания упругих деформаций конструкции используются уравнения линейной теории упругости или строительной механики, которые справедливы лишь при малых деформациях конструкции. Свободная конструкция может совершать свободные пространственные перемещения и развороты. Если описывать движение упругой конструкции в инерциальной (неподвижной) системе координат, то смещение точек конструкции будет складываться из смещений конструкции как твердого тела (больших смещений) и упругих перемещений, которые определяются как интегралы малых деформаций конструкции. Если размеры упругой конструкции не очень велики, то вследствие малости деформаций будут малы и упругие перемещения элементов конструкции по отношению к недеформированному состоянию. Поэтому, как правило, со свободной конструкцией связывается некоторая система координат [9], в которой отсчитываются упругие деформации, а большие перемещения и развороты конструкции как твердого тела описываются как движение связанной системы координат. Свободное твердое тело при пространственном движении имеет шесть степеней свободы — три смещения и три поворота вокруг координатных осей. Поэтому необходимо шесть соотношений, определяющих положение связанной системы координат по отношению к упругой конструкции. Эти шесть соотношений определяют способ введения связанной системы координат. Общее количество степеней свободы конструкции от выбора связанной системы координат не зависит, так что добавление шести уравнений, характеризующих движение связанной системы координат, компенсируется шестью уравнениями связи.
Шесть дополнительных соотношений должны однозначно определять положение связанной системы координат (начала связанной системы и ориентацию ее осей) по отношению к упругой свободной конструкции. Эти шесть уравнений связи могут иметь вид дифференциальных, интегральных или конечных соотношений. Выбор того или иного способа введения связанной системы координат определяется характером рассматриваемой динамической задачи, и во многих случаях носит чисто субъективный характер. В некоторых случаях за счет рационального способа введения связанной системы координат удается достичь значительного упрощения уравнений динамики свободной упругой конструкции и, тем самым значительно сократить вычислительные затраты.
Известен ряд способов задания связанной системы координат [4], [5]: это локально связанная система координат, главная система, система Тиссерана, средняя система. Кратко остановимся на способах задания этих систем координат.
Первый способ — это локально связанная система координат. Если свободная конструкция может быть условно представлена в виде недеформируемого несущего тела и соединенных с ним гибких носимых тел, то оси локально связанной системы координат фиксируются относительно центрального несущего тела [4]. Введение такой системы координат удобно при анализе динамики летательных аппаратов, когда начало и ориентацию осей полезно связать с приборным отсеком, в котором расположена навигационная аппаратура. Однако не всегда удается в реальной конструкции найти ее недеформируемый элемент, с которым было бы возможно связать систему координат. При построении математической модели свободной конструкции при использовании метода конечных элементов такой проблемы не возникает.
Если для описания упругого поведения конструкции используются уравнения строительной механики и соответствующие конечные элементы (балка, пластинки и т. д.), то следует выбрать узел с шестью степенями свободы и исключить их из расчетной модели. Начало связанной системы координат помещается в этот узел, и при выбранной ориентации осей строится конечноэлементная расчетная модель конструкции. При таком способе введения связанной системы координат шесть уравнений связи,
о которых речь шла выше, учитываются вычеркиванием соответствующих строк и столбцов в матрицах системы уравнений (2.14), (2.15). В случае, если при анализе деформаций конструкции используются уравнения линейной теории упругости, т. е. в узлах расчетной модели допускаются три степени свободы — три смещения, то при введении системы координат поступают следующим образом. Выбираются три узла расчетной модели конструкции, не лежащие на одной
прямой. В одном из узлов размещается начало связанной системы координат, т. е. запрещаются упругие смещения в этом узле. Ось 0x2 перпендикулярна оси Охі, а ось Охз перпендикулярна плоскости 0х\х2 и образует правую систему координат. При таком введении осей связанной системы координат запрещается упругое смещение вдоль оси Охз во втором и третьем узлах, а также смещение в направлении оси 0х2 во втором узле, т. е. выполняется шесть дополнительных соотношений, которые также сводятся к вычеркиванию соответствующих строк и столбцов в матрицах соотношений (2.14), (2,15).
При использовании МКЭ для анализа упругого поведения свободной конструкции описанный способ введения локально связанной системы координат является наиболее простым и естественным.
Следующая группа способов введения связанной системы координат отличается от описанного выше тем, что положение начала и ориентация осей не связывается жестко с упругой конструкцией. Такие системы координат называются «плавающими». При движении и деформировании упругой конструкции система координат движется вместе с ней, смещаясь относительно ее элементов тем или иным образом. Смысл введения плавающих систем координат состоит в том, что посредством выбора соотно-
шений, определяющих положение начала и ориентацию осей системы, можно несколько упростить уравнения движения упругой конструкции (1.33)—(1.35).
Так, например, в работе [6] обсуждается способ введения так называемой главной системы координат. Начало связанной системы координат помещается в центре масс конструкции, а ориентация осей системы выбирается в направлении главных осей инерции конструкции. При такой ориентации связанной системы координат матрица инерции конструкции становится диагональной, т. е.
что, в силу симметрии матрицы инерции [7] дает три соотношения связи. Еще три соотношения связи следуют из того, что начало связанной подвижной системы координат должно располагаться в центре масс конструкции
О
Третий способ введения связанной подвижной системы координат описан Тиссераном [6]. Положение начала и ориентация осей этой системы определяется из условия равенства нулю векторов относительного количества движения и относительного кинетического момента
Два векторных соотношения (3.3) дают шесть уравнений связи. Эти соотношения носят нелинейный характер, и их выполнение в процессе интегрирования по времени системы уравнений (1.33)—(1.35) порождает определенные трудности.
Линеаризация соотношений (3.3) предложена в работе [7]. Положение начала и ориентация связанной системы координат определяются из условий
Введенная таким образом система координат называется средней. В средней системе координат первое векторное равенство (3.4) указывает на то, что ее начало совпадает с центром масс конструкции. Средняя система координат часто используется при анализе динамики свободных упругих конструкций в силу того, что ее использование во многих случаях позволяет существенно упростить вычислительную процедуру. Эффективность применения средней системы координат подробно обсуждается в работе [8].
Работа поддержана грантом Российского фонда фундаментальных исследований № 98 -01-
= 0, і * І (і, І = 1, 2, 3),
(3.1)
(3.2)
(3.3)
О
О
(3.4)
О
О
ЛИТЕРАТУРА
1. Бранец В. Н., Шмыглевский И. П. Применение кватернионов в задачах ориентации твердого тела.— М.: Наука.— 1973.
2. Зенкевич О. Метод конечных элементов в технике.— М.: Мир.— 1975.
3. Бате К., Вилсон Е. Численные методы анализа и метод конеч-
ных элементов.— М.: Стройиздат.— 1982.
4. Лурье А. И. Аналитическая механика.— М.: Физматгиз.— 1961.
5. Canavin J. R., Likins P. W. Flotating reference frames for flexible spacecraft//! Spacecraft and
Rockets.— 1977. V. 14, N 12.
6. Докучаев Л. В. Нелинейная динамика летательных аппаратов с деформируемыми элементами.— М.: Машиностроение.— 1987.
7. М а р к е е в А. П. К динамике упругого тела в гравитационном поле//Космические исследования.— 1989. Т. 27. Вып. 2.
8. Баничук Н. В., Карпов И. И., Климов Д. М., Марке -
ев А. П., Соколов Б. Н., Шаранюк А. В. Механика больших космических конструкций.— М.: Изд-во «Факториал».— 1997.
9. Шиманский Ю. А. Динамический расчет судовых конструкций.— Л.: Судпромгиз.— 1964.
10. Argyris J., Mlejnek H. P. Dynamics of Structures. North Hol-
land.— 1991.
Рукопись поступила 27/IX1999 г.
