Сравнительный анализ эффективности использования конечных элементов различной мерности при анализе НДС тонких оболочек

Полный текст

Введение Тонкостенные конструкции из пластин и оболочек являются неотъемлемой частью современных зданий, сооружений и архитектурных форм [1-3]. В настоящее время при анализе НДС конструкций из оболочек на первый план выдвинулись преимущественно численные методы расчета [4-10], в частности МКЭ [11-18]. В наши дни матрицы жесткости конечных элементов формируются на основе двух теорий: теории тонких оболочек с использованием теории прямой нормали для сведения трехмерного НДС к двумерному состоянию, а также соотношений теории упругости без дополнительных гипотез о деформировании нормального элемента. В последнее время многими исследователями, занимающимися данной проблематикой, высказывается мнение о том, что при построении конечноэлементных моделей оболочечных конструкций предпочтение следует отдавать трехмерным конечным элементам как наиболее универсальным. При этом в качестве основных геометрических соотношений рекомендуется использовать соотношения теории упругости. Данный подход является вполне обоснованным при исследовании НДС оболочечных конструкций средней толщины и толстостенных сосудов. Однако, использование трехмерных конечных элементов при анализе НДС тонких оболочек является не вполне оправданным и требует детального обоснования. В связи с этим достаточно актуальной остается задача сравнитель- ного анализа эффективности использования дву- {r0 r0 r0 }; {r0 }T {r0 r0 r0 }; и трехмерных конечных элементов при определении НДС тонких оболочек. = a1,x a2,x a,x a,θ 1´3 = a1,θ a2,θ a,θ В настоящей работе на примере расчета жесткозащемленной по торцам цилиндрической оболочки выполнен сравнительный анализ точности конечноэлементных решений, полученных при использовании двуи трехмерных элементов дискретизации. Геометрия оболочки матрица [m0] - нулевая, а в матрице [n0] ненулевыми элементами являются n22, n23 и n32. Положение точки оболочки, отстоящей от срединной поверхности на расстоянии t в исходном M 0t и деформированном M t состояниях, определяется соответствующими радиус-векторами Радиус-вектор, задающий срединную поверх- R r0t R = r0 a +tr0; t R R r = r0t r +V, (6) ность эллиптического цилиндра, может быть представлен выражением где V - вектор перемещения точки из положе- 0t t r r r r ния M в положение M . R0 = xi + r(θ)sinθj + r(θ)cosθk, (1) Дифференцированием (6) по x, θ и t могут быть получены ковариантные векторы базиса исгде x - осевая координата, θ - угол, отсчитываемый от оси Oz против хода часовой стрелки в ходного и деформированного состояний r r r r r плоскости, перпендикулярной оси Ox. 0 0t 0 0 0t 0 ( ) Входящие в (1) функция r(θ) определяется формулой g1 = R,x = a1 ; g2 = R,θ = a2 1 + tn32 ; r(x,θ) = bc , b2 cos2 θ + c2 sin2 θ (2) r 0 g3 = R0t = r 0 ; ,t a r g1 = Rt = r0 + ,x g1 r V,x ; r rt g2 = R,θ = r0 r g2 +V,θ ; r t r0 r где b и c - параметры эллипса. g3 = R,t = g3 +V,t . (7) Ковариантные векторы локального базиса точки срединной поверхности определяются дифференцированием (1) по x и θ ; r0 r0 r Входящие в (7) производные вектора перемещения точки M 0t могут быть представлены компонентами, отнесенными к базису отсчетной поверхности r0 0 a2 = R,θ = a1 = R,x = i (r,θ (θ)sinθ + r ) r (θ)cosθ j + V,x = t1 0 2 0 0 r r r 1 a1 + t1 a2 + t1 a ; 1 r 0 2 r 0 r 0 +(r,θ (θ)cosθ - r (θ)sinθ)k . (3) V,θ r = t2 a1 r · t2 a2 r · t2 a ; r 1 0 2 0 3 0 . Орт нормали к срединной поверхности опре- V,t = t3a1 · t3 a2 + t a (8) деляется векторным произведением r r r r r Ковариантные компоненты тензора деформаций могут быть получены из известного соотноa 0 = a 0 ´ a 0 / a 0 ´ a 0 . (4) шения механики сплошной среды 1 2 1 2 Произведение векторов локального базиса (3), εt =(g - g0 )/ 2. (9) (4) по x, θ могут быть представлены в матричном виде следующим образом r r r r mn mn mn Входящие в (9) ковариантные компоненты метрического тензора определяются скалярными {a 0 } = ëém0 ûù {a 0 }; {a 0 } = éën0 ùû {a 0 }, (5) произведениями (7) ,x ,θ 3´1 3´3 3´1 3´1 3´3 3´1 r r r r где g0 = g0 × g0; gmn = gm × gn. {r0}T {r0 r0 r0 }; {r0 }T mn m n (10) a = a1 a2 a a,x = Соотношения (8) в развернутом виде для трех- 1´3 1´3 мерной формулировки будут иметь следующую структуру t 1 t 1 0 2 0 2 ü ε11 = t1 ; 2ε12 = t2 + a22t1 + a22n32t × t1 ; { G}Т = ïì{1G}Т { 2G}Т { G}Т ï, 2ε13 = t3 + t1 ; Uy 1´96 í vу ïî1´32 vу 1´32 vу ý 1´32 ïþ (16) t 1 t 0 2 0 2 L Т i j k l m n p h i h i h i h ε22 = a22t2 + a22n32t × t2 ; где {q } ={q q q q q q q q q ...qξqη ...qηqς ...qς}; у ,ξ 1´32 , , , , , t 0 2 0 2 t G Т i j k l m n p h i h i h i h 2ε23 = a22t3 + a22n32t × t3 + t2 ; ε33 = t3 , (11) {qу } ={q q q q q q q q q,x ...q,xq,θ ...q,θq,t ...q,t }; r r 1´32 где a0 = a0 × a0. q - компонента вектора перемещения V точки 22 2 2 Элементы дискретизации Элементом дискретизации в двумерной постановке выбирается четырехугольный фрагмент срединной поверхности с узлами i, j, k, l, расположенными в его вершинах. Столбцы узловых не- M 0t, отстоящей от отсчетной поверхности на расстояние t. В трехмерном элементе дискретизации функции формы представляют собой триадные произведения полиномов Эрмита третьей степени: q = {ψ}Т {qL}Т . известных в локальной ξ, η и глобальной x, θ системах координат имеют следующий вид [19]: 1´32 y 32´1 (17) { L}Т ì{ 1L}Т { 2L}Т { L}Т ü; Матрицы жесткости и столбцы узловых двуи трехмерного элементов дискретизации форми- Uy = í vу vу vу ý (12) ровались стандартным образом путем минимиза- 1´36 î 1´12 1´12 1´12 þ ции функционала Лагранжа [12; 17; 21]. { G}Т ì{ 1G}Т { 2G}Т { G}Т ü, Пример расчета Uy = í vу vу vу ý (13) 1´36 î 1´12 1´12 1´12 þ Было исследовано НДС жестко защемленного L Т i j k l i l i l где {qу } 1´12 = {q q q q q,ξ ...q,ξ q,η ...q,η}; по торцам цилиндра, нагруженного внутренним давлением интенсивности q. Приняты следующие G Т i j k l i l i l исходные данные: радиус срединной поверхности {qу } 1´12 ={q q q q q,x ...q,x q,θ ...q,θ}; R = 1,0 м; длина образующей L = 1,0 м; толщина стенки h = 0,02 м; q = 5 МПа; E = 2·105 МПа; q - компонента вектора перемещения точки срединной поверхности. Элементы матрицы-строки функций формы {φ}T для двумерного элемента дискретизации представлены диадными произведениями полиномов Эрмита третьей степени v = 0,3. Расчеты выполнялись по двум вариантам: в первом в качестве элемента дискретизации использовался двумерный конечный элемент 36×36 с девятью неизвестными в узле (12), (13); во втором применялся объемный конечный элемент 96×96 с двенадцатью узловыми варьируемыми параметq = {j}Т { y q L }Т . (14) рами. В первом варианте при интерполировании 1´12 12´1 В качестве объемного элемента дискретизации выбирается восьмиузловой шестигранник с узлами i, j, k, l, m, n, p, h, расположенными в его вершинах [20]. Столбцы узловых неизвестных трехмерного элемента дискретизации в локальной ξ, η, ς и глобальной x, θ, t системах координат имеют следующую структуру: по площади элемента использовалась квадратура Гаусса с 6×6 точками интегрирования, а по толщине - формула Симпсона с семью точками интегрирования. Во втором варианте при интегрировании по объему элемента также использовалась квадратура Гаусса с 6×6×6 точками интегрирования. В силу наличия плоскостей симметрии рассчитывалась 1/8 часть цилиндра. Результаты повариантных расчетов представлены в табл. 1, в которой приведены значения в L Т ïì L Т L Т L Т ïü меридиональных напряжений на внутренней σ у {Uy } = í{v1 у } {v2 } {vу } ý; (15) и наружной σн поверхностях цилиндра в опорном 1´96 ïî1´32 1´32 1´32 ïþ и пролетном сечениях в зависимости от густоты сетки дискретизации, причем во втором варианте толщина сетки цилиндра моделировалась только одним рядом объемных конечных элементов. Анализ данных, представленных в таб. 1, показывает, что и в первом, и во втором вариантах расчета наблюдается устойчивая сходимость вычислительного процесса. Однако значения напряжений в опорном сечении во втором варианте оказались существенно меньше (примерно на 100 МПа) по сравнению с первым вариантом расчета. Численные значения напряжений при использовании двуи трехмерного конечных элементов [Table 1. Numerical stress values using twoand three-dimensional finite elements] Таблица 1 Вариант расчета [Variant of calculation] Сетка дискретизации [Sampling grid] Сечение [Section] Общее число неизвестных [Total number of unknowns] Опорное [Reference] Пролетное [Span] Напряжения, МПа [Stress, MPa] σв σн σв σн I 2×9 460,1 -339,9 68,45 50,93 78 2×17 475,5 -356,0 68,01 51,38 158 2×33 479,6 -360,2 67,90 51,48 318 2×65 480,7 -361,3 67,88 51,51 638 2×97 480,9 -361,5 67,87 51,51 958 II 2×9×2 341,6 -230,5 63,85 53,60 216 2×17×2 375,6 -263,6 64,45 53,52 440 2×33×2 381,0 -267,5 64,57 53,55 888 2×65×2 381,8 -266,7 64,61 53,57 1784 2×97×2 382,4 -266,5 64,62 53,58 2680 2×129×2 383,0 -266,7 64,62 53,58 3576 Численные значения напряжений при использовании трехмерного конечного элемента с дополнительной дискретизацией цилиндра по толщине Таблица 2 [Table 2. Numerical stress values using a three-dimensional finite element with additional thickness discretization of the cylinder] Число рядов элементов по толщине цилиндра [The number of rows of elements through the thickness of the cylinder] Сетка дискретизации [Sampling grid] Сечение [Section] Общее число неизвестных [Total number of unknowns] Опорное [Reference] Пролетное [Span] Напряжения, МПа [Stress, MPa] σв σн σв σн 2 2×9×3 370,3 -254,2 65,27 52,76 324 2×17×3 423,4 -299,5 65,27 52,76 660 2×33×3 439,4 -305,8 65,53 52,62 1332 2×65×3 445,8 -308,2 65,61 52,59 2676 2×97×3 448,0 -310,4 65,63 52,58 4020 2×129×3 449,2 -311,9 65,64 52,58 5364 3 2×9×4 378,0 -260,7 64,57 53,11 432 2×17×4 442,0 -314,3 65,43 52,63 880 2×33×4 471,8 -330,4 65,69 52,47 1776 2×65×4 487,9 -340,1 65,78 52,43 3568 2×97×4 491,7 -343,6 65,80 52,43 5360 2×129×4 492,8 -345,0 65,81 52,43 7152 В табл. 2 приведены результаты второго варианта расчета при дополнительном разбиении стенки цилиндра по толщине на 2 и 3 элемента. Анализ значений напряжений, представленных в табл. 2, показывает, что и при разбиении сетки цилиндра по толщине на 2 ряда объемных конечных элемен- тов значения напряжений в опорном сечении оказались заниженными примерно на 10 % по сравнению с первым вариантом расчета. И только при разбиении стенки цилиндра по толщине на 3 ряда объемных элементов значения напряжений в опорном сечении достигли уровня значений напряжений первого варианта расчета (при сетке дискретизации по площади 2×65). Значения напряжений в пролетном сечении оказались примерно одинаковыми в обоих вариантах расчета при любой сетке дискретизации. Заключение Учитывая, что общее число искомых неизвестных при сетке 2×65 в первом варианте расчета в 5,6 раза меньше, чем во втором варианте при аналогичной сетке дискретизации (при разбиении стенки цилиндра по толщине на 3 ряда элементов), можно сделать вывод о том, что использование двумерных элементов дискретизации при анализе НДС тонких оболочек вполне обосновано и является более целесообразным с точки зрения затрат машинного времени и ресурсов применяемой компьютерной техники.

Об авторах

Юрий Васильевич Клочков

Волгоградский государственный аграрный университет

Автор, ответственный за переписку.
Email: klotchkov@bk.ru

доктор технических наук, профессор, заведующий кафедрой высшей математики

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Анатолий Петрович Николаев

Волгоградский государственный аграрный университет

Email: anpetr40@yandex.ru

доктор технических наук, профессор, профессор кафедры прикладной геодезии, природообустройства и водопользования

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Татьяна Алексеевна Соболевская

Волгоградский государственный аграрный университет

Email: moonway13@rambler.ru

кандидат технических наук, доцент кафедры высшей математики

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Михаил Юрьевич Клочков

Московский государственный университет им. М.В. Ломоносова

Email: m.klo4koff@yandex.ru

студент третьего курса физического факультета

Российская Федерация, 119991, Москва, Ленинские горы, 1

Список литературы