Эффективные методы численного моделирования динамики нелинейных систем абсолютно твёрдых и деформируемых тел (1105385), страница 12
Текст из файла (страница 12)
A(φ) – этоматрица поворота тела⎡cos ϕA(ϕ ) = ⎢⎣ sin ϕ− sin ϕ ⎤.cos ϕ ⎥⎦(3.6)Таким образом, после того, как мы присоединили тело к балке, мы получим дополнительно (в скалярной форме) три линейных дифференциальныхуравнения (3.3) и три нелинейных уравнения связей (3.4), (3.5). Уравнениядвижения полученной гибриднЮй системы могут быть записаны в виде системы дифференциально-алгебраических уравнений [51, 19, 15] индекса 1⎡M e⎢⎢ 0⎢ Ge⎣0MrGrG e T ⎤ ⎧&e& ⎫ ⎧Q g − Q e ⎫⎥⎪ ⎪ ⎪ g⎪&& ⎬ = ⎨ f − f i ⎬ ,G r T ⎥ ⎨u0 ⎥⎦ ⎪⎩ λ ⎪⎭ ⎪⎩ h ′′ ⎪⎭где λ – столбец множителей Лагранжа, Ge, Gr – матрицы Якоби уравненийсвязи по координатам e и u.Есть возможность уменьшить число дополнительных уравнений путёмзаписи уравнений (3.3) относительно точки прикрепления O.
Тогда матрицыи векторы в этих уравнениях становятся непостоянными:~⎡ mII A ρC ⎤⎧r ⎫ru = ⎨ ⎬ , M (ϕ ) = ⎢ ~,T2⎥ϕ+(IAρ)Jmρ⎩ ⎭CCC⎦⎣⎧− ϕ& 2 A ρC ⎫⎧m g⎫i&f (ϕ ) = ⎨ T ~⎬.⎬ , f (ϕ , ϕ ) = ⎨0⎩m g I A ρ C ⎭⎩⎭gВ этом случае уравнение связи (3.4) превращается в rℓ – r = 0 и можетбыть исключено с применением процедуры сборки, см. приложение 6.3.
Вэтом случае получим дополнительно одно нелинейное дифференциальноеуравнение (после сборки) и одно нелинейное уравнение связи (3.5).Как и следовало ожидать, использование угла поворота тела в качествеобобщённой координаты в формализме абсолютных узловых координат приводит к появлению ДАУ, поскольку в указанном формализме используются76компоненты касательных векторов. ДАУ вносят дополнительные трудности впроцесс интегрирования, например, возникают проблемы «ухода» системысо связей и другие.
Для решения ДАУ разработаны, тем не менее, эффективные численные методы [25, 46, 42, 39].Но несомненно, что если можно избежать появления ДАУ и свести задачу к ОДУ, то такой подход предпочтительнее, поскольку он менее трудоёмок. В нашем случае такая возможность существует.3.1.4.2.Использование абсолютных узловых координат в качествеобобщённыхКак показано выше, если угол поворота используется в качестве обобщённой координаты для тела, присоединённого к концу балки, то возникаетнеобходимость добавления уравнений связи.
Естественно отказаться от этогои использовать для абсолютно твёрдого тела тот же набор узловых координат, что и для концевого узла балки:x = {r T , τ T } T .(3.7)Этот подход позволяет избежать появления уравнений связей. Покажемэто и выведем уравнения движения для тела, показанного на Рис. 3.7.yvЗащемлениеωТелоA ρCCAρτЭлемент балкиxOϕГоризонталь-rная осьРис. 3.7. Кинематика тела, присоединённого к балкеАбсолютная скорость произвольной точки, заданной радиус-вектором ρ,определяется выражением v = r& + ω × A ρ .
В двумерном случае можно запи-77сать это соотношение в виде~v = r& + ω I A ρ~с использованием кососимметричной матрицы I из уравнения (2.11) и матрица поворота A (выражение (3.6)).Определим угол поворота тела следующим образом: ϕ = arctg(τ 2 τ 1 ) .Угловая скорость тела ω – это производная от угла поворота~τ&2τ 1 − τ&1τ 2 τ& T I τω = ϕ& = 2 2 = 2 , τ = τ = τ 12 + τ 22 .τ1 + τ 2τТогда вектор v скорости произвольной точки представим в виде v = Φ x& ,~∂ω I τ~T.где введена матрица Якоби Φ = I I A ρ ω τ& , здесь ω τ& ==∂τ& τ 2[]Теперь мы готовы записать уравнения движения тела.
Проще всего этосделать на основе общего уравнения динамики∫V δrTµ (a − g ) dV = 0 ,(3.8)где введены векторы возможного перемещения δr = Φ δx и ускорения& x& произвольной точки тела. g – это вектор ускорения силыa = v& = Φ &x& + Φтяжести, µ – плотность материала в кг/м3. Интегрирование по объёму V тела.После довольно трудоёмких вычислений мы найдём, что⎧ ω 2A ρ ⎫⎪& x& = −Φ ⎪⎨ 2 TΦ⎬.&&()τττ⎪⎭⎪⎩τ 2Поскольку вариации обобщённых координат δx произвольные, на нихможно «сократить» и уравнение (3.8) превращается в∫V µ ΦилиT& x& − g ) dV = 0 ,(Φ &x& + ΦM &x& + f i = f g ,(3.9)с матрицей масс M и векторами обобщённых сил инерции f i и тяжести f g:~⎡mIm I A ρC ωTτ& ⎤T,M( x ) = ∫ µ Φ Φ dV = ⎢~TT ⎥mω(IAρ)Jωωτ&CO τ& τ& ⎦⎣V78⎧ ω 2 A ρC ⎫⎪& dV x& = − M ⎪⎨ 2f ( x, x& ) = ∫ µ Φ Φ⎬,T&&τττ()⎪⎩τ 2⎪⎭Vmg⎧⎫f g ( x ) = ∫ µ ΦT dV g = ⎨~⎬.TωgIAρm&⎩⎭τCViTВ этих выражениях m – масса тела, ρC – радиус-вектор центра масс тела,J O = J C + mρ C2 – момент инерции тела относительно точки прикрепления O,JC – момент инерции тела относительно центра масс.Заметим, что все члены уравнений (3.9) нелинейны.
Но это не являетсянедостатком описанного метода, поскольку уравнения движения данной системы будут нелинейными и в других формулировках, указанных в предыдущем пункте. Основное преимущество предложенного метода заключается втом, что он приводит к обыкновенным дифференциальным уравнениям, а нек дифференциально-алгебраическим.3.1.5.Сравнение экспериментальных данных и расчёта3.1.5.1.Сравнение частот малых колебанийЦелью данного исследовния было сравнение частот малых колебанийконсольной балки, полученных с использованием конечно-элементой расчётной схемы, с аналитическими решениями для этой простой задачи, а также сданными, полученными из экспериментов (см. пункт 3.1.2).Постановка проблемы собственных колебаний описана в п. 2.5.3.Табл.
3.4. Расчёт частоты малых свободных колебаний консольной балкиБалка № 1 (d = 1,0 мм)Частота, рад/сЭлементыω2ω1127,740 273,31227,622 174,49427,610 173,22527,609 173,11827,609 173,04∞ (точно) 27,609 173,02Балка № 2 (d = 0,6 мм)Частота, рад/сЭлементыω2ω1116,604 163,60216,534 104,44416,526 103,69516,526 103,62816,526 103,57∞ (точно) 16,526 103,56Как видно из таблицы, рассчитываемые значения первых частот при измельчении конечноэлементой сетки стремятся к точным значениям.
Это свидетельствует о корректности модели конечного элемента.3.1.5.2.Сходимость результатов численного моделированияОдним из критериев, определяющих корректность расчётов по методуконечных элементов, является сходимость результатов вычислений при измельчении сетки конечных элементов. Результаты подобного теста на сходимость для балки диаметром 1 мм, с грузом 20 г приведены на Рис. 3.8. На нёмпоказаны графики вертикального перемещения конца балки при различномчисле конечных элементов: 5, 10 и 15.0.00.51.0 время, с 1.52.02.53.00перемещение, мм-505 КЭ-100-150-2 00-25010 и 15 КЭ-300-350Рис. 3.8.
Вертикальное перемещение конца балки во времени при различномчисле конечных элементов. Балка: d = 1 мм; груз m = 20 г.Графики очень близки, и при числе элементов, большем 10, практическисовпадают. Поэтому все последующие численные эксперименты выполнялись при числе элементов 10.3.1.5.3.Учёт затухания колебанийДля моделирования эффектов внутреннего трения в материале и сопротивления воздуха использовалась линейная модель сил демпфированияQ d = D e& .Применяя классическую модель демпфирования по Рэлею [19, 15], мат-80рицу демпфирования D определим формулойD = α M + β C,(3.10)в которую входят матрица масс M и касательная матрица жёсткости C (матрица Якоби сил) с коэффициентамиα=2ω1ω 2 (ζ 1ω 2 − ζ 2ω1 ),ω 22 − ω12β=2 (ζ 2ω 2 − ζ 1ω1 ).ω 22 − ω12(3.11)Они зависят от частот ω1, ω2 и параметров демпфирования ζ1, ζ2 первых двухформ колебаний и входят в уравнения движения в главных координатах&x&i + 2ζ iω i x&i + ω i2 xi = 0 , i = 1,2.Величины ζ1 и ζ2 являются долями демпфирования от критического идолжны быть определены перед началом каждого численного моделированияиз экспериментальных данных.
Формула для вычисления ζ1 имеет вид [2]δζ1 = 12π⎛ δ1 ⎞⎜1 +⎟⎝ 2π ⎠−12≈δ1.2π(3.12)где δ1 – логарифмический декремент затухания. Приближённое выражениесправедливо в случае малого затухания колебаний, как в нашем случае.Логарифмический декре-xA1мент δ1 можно найти по экспе-A2Anриментальной последовательности амплитуд A1,…,An (см.Рис. 3.9), используя простуюtOформулуРис. 3.9. Затухающие колебания1n −1(3.13).δ 1 = − ln ∆1 , ∆1 ≈ ⎛⎜ An A ⎞⎟1⎠⎝Таким образом, алгоритм нахождения коэффициентов α и β следующий:1) Вычислить или измерить первую и вторую частоту ω1, ω2;2) Определить величину декремента колебаний ∆1 и логарифмическийдекремент δ1 по экспериментальным осциллограммам и формулам (3.13);3) Вычислить коэффициент демпфирования ζ1, используя соотношение81(3.12); мы также полагаем, что ζ1 = ζ2;4) Наконец, получить значения α и β по формулам (3.11).3.1.5.4.Колебания свободной балки без грузаПроведём вычисления согласно описанному алгоритму для первого эксперимента: балка d = 1 мм, масса груза (бумажной метки) m0 = 0,023 г.1) ω1 = 27,6 с–1, ω2 = 173 с–1; 2) ∆1 = 0,987; 3) ζ1 = ζ2 = 0,002;4) α = 0,1 с–1, β = 2·10–5 с.Можно заметить, что величина β намного меньше α.
Поэтому естественно попытаться отбросить часть сил демпфирования (3.10), пропорциональную матрице жёсткости и использовать матрицу демпфированияD = α M.(3.14)Численное интегрирование показало, что результаты, полученные с использованием двух моделей сил (3.10) и (3.14) действительно различаютсятолько в 4-5-й значащей цифре. Однако шаг интегрирования в случае полноймодели матрицы демпфирования (3.10) примерно в 20 раз меньше, посколькууравнения движения получаются намного жёстче. Поэтому в дальнейшемиспользуется упрощённая матрица демпфирования (3.14).На Рис. 3.10 показаны графики вертикального перемещения конца балки, найденные из эксперимента (сплошная линия) и расчётные (пунктирная).0123 время, с 456780перемещение, мм-5-10-15-20-25-30-3 5-4 0-4 5Рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
