Эффективные методы численного моделирования динамики нелинейных систем абсолютно твёрдых и деформируемых тел (1105385), страница 7
Текст из файла (страница 7)
Осевая линия балкиПервое слагаемое r* представляет начальную конфигурацию осевой линии(прямую), а вектор u определяет перемещения балки по её длине. Соотношения (1.42) приводят к выражению для продольной деформации40ε = u1′ + 12 (u1′2 + u2′2 ) ,которое содержит как линейные, так и квадратичные члены.Поперечные силы вычисляются аналогично продольным:∂П κQi == EJ∂e iκl∂κ∫ κ ∂e i dp .(1.43)0Кривизна осевой линии балки вычисляется через первую r' и вторую r″производные от радиуса-вектора r по дуговой координате p:κ=r ′ × r ′′r′3=r1′ r2′′ − r2′ r1′′.( r1′2 + r2′2 ) 3 2Итак, формулы (1.41) и (1.43) дают способ вычисления обобщённых упругих сил.
Они довольно громоздки, но для их непосредственного вычисления используются разного рода упрощения и допущения. Различные моделиэтих сил были предложены в работах [17, 60, 69]. Все они имеют видQe = K(e) e, где K – нелинейная матрица жёсткости. Подробнее см. § 2.2.1.2.4.Другие модели балочных элементов, а также пластинКроме описанного конечного элемента балки Эйлера-Бернулли в плоской постановке были предложены и более сложные – например, элементплоской балки Тимошенко с учётом сдвиговых поперечных деформацийосевой линии [37], а также элемент балки в пространственной постановкеШабáны и Якýба [53], который изображён на Рис.
1.16. В качестве обобщёнОсевая линияr0;3x3rl;3r0;1r0r0;2x1 rl;1X3X1rlных координат они использовали ради-x2ус-векторы двух концов r0 и rℓ, а такжеrl;2касательные векторы rα;k = ∂rα/∂xk в этихточках,xk–криволинейныекоординатные линии, связанные с поперечнымиX2гдесечениямибалки.элемент имеет 24 степени свободы.Рис. 1.16. Элемент пространственнойбалки Шабáны и ЯкýбаЭтот41Миккола и Шабана предложили и прямоугольный элемент пластины впространстве [36]. Обобщёнными координатами, как и в случае балки, такжеслужат радиус-векторы rα четырёх его углов, α = 1,…,4, и касательныевекоры rα;k = ∂rα/∂xk, k = 1,…,3.
Элемент имеет 48 степеней свободы.Авторы используют трёхмерные функции форм, и поэтому, фактически,их элементы балки пластины представляют собой объёмные тела. Для построения выражения для потенциальной энергии используются соотношениямеханики сплошной среды, а вычисление обобщённых сил сводится к численному интегрированию по объёму элементов, что довольно трудоёмко.1.3.
ПЕРСПЕКТИВЫ РАЗВИТИЯ МЕТОДОВ МОДЕЛИРОВАНИЯВ заключение данной обзорной главы отметим основные недостаткисуществующих методов моделирования деформируемых тел. Использованиетвёрдотельных конечных элементов, по существу, хорошо отражает лишьстатику деформируемых тел – особенно это касается трёхмерных задач. Использование формализма конечных углов поворота, инкрементных методов идругих, связанных с введением подвижной системы координат, в случаебольших перемещений всей конструкции как твёрдого тела приводит к появлению уравнений движения, в которых все элементы отличаются сильнойнелинейностью. В этой связи перспективным, по мнению автора, являетсяметод абсолютных узловых координат, который свободен от указанных недостатков.
Поэтому глава 2 посвящена развитию этого формализма.422. РАЗВИТИЕ МЕТОДОВ МОДЕЛИРОВАНИЯДЕФОРМИРУЕМЫХ ТЕЛ2.1. НОВАЯ ТРАКТОВКА ФОРМАЛИЗМА АБСОЛЮТНЫХ УЗЛОВЫХКООРДИНАТ КАК ОБОБЩЕНИЯ МЕТОДА КОНЕЧНЫХЭЛЕМЕНТОВВ этой главе мы рассмотрим использование формализма абсолютныхузловых координат как обобщения узловых координат и полей перемещений,используемых при построении обычных конечных элементов. Сначала опишем его для двумерной балки, а затем обобщим для пластин в главе 2.4.1.Рассмотрим элемент балки длиной ℓ на Рис. 2.1, используемый в прило-yжениях метода конечных элементов [67].y = y (x)y0′двух его узлов предполагаются малыми.Под наклонами понимаются тангенсы уг-yly0OПеремещения y0, yℓ и наклоны y0′ , y l′yl′lлов наклона αχ осевой линии к оси x:xy ′χ = (dy dx ) x = χ = tg α χ .Рис.
2.1. Стандартный балочный элементПоле перемещений внутри элемента можно интерполировать в видеy ( x ) = s1 ( x ) y0 + s2 ( x ) y0′ + s3 ( x ) y l + s4 ( x ) y l′ ,где s1(x), …, s4(x) – функции Эрмита [9]:s1 ( p ) = s3 ( l − p ) = 1 − 3ξ 2 + 2ξ 3 ,s3 ( p ) = 3ξ 2 − 2ξ 3 ,s2 ( p ) = − s4 ( l − p ) = l (ξ − 2ξ 2 + ξ 3 ) , s4 ( p ) = l (ξ 3 − ξ 2 ) ,ξ=p.lЧтобы описать произвольные перемещения точек срединной линии, мыможем параметризовать её, используя дуговую координату p ∈ [0, l] и вводядва поля перемещений x(p) и y(p) вместо y(x), как показано на Рис.
2.2:x ( p ) = s1 ( p ) x0 + s2 ( p ) x0′ + s3 ( p ) x l + s4 ( p ) xl′ ,y ( p ) = s1 ( p ) y 0 + s2 ( p ) y 0′ + s3 ( p ) y l + s4 ( p ) y l′ .43x, yy0y = y( p)xl′y 0′y' = dy/dp имеют смысл косинусов угловy l′ наклона осевой линии к осям координат иylx = x( p)Ox0Отметим, что теперь наклоны x' = dx/dp иxlx0′поэтому являются компонентами каса-plтельного вектора к осевой линии.Рис. 2.2.
Параметризация балкиВведя векторы перемещений и касательных в концевых точках балки⎧x ⎫⎧ x′ ⎫⎧x ⎫⎧ x′ ⎫r0 = ⎨ 0 ⎬ = e1 , r0′ = ⎨ 0 ⎬ = e 2 , rl = ⎨ l ⎬ = e 3 , rl′ = ⎨ l ⎬ = e 4 ,⎩ y0 ⎭⎩ y 0′ ⎭⎩ yl ⎭⎩ y l′ ⎭получим параметризацию осевой линии, Рис. 2.3:⎧r0 ⎫⎪r ′ ⎪ 4⎪ ⎪r ( p ) = [s1I s2 I s3I s4 I] ⎨ 0 ⎬ = ∑ sk e k ,⎪rl ⎪ k =1⎪⎩rl′ ⎪⎭(2.1)здесь I = diag(1,1).
Знаки суммирования часто будем опускать.r0′yr( p)p=lp=0rlr0Orl′xРис. 2.3. Полученный элементЗаметим, что теперь компоненты векторов узловых перемещений и касательных векторов не обязаны быть малыми. Касательные векторы, крометого, могут иметь неединичную длину, которая отражает продольную деформацию балки. В работе [17] было показано, что балочный элемент, построенный на поле перемещений (2.1), может представлять произвольныедеформации, а также произвольные движения балки как твёрдого тела.442.2. ДЕТАЛИЗАЦИЯ УРАВНЕНИЙ ДЛЯ БАЛОЧНОГО ЭЛЕМЕНТАВ п.
1.2.3.2 уравнения движения были получены в видеM &e& + Q e = Q g .(2.2)Здесь мы приведём их явный вид, используя определение радиусвектора (2.1). Так, матрица масс будет иметь видсимм.⎤⎡ 156 I⎢⎥4l 2 Iµ l ⎢ 22l I⎥,M=13l I156 I420 ⎢ 54 I⎥⎢ − 13l I − 3l 2 I − 22l I 4l 2 I⎥⎣⎦(2.3)а столбцы обобщённых сил –⎧Q1e ⎫⎪ e⎪⎪Q ⎪eQ = ⎨ 2e ⎬ ,⎪Q 3 ⎪⎪⎩Q e4 ⎪⎭⎧µ g l 2⎫⎪ µ g l 2 12 ⎪⎪⎪Qg = ⎨⎬.gµl2⎪⎪2⎪⎩− µ g l 12 ⎪⎭Блоки Mij матрицы масс и вектора обобщённых сил тяжести Qig вычисляются по формуламllδW∂ 2Tg= µ g ∫ si dp , i, j = 1…4.M ij == µ ∫ si s j dp I = M ij I , Q i =δe∂e& i ∂e& Tji00Элементы вектора обобщённых упругих сил Q ie = ∂П ∂e i являются наиболее громоздкими, как сказано в п.
1.2.3.3, из-за сложности выражения дляпотенциальной энергии деформации (1.40). Тем не менее, в работах [17, 37,60], а также [69] разработаны различные модели обобщённых упругих сил.2.2.1.Модели обобщённых продольных силВектор обобщённых сил, вызванных продольной деформацией балки,является вектором-градиентом соответствующей потенциальной энергии Пε:∂П ε∂εdp ,Qi == EA ∫ ε∂e i∂ei0εl(2.4)куда входит продольная деформация осевой линииε = r ′T r ′ − 1 ≈ 12 (r ′T r ′ − 1) , где r' = dr/dp.(2.5)45Эта деформация и производные от неё с учётом зависимости (2.1) вычисляется следующим образом:2.2.1.1.⎞1⎛ 4 4ε ≈ ⎜⎜ ∑ ∑ sm′ sn′ e Tm e n − 1⎟⎟ ,2 ⎝ m =1 n =1⎠(2.6)4∂ε= ∑ si′sk′ e k .∂e i k =1(2.7)Модель L3Прямая подстановка выражений (2.6) и (2.7) в (2.4) приводит к обобщённым силам, которые являются кубическими функциями координат e:Qεi =4⎛ 4 4 1111 T111 ⎞⎟⎜EASeeSe=−∑ 2 ⎜ ∑ ∑ ikmn m n ik ⎟ k ∑ K ikε e k ,k =1=1 n =1k =1⎝ m4144442444443⎠4(2.8)εK ikll001111Выражения для Sikmn= ∫ si′sk′ sm′ sn′ dp и Sik11 = ∫ si′sk′ dp см.
в приложении 6.7.Выражение (2.8) может быть представлено в матричной форме:ε⎡ K 11сим.⎤ ⎧ e1 ⎫⎢ ε⎥⎪ ⎪εKKεε22⎥ ⎪⎨e 2 ⎪⎬ ,Q = K e = ⎢ ε21K εij = K ijε I ,εε⎢ K 31 K 32 K 33⎥ ⎪e 3 ⎪⎢ εεεε ⎥⎣⎢K 41 K 42 K 43 K 44 ⎦⎥ ⎪⎩e 4 ⎪⎭где Kε – матрица жёсткости.Эта модель полностью идентична модели L2, разработанной Берзéри иШабáной [17], с той лишь разницей, что у этих авторов элементы матрицыжёсткости получены непосредственно интегрированием матричных зависи1111в нашем случаемостей. Использование многоиндексных символов Sikmnможет показаться слишком трудоёмким, но можно показать, что число операций в формуле (2.8) является в точности таким же, как и в работе указанных авторов. Кроме того, последний подход позволяет вычислить не только1111.обобщённые силы, но и их матрицы Якоби, используя те же символы SikmnДля решения различных задач, например, поиск положений равновесиябалки (п. 2.3.1), интегрирование жёстких уравнений движения неявнымиметодами [39] или анализ собственных частот колебаний, необходимо вы-46числять матрицу Якоби от обобщённых упругих сил (2.8):εCij =2.2.1.2.∂Qεi∂e Tj4= K ij I + EA ∑ε41111e m e Tn .∑ Sijmnm =1 n =1Модель L2Вычисление обобщённых сил может быть упрощено, если предположить, что продольная деформация ε является постоянной по длине элемента:l⎞ 1 ⎛ 1 4 4 11 T⎞1 ⎛⎜ 1T⎟ ≈ ⎜ ∑ ∑ S mn e m e n − 1⎟ .′′ε =rrdp1−⎟⎟ 2⎜l2 ⎜⎝ l ∫0⎝ m =1 n =1⎠⎠Использование этого значения осреднённой деформации и его градиента(по-прежнему вычисляемого по формуле (2.7)) ведёт к следующей модели:4Qi = ∑εk =1EAε Sik11 e k14243εK ik4с матрицей Якоби Cεij = K ijε I + EA∑ Sik11e kk =14= ∑ K ikε e k(2.9)k =11 4 11 T∑ S jm e m .l m =1Эта промежуточная модель сил не имеет аналогов в литературе.