Эффективные методы численного моделирования динамики нелинейных систем абсолютно твёрдых и деформируемых тел (1105385), страница 9
Текст из файла (страница 9)
Узловые векторы конечного элемента пластины.Далее для элементов вектора e используются обозначения e11, e12, e13, e14,…, …, e41, e42, e43, e44. Представленный набор двумерных функций формыэлемента пластины является декартовым произведением одномерных балочных функций.
Тот факт, что выражения разделены относительно p1, p2, имеетопределённое преимущество для нашего исследования: это позволяет свестидвойное интегрирование по поверхности пластины к однократным интегралам при вычислении обобщённых сил. Однако есть и некоторые неудобствапри использовании этих функций формы из-за присутствия вторых производных, которые увеличивают число степеней свободы элемента с 36 до 48 ине имеют ясного геометрического смысла. Ниже в п.
2.6.2 показана возможность исключения этих переменных из набора узловых координат (2.16).2.4.2.Матрица масс элемента пластиныВыражение для матрицы масс пластины имеет в блочной форме вид⎡ M11⎢MM = ⎢ 21⎢M 31⎢M⎣ 41M12M 22M 32M 42M13M 23M 33M 43M14 ⎤M 24 ⎥⎥,M 34 ⎥M 44 ⎥⎦где каждый блок Mij – также блочная матрица99Соглашения об использовании многоиндексных обозначений описаны в приложении 6.1.55⎡ M ij11⎢Mij 21M ij = ⎢⎢M ij 31⎢M⎣ ij 41M ij12M ij 22M ij 32M ij 42M ij13M ij 23M ij 33M ij 43M ijkl = M ijkl I ,M ij14 ⎤M ij 24 ⎥⎥,M ij 34 ⎥M ij 44 ⎥⎦(2.18)M ijkl = ∫∫ µ S ik S jl dP = µ ∫∫ sˆi sˆˆk sˆ j sˆˆl dPPPabˆ= µ ∫ sˆi sˆ j dp1 ∫ sˆˆk sˆˆl dp2 = µ Sˆij00 Sˆ kl00 .00Элементы матриц с крышками пропорциональны элементам матрицымасс для балки (2.3).Матрица S также легко вычисляется:S = S11I, S12 I, S13 I, S14 I;...;...; S 41I, S 42 I, S 43 I, S 44 I ,abˆS ij = ∫∫ S ij dP = ∫∫ sˆi sˆˆ j dP = ∫ sˆi dp1 ∫ sˆˆ j dp2 = Sˆi0 Sˆ 0j ,[]PP0(2.19)0причём символы с крышками по-прежнему пропорциональны балочным.Таким образом, получены почти все компоненты уравнений движенияэлемента пластины, за исключением обобщённых упругих сил Qe, которыенаиболее трудоёмки из-за сложности выражения для потенциальной энергии.2.4.3.Энергия деформации пластиныСледуя теории пластин Кирхгофа, энергию деформации ортотропнойпластины мы разделим на энергию деформаций пластины в срединной поверхности и энергию изгиба [9, 62]:П = Пε + Пκ ,⎛ 2 2⎞11⎜ ∑ ∑ Dijε ij2 + 2 D22⎟εε11 22 ⎟ dP ,∫∫ ⎜P ⎝ i =1 j =1⎠⎞1 ⎜⎛ 2 2κ11П = ∫∫ ∑ ∑ Dijκ ij2 + 2 D22κ 11κ 22 ⎟ dP .⎟2 P ⎜⎝ i =1 j =1⎠6П = 2hε(2.20)(2.21)Эти формулы содержат, во-первых, параметры упругости материала: цилиндрические жёсткости D11, D22 и жёсткость на кручение D12:E11h 3E22 h 3E12 h 3D11 =, D22 =, D12 = D21 =,12 (1 − ν 12ν 21 )12 (1 − ν 12ν 21 )6(2.22)5611= 0,5 ( D11ν 21 + D22ν 12 ) .
Эти жёсткостиа также коэффициент жёсткости D22зависят от модулей Юнга на растяжение E11, E22 и на сдвиг E12, а также откоэффициентов Пуассона ν12 и ν21, причём E11ν21 = E22ν12.Во-вторых, энергия деформации содержит геометрические характеристики: продольные деформации ε11, ε22, сдвиговую деформацию ε12 = ε21, атакже кривизны κ11, κ22 и кручение срединной поверхности κ12 = κ21.Используя соотношения из дифференциальной геометрии поверхностей[10], вычислим для нашей параметризованной срединной поверхностиr = r(p1, p2) деформацииε ij = 12 (riT r j − δ ij )(2.23)с символами Кронекера δij, а также кривизныκ ij = rijT n n3(2.24)с вектором нормали n = r1 × r2 . Другие использованные обозначения:∂r∂ 2riijri == Smne mn , rij == S mne mn .∂pi∂pi ∂p j(2.25)(Подразумевается суммирование по m и n).Вычислив градиенты потенциальных энергий (2.20) и (2.21), получимвекторы обобщённых продольных и поперечных сил. Эти громоздкие вычисления в различных вариантах были проведены в работе [73], см.
п. 2.4.4.Покажем, что уравнение (2.23) выражает нелинейные соотношения Грина между перемещениями и деформациями [8, 9]. Действительно, представим, что деформированное состояние пластины задано уравнениями⎧ p1 ⎫ ⎧u1 ( p1 , p2 ) ⎫⎪ ⎪ ⎪⎪r = u∗ + u = ⎨ p2 ⎬ + ⎨u2 ( p1 , p2 )⎬ .⎪ 0 ⎪ ⎪u ( p , p ) ⎪⎩ ⎭ ⎩ 3 1 2 ⎭Здесь u* определяет начальную (плоскую) форму пластины, а u соответствует отклонениям от неё. Если мы применим зависимости (2.23), то получим компоненты тензора деформации поверхности, содержащие как линейные, так и квадратичные члены, что и требовалось показать:571 ⎛ ∂u3∂u∂u ∂u ⎞ε ij = ⎜⎜ j + i + ∑ k k ⎟⎟ ,2 ⎝ ∂pi ∂p j k =1 ∂pi ∂p j ⎠i, j = 1…2.Точно также можно убедиться [72, 73], что выражение (2.24) даётправильное значение компонентов тензора кривизны поверхности, линейныечлены которых совпадают с используемыми в теории малых деформаций.2.4.4.Модели обобщённых сил от деформаций всрединной поверхности пластиныИскомые обобщённые силы – это градиенты энергии деформации (2.20):∂U ε 12Qkl ==∂e kl h 2 ∫∫Pε∂ε⎛∂ε 2211 ∂ε11⎜ Dijε ij ij + D22εε+2211⎜∂e kl∂e kl∂e kl⎝⎞⎟ dP .⎟⎠(2.26)Продольные и сдвиговые деформации (2.23) и их градиенты с учётомсоотношения (2.14) принимают вид()1 T i⊗ je mn S mnpq e pq − δ ij ,(2.27)2∂ε iji⊗ j= S klrse rs ,(2.28)∂e kl1 i ji⊗ jjгде введены символы S mnpq= S mnS pq + S mnS ipq .2Несколько моделей обобщённых сил вида (2.26), имеющие различнуюε ij =()сложность, были получены с учётом различных допущений.
Приведём однуиз них, модель L3, как наиболее зарекомендовавшую себя в расчётах.В модели L3 не используются никакие упрощения для вычисления компонент деформации срединной поверхности (2.27). После интегрированиявыражения для обобщённых сил (2.26) будут кубическими по координатам e:Qεkl = K εklmn e mn ,(2.29)с использованием следующей матрицы жёсткости:()pqrs T′K εklmn = K klmne pq e rs − K klmnI,гдеmnpqK klrs=6h2(D Si⊗ j; i⊗ jij klrs;mnpq′ =K klrs6h2())1⊗1; 2⊗22⊗2; 1⊗111+ D22S klrs;mnpq + S klrs;mnpq ,()i⊗ j11.δ ij Dij + D22S klrs(2.30)58i ⊗ j; i ⊗ ji⊗ ji⊗ ji⊗ j i⊗ jЗдесь символы S klrs= ∫∫ S klrsdP и S klrs;mnpq = ∫∫ S klrs S mnpq dP имеютPPзначения, приведенные в приложении 6.8.Как и в случае с балочными элементами (п.
2.2.1), для эффективногорешения жёстких уравнений движения необходимо вычислять матрицы Якоби от обобщённых сил. Это производится достаточно просто:Cεklmn =2.4.5.∂Qεkl∂e Tmnmnpq= K εklmn + 2 K klrse rs e Τpq .Модели обобщённых сил от поперечных деформацийДанная часть обобщённых сил является градиентом энергии (2.21):Qκkl =∂U κ=∂e kl ∫∫P∂κ⎛∂κ 22 ∂κ1111⎜ Dijκ ij ij + D22κκ 22+11⎜∂∂eee∂klklkl⎝⎞⎟ dP.⎟⎠(2.31)Вычислим кривúзны (2.24) срединной поверхности и их производные10:κ ij = rijT n f 3 ,∂κ ij∂e klf = n = nT n ,TT⎞ 3⎡ ∂n ⎤∂f1 ⎛⎜ ⎡ ∂rij ⎤,= 3 ⎢ T ⎥ n + ⎢ T ⎥ rij ⎟ − 4 (rijT n)⎟ f∂ee∂f ⎜ ⎣ ∂e kl ⎦kl⎣ kl ⎦⎝⎠(2.32)T1 ⎡ ∂n ⎤∂f= ⎢ T ⎥ n.f ⎣ ∂e kl ⎦∂e klРазличные модели поперечных сил были разработаны в работе [73].Приведём здесь описание одной из них, модели T2. В ней реализована идеяусреднить вектор нормали n = r1 × r2 по поверхности пластины, чтобы упростить подынтегральное выражение в (2.31).
Примем, что вектор среднейнормали к пластине равен среднему арифметическому нормалей в её углах11:n=14(e 21 × e12 + e 41 × e 32 + e 23 × e14 + e 43 × e 34 ) ,f = n = nT n .Тогда производные в формулых (2.32) принимают явный вид10Здесь используются операции дифференцирования вектора по строке координат, их матрицыЯкоби, а также другие абстракции, описанные в приложении 6.1.11См. Рис. 2.9 и примечание после формулы (2.16).59T⎡ ∂n ⎤⎢ T ⎥ rij⎣ ∂e kl ⎦⎧+ 14 e k −1, l +1⎪b kl = ⎨− 14 e k +1, l −1⎪0⎩∂κ ij∂f1= b kl × n ,∂e klf= b kl × rij ,for {k , l} ∈ {{2,1},{4,1}, {2,3}, {4,3}},for {k , l} ∈ {{1,2},{3,2}, {1,4}, {3,4}},for the rest combinations of k , l ,()()13 ij TijijSn+Sb×e+S e n n × b kl .mnklmnkl35 mn mn∂e klffПодставляя вычисленные кривизны и их градиенты в выражение (2.31),и, наконец,=получим обобщённые силы, которые квадратичны по e:Qκkl =⎛ ∗⎞1 T3∗∗(e pq n ) ⎜⎜ Sklpqn + Smnpqb kl × e mn + 2 (e Tmn n ) Smnpqn × b kl ⎟⎟ ,3ff⎝⎠()∗ijij1111222211,S mnpq= Dij S mnpq+ D22S mnpq+ S mnpqijijнаходятся в приложении 6.8.явные выражения для символов S mnpqЭти силы могут быть записаны в виде Qκkl = K κklpq e pq , как в случае продольных сил (2.29).
Кроме того, расчёты показывают12, что эти выражениямогут быть значительно упрощены удерживанием лишь первого слагаемого:Qκkl ≈1 T∗(e pq n ) Sklpqn.3fТакже при моделировании используется следующее упрощённое выражение для матрицы Якоби этих сил:κCklmn∂Qκkl1 ∗= T ≈ 3 Sklmnn nT .∂e mn f2.5. ПРИМЕРЫ МОДЕЛИРОВАНИЯ МЕМБРАН И ПЛАСТИНВ этом пункте приведены результаты решения тестовых задач, в которых проверялась адекватность полученных моделей пластинчатых элементов.
Сравнение с результатами натурных экспериментов проведено в главе 3.Численные эксперименты выполнялись на программном комплексе УМ13, гдереализованы описанные конечные элементы пластины.12Результаты вычислений с использованием полного и сокращённого выражений с высокой точностью совпадают.13Ссылку на программу см. в п. 2.3.602.5.1.Статические деформации тяжёлой мембраныЦелью данного теста была проверка как возможностей функций форм,так и сходимости результатов. Первая модель представляла собой тяжёлуюэластичную квадратную мембрану с параметрами: размеры a×b×h = 1×1×0,01м, плотность материала µ = 1000 кг/м3, модуль Юнга E = 105 Па, коэффициент Пуассона ν = 0,3. Мембрана подвешена за три своих угла, так что четвёртый свободно свисает – точка E на Рис.
2.10. Для вычисления обобщённыхупругих сил использовалась модель L3.На рисунке приведены положения мембраны в положении равновесияпри использовании различного числа конечных элементов.Рис. 2.10. Тяжёлая мембрана, подвешенная за три угла:12, 22, 32, 42, 62 и 82 конечных элементовКак видно из рисунка, форма положения равновесия сходится плохо.Это объясняется тем, что её предельноеКоордината точки E , м1состояние соответствует излому по-минус zE0.9верхности мембраны из-за отсутствия0.80.7изгибной жёсткости, а гладкие функ-0.6ции формы не могут описать этого.xE и yE0.5Иначе говоря, данная задача плохообусловлена.
Тем не менее, значения0.412345678Число элементов на сторону пластиныкоординат точки E сходятся хорошо,Рис. 2.11.Рис. 2.11. Сходимость координат точки E61На Рис. 2.12 мембрана размерами 2×1×0,01 м подвешена за все свои четыре угла. Эта хорошо обусловленная задача показывает хорошую сходимость как формы равнвесия, так и координат точки E (zE → –0,58 м).Рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
