Эффективные методы численного моделирования динамики нелинейных систем абсолютно твёрдых и деформируемых тел (1105385), страница 8
Текст из файла (страница 8)
Она является почти столь же точной, как и модель L3, и почти столь же быстрой, какпростейшая модель L1, описываемая ниже.2.2.1.3.Модель L1Простейшая модель продольных сил получена в работе [17] с использованием следующего способа осреднения продольной деформации (Рис. 2.3):r0 − rl( e1 − e 3 ) T ( e 1 − e 3 )ε ≈−1 =− 1.llВ этом случае продольные силы являются линейными по e (если не учитывать зависимость ε от e) и принимают вид (2.9). Их матрица Якоби:⎧(e1 − e 3 ) ε l 2 ,T⎪4⎛⎞∂ε0∂ε⎪11 ⎜εε⎟ , гдеCij = K ij I + EA ∑ Sik e k=⎨⎜ ∂e ⎟∂e j ⎪ (e 3 − e1 ) ε l 2k =1⎝ j⎠⎪⎩0еслиеслиеслиеслиj = 1;j = 2;j = 3;j = 4.472.2.2.Модели обобщённых поперечных силОбобщённые силы, вызванные изгибной деформацией осевой линиибалки, вычисляются подобно продольным силам:∂П κ∂κdp .Qi == EJ ∫ κ∂e i∂ei0lκ(2.10)В этом выражении кривизна осевой линии κ зависит от первой и второйпроизводных r' и r″ радиус вектора r по дуговой координате p:κ=r ′ × r ′′r′3=r1′r2′′ − r2′ r1′′f3=⎡ 0 − 1⎤ ⎧ r1′ ⎫ 1 T ~1′′′′{}rr12 ⎢⎥ ⎨r ′ ⎬ = 3 r ′′ I r ′ ,+10f3⎣⎦⎩ 2⎭ f(2.11)~где f = r ′T r ′ – градиент деформации, а I – кососимметрическая матрица,значение которой приведено в формуле (2.11).
Далее приведен один из способов вычисления обобщённых сил.Модель T3. Кривизна осевой линии балки может быть выражена черезфункции форм и узловые координаты согласно формуле (2.1):κ=1 4 4 1 ⊗ T~S e Ien ,3 ∑ ∑ 2 mn mf m =1 n =1(2.12)⊗где введены антисимметричные символы Smn= sm′′ sn′ − sm′ sn′′ .Производные от кривизны вычисляются аналогично:4 4∂κ1 4 ⊗~1 4 ⊗~1 ⊗ T ~ ⎛ − 3⎞ 1 4= 3 ∑ Sik I e k + ∑ ∑ S mn e m I e n ⎜⎜ 4 ⎟⎟ ∑ si′sk′ e k ≈ 3 ∑ Sik I e k .∂e if =12f k =114n =424443 ⎝ f ⎠ 1k4243 f k =11442443 m1=14оставим лишь это(2.13)∂f ∂e if 3κПодстановка кривизны (2.12) и её градиентов (2.13) в (2.10) даётQκi =(4)4EJ 4 4 1 ⊗⊗ T ~~~SeIeIe=∑ f 6 ∑ ∑ 2 ikmn m n k ∑ K ikκ I e k ,k =1 14m4=1n4=14244443k =1κK ikгде f =l1r ′T r ′ dp =∫l0l1 4 4 11 T∑ ∑ Smn e m e n – усреднённая деформация, иl m =1 n =1⊗⊗⊗2211221122112211Sikmn= ∫ Sik⊗ S mndp = Simkn− S kmin− Sinkm+ S knim,0l2211= ∫ si′′sm′′ sk′ sn′ dp .Simkn0Явные значения этих символов приведены в приложении 6.7.48Эта модель сил является новой и названа T3, поскольку она кубичная покоординатам e.
Можно задать вопрос – а в чём состоят достоинства и недостатки этой модели по отношению к известным моделям [17, 37, 60]? С однойстороны, предложенная модель является более сложной, чем модель T2 вработе [17], но, с другой стороны, она имеет тот же порядок, что и описаннаявыше модель L3 (L2 в обозначениях работы [17]), и поэтому можно ожидать,что она будет более точной при совместном использовании с моделью L3.Приближённое значение матрицы Якоби от полученных сил имеет вид∂Qκi~ EJ 4 4 ⊗⊗ ~~κI e m ( I e n )T .C ij = T ≈ K ijκ I + 6 ∑ ∑ S imjn∂e jf m =1 n =12.3. ПРИМЕРЫ МОДЕЛИРОВАНИЯ БАЛОК И СРАВНЕНИЕРАЗЛИЧНЫХ ПОДХОДОВНиже приведены несколько решений тестовых задач с плоскими балками с целью сравнения результатов либо с известными аналитическими иличисленными решениями, а также с работами других авторов. В качестве инструмента для моделирования был использован программный комплекс«Универсальный механизм» (УМ), [42], www.umlab.ru. В составе комплексаавтором были реализованы конечные элементы, описанные выше.
Дальнейшая проверка корректности балочного элемента проведена в главе 3.2.3.1.Изгиб консольной балки сосредоточенной силойНа Рис. 2.4 изображена консольнаябалка, нагруженная сосредоточенной вертикальной силой P на свободном конце.Сила P вызывает большие перемещения концевого сечения: угол поворота θ,вертикальное δв и горизонтальное δг перемещения.Рис. 2.4. Изгиб консольной балки49В Табл. 2.1 дано сравнение решений этой задачи с использованием конечных углов поворота (п 1.2.2), а также формализма абсолютных узловыхкоординат (п. 1.2.3), и с точным решением при различном числе конечныхэлементов n. Точное, а вернее численное решение аналитической задачи обэластике, берётся из работы [11].Табл.
2.1. Большие перемещения торца консольной балкиPL2EJθπ 2ТочноерешениеδвδгLL0,25 0,079 0,083 0,0040,5 0,156 0,162 0,01610,294 0,302 0,05620,498 0,494 0,16050,774 0,714 0,388100,911 0,811 0,555Формализм конечныхуглов поворотаnθπ 231617112222101260,0790,0800,1560,1590,2940,3180,4980,5230,7740,7770,9140,921δвδгL0,0830,0830,1620,1670,3020,3330,4940,5210,7140,7170,8140,822L0,0040,0000,0160,0000,0560,0000,1600,1320,3880,3870,5540,551Формализм абсолютныхузловых координатNθπ 221314152611260,0790,0760,1560,1430,2940,2760,4980,4930,7740,8580,9110,911δвδгL0,0830,0780,1620,1380,3020,2390,4940,4790,7140,5900,8110,807L0,0040,0040,0160,0120,0560,0370,1610,1540,3880,2820,5550,553Как видно из таблицы, при использовании абсолютных координат требуется примерно в два раза меньше конечных элементов для получения одинаковой точности.
Это неудивительно, так как при использовании конечныхуглов поворота используется соотношения линейной теории (1.31), а в случае абсолютных координат применяется нелинейное соотношение (1.42). Ноесть и ещё более важное преимущество абсолютных координат. Уравненияравновесия в последнем случае получаются алгебраическими, в отличие оттрансцендентных уравнений в методе смешанных координат. Для устойчивости численного решения трансцендентных уравнений приходится применятьсхему релаксации, которая замедляет процесс решения.
Алгебраические жеуравнения решаются устойчиво и без релаксации.502.3.2.Сжатие консольной балки закритической силой спотерей устойчивостиНа Рис. 2.5 изображеныформы потери устойчивостиконсольной балки при приложении к её свободномуторцу сжимающей силы P,большей критической силыЭйлераPкр = π 2 EJ 4 L2 .Рис. 2.5. Потеря устойчивости балкиТабл. 2.2.
Углы поворота сечения консольной балки при потере устойчивостиP Pкр1,015 1,063 1,152 1,293θ˚точно2019,8619,6618,1714,93θθθабс . к − ты20 КЭабс . к − ты5 КЭотн . к − ты20 КЭθ 5отКЭн. к −ты4039,8939,5639,0926,436060,1459,7559,6452,971,5181,8842,5414,0299,1168010080,00 100,0879,61 99,7279,65 99,8475,53 97,37120120,06119,78119,92118,73140140,04139,85140,01139,94160160,03159,91160,09161,10180176,04175,90176,11177,60В Табл. 2.2 приведены значения углов θ поворота свободного торца, полученные двумя методами (точные решения этой задачи, а также решения вабсолютных координатах приведены также в работе [18]).2.3.3.Движение маятника в виде гибкой балкиТестовый пример движения маятника, составленного из 100 КЭ, показанна Рис. 2.6.
Он приведен для сравнения с данными из работы [18].Параметры модели: длина балки L = 1,2 м; плотность материалаµ = 5540 кг/м3; площадь поперечного сечения F = 18·10-4 м2; момент инерцииплощади сечения J = 1,215·10-8 м4; модуль Юнга E = 0,7·106 Па. Параметрыпроцесса численного интегрирования: шаг интегрирования h = 10-5 c, времяинтегрирования T = 600 с; процессор Pentium III, 650 МГц.51Рис. 2.6. Положения маятника во время движения2.3.4.Движение гибкой линейки эллипсографа с маятникомНа Рис.
2.7 приведен пример движения простейшей гибридной системы,состоящей из гибкой линейки эллипсографа (20 КЭ) и подвешенного к еёсередине маятника (твёрдое тело).Параметры:линейки:L = 1 м;µ = 7800 кг/м3;F = 10-4 м2;J = 10-8 м4;E = 108 Па;маятника: L1 = 0,5 м;m1 = 0,2 кг;J1 = 0,1 кг·м2;интегрирования: h = 10-4 с;T = 180 с.Рис. 2.7. Движение эллипсографа с маятникомВ задачах статики систем с большими перемещениями метод абсолютных координат имеет преимущества перед методом, использующим конечные углы поворота в качестве координат. При моделировании динамики нелинейных систем это преимущество сохраняется за счёт постоянства матрицмасс конечных элементов.522.4.
НОВЫЙ ПЛАСТИНЧАТЫЙ ЭЛЕМЕНТ НА ОСНОВЕОБОБЩЕНИЯ ФОРМАЛИЗМА АБСОЛЮТНЫХ УЗЛОВЫХКООРДИНАТКак уже упоминалось в п. 1.2.3, ранее была предложена реализация конечного элемента толстой пластины в формализме абсолютных координат(Шабáна и Мúккола [36]). В этой работе были использованы трёхмерныефункции форм, и их элемент был, фактически, массивным телом. Авторыиспользовали соотношения механики сплошной среды в трёхмерном случаедля вывода выражений для потенциальной энергии деформации пластины.Здесь рассматриваются несколько предложенных автором моделей элементов-пластин [74, 73]: прямоугольные с 48-ю и 36-ю, а также треугольныйэлемент с 27-ю степенями свободы.
Теория пластин Кирхгофа и нелинейныесоотношения между перемещениями и деформациями используются дляполучения обобщённых упругих сил, а также дифференциальная геометрияповерхностей для вычисления деформаций и кривизн.Предложенный элемент способен описывать большие движения и допускает большие относительные деформации.2.4.1.Узловые векторы и функции форм конечногоэлемента тонкой пластиныНачнём рассмотрение пластинчатых элементов с прямоугольного эрмитова элемента с размерами a, b, h (длина, ширина, толщина), см. Рис.
2.8.zyz3z4az1′1z ′xyOz ′y1bz2bz ′x1aРис. 2.8. Стандартный элемент пластиныx53В каждом из четырёх его узлов вводятся четыре степени свободы: например, мы имеем вертикальное перемещение z1 для узла 1, два наклоненияz ′x1 = (∂z ∂x )1иz ′y1 = (∂z ∂y )1 ,атакжевторуюпроизводную′ 1 = (∂ 2 z ∂x ∂y )1 . Поле перемещений пластины выражается через узловыеz ′xyперемещения и функции форм следующим образом [67]:z ( x, y ) = { S11 , S12 , S13 , S14 ; ...; ...; S 41 , S 42 , S 43 , S 44 } qс использованием функций Эрмита S ij ( x, y ) = si ( x ) s j ( y ) , построенных избалочных функций. q – это вектор узловых степеней свободы:′ 1 , z ′x 3 , z ′xy′ 3 , z 2 , z ′y 2 , z 4 , z ′y 4 , z ′x 2 , z ′xy′ 2 , z ′x 4 , z ′xy′ 4 }Tq = { z1 , z ′y1 , z3 , z ′y 3 , z ′x1 , z ′xyСледуя тем же путём, как и в пункте 2.1, мы заменяем поле перемещений z(x, y) тремя функциями x(p1, p2), y(p1, p2), z(p1, p2), описывающими параметризованную поверхность r(p1, p2):r =Se(2.14)с матрицей глобальных функций формS = [S11I, S12 I, S13I, S14 I; ...;...; S 41I, S 42 I, S 43I, S 44 I](2.15)S ij = si ( p1 , a ) s j ( p 2 , b) = sˆi sˆˆ j .Здесь введена единичная матрица I размером 3 × 3, а также балочные функции форм (ξ = p l)s1 ( p, l) = 1 − 3ξ 2 + 2ξ 3 ,s3 ( p, l) = 3ξ 2 − 2ξ 3 ,s2 ( p, l) = l (ξ − 2ξ 2 + ξ 3 ) , s4 ( p, l) = l (ξ 3 − ξ 2 ) .Вектор узловых координат, используемый в выражении (2.14) имеет вид{00 01 00 01 10 11 10 11 00 01 00 01 10 11 10 11e = r00r00 r0b r0b r00 r00 r0b r0b ra 0 ra 0 rab rab ra 0 ra 0 rab rab}T ,(2.16)ijij.
Элементы ruv– это векторыгде опущены знаки транспонирования над ruvijruv=∂i + jr∂p1i ∂p2jp1 = u ,p2 = v(2.17)являющиеся либо векторами узловых перемещений углов пластины (когдаi = j = 0), либо касательными векторами (когда i + j = 1), либо векторами вторых производных (когда i + j = 2). Их геометрический смысл ясен из Рис. 2.9.5411rab01rab00rab11r0b01r0bbra010p1ra110ara100ra00010r0b00r0bp210rab11r0001r0000r0010r00Рис. 2.9.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
