Using Different FEM Formulations in Calculations of Thin-Walled Structures

Full Text

1. Введение К настоящему времени тонкостенные конструкции, состоящие из тонких оболочек или их фрагментов, получают все более широкое распространение. К таковым относятся арочные конструкции, трубопроводы различного назначения, ангары, газгольдеры, резервуары, бункеры и многие другие. Прочностные расчеты подобного рода оболочечных конструкций в настоящее время базируются на основе использования численных методов, преобладающим из которых является метод конечных элементов (МКЭ) [1-5]. Во многих современных вычислительных комплексах МКЭ реализуется в форме метода перемещений [6-15]. Однако применение данного вида формулировки МКЭ к расчету оболочечных конструкций требует при использовании теории тонких оболочек [16] определения производных компонент вектора перемещения вплоть до второго порядка включительно, что сопряжено с дополнительными вычислительными трудностями. К ним относится необходимость использования интерполяционных полиномов высокого порядка с дополнительным включением в число искомых неизвестных производных нормального перемещения второго порядка [17]. Для вычисления искомых силовых неизвестных в виде продольных сил и изгибающих моментов необходимо предварительно получить деформации и искривления срединной поверхности как функций от перемещений и их производных. Альтернативой данному подходу может служить использование МКЭ в смешанной формулировке [18-26], при которой открывается возможность одновременно получать как кинематические (перемещения и их производные), так и силовые (продольные силы, моменты, деформации) искомые прочностные параметры. Естественным образом возникает вопрос о том, какая же формулировка МКЭ является наиболее предпочтительной. В данной работе на примере расчета фрагмента эллиптического цилиндра представлен сравнительный анализ эффективности использования МКЭ в форме метода перемещений и смешанного варианта МКЭ, базирующегося на применении четырехузлового элемента дискретизации с варьируемыми параметрами в виде компонент вектора перемещения и их производных первого порядка. 2. Методы 2.1. Геометрические соотношения Срединная поверхность тонкостенной конструкции эллипсоидального типа может быть задана радиус-вектором Rr0=xir+y x t j z x t k( , )r+ ( , ) ,r (1) где t - параметр эллипса поперечного сечения оболочки плоскостью, перпендикулярной оси Ox . Ковариантные векторы базиса точки M0 срединной поверхности оболочки в недеформированном состоянии определяются формулами arρ0 = Rr,0ρ (ρ=1,2 ,) (2) а орт нормали в точке M0 - векторным произведением a a a ar0 = r10×r20 0 , (3) где a0 =(ar r10 ⋅a10)(ar r20 ⋅a20)-(ar r10 ⋅a20)2. При приложении к оболочечной конструкции внешней поверхностной нагрузки точкаM0 займет новое положение M , определяемое радиус-вектором Rr r= R0 +vr = Rr0 +(v aρrρ0 +var0 ), (4) а точка M0ζ , находящаяся в произвольном слое оболочки, переместится в точку Mζ , определяемую радиус-вектором Rrζ = Rr0ζ + =Vr (Rr0 +ζar0 )+Vr. (5) Входящий в (5) вектор перемещения точки M0ζ может быть определен с привлечением гипотезы о прямой нормали [16] Vr =vr +ζ(ar r-a0 ), (6) где a a ar r r= 1× 2 a - орт нормали в точке M деформированной срединной поверхности; a Rrρ = r,ρ =(R vr0 +r) ; ar =(a ar r1 1⋅ )(a ar r2 2⋅ )-(a ar r1 2⋅ )2. ,ρ r r Дифференцированием R0ζ и Rζ по x и t определяются векторы базиса исходного и деформированного состояний grρ0 = Rr,0ρζ; grρ = Rr,ζρ, (7) скалярные произведения которых, в свою очередь, определяют компоненты метрического тензора в данных состояниях: gργ0 = gr rρ0 ⋅ gγ0; gργ = gr rρ ⋅ gγ. (8) Деформации произвольного слоя оболочки, отстоящего на расстоянии ζ от ее срединной поверхности, могут быть получены посредством использования соотношений механики сплошных сред [27]. εργζ =0.5(gργ - gργ0 ). (9) 2.2. Четырехузловой элемент дискретизации Тонкостенная конструкция моделируется четырехугольными фрагментами срединной поверхности с узлами i j k l, , , , расположенными в их вершинах. Узловыми искомыми неизвестными выбираются компоненты вектора перемещения и их производные первого порядка при реализации МКЭ в форме метода перемещений. При использовании смешанной формулировки МКЭ в качестве искомых узловых неизвестных наряду с перемещениями и их производными дополнительно привлекаются силовые параметры в виде деформаций и искривлений срединной поверхности тонкостенной конструкции. Таким образом, для используемого элемента дискретизации вводятся два столбца узловых неизвестных {UL} = {v } {v2L} { }vL ; 1 36× 112× 112× 112× (10) T 1L T T T {εℵ =T { }εργ Ty {ℵργ}Ty , (11) y 1 24× 1 12× 1 12× T где qL} ={q q q q qi j k l ,iξ...q q,lξ ,iη...q,lη}; под q понимается компонента вектора перемещения vρ 1 12× или v ; -1≤ξ, η≤1 - локальные координаты; { }εργ Ty = { } {ε11 Ty ε22} {Ty 2ε12}Ty ; 112× 1 4× 1 4× 1 4× ℵργy = ℵ { 11} {y ℵ22} {y ℵ12}y . T T T T 2 112× 1 4× 1 4× 1 4× Для интерполяции компонент вектора перемещения привлекаются произведения полиномов Эрмита третьего порядка q={ }φ T{qL}, (12) 112× 12 1× а для деформаций и искривлений - билинейные функции локальных координат -1≤ξ, η≤1, например: ε11 ={ }ψ εT { 11}y; ℵ =11 { }ψ T {ℵ11}y . (13) 1 4× 4 1× 1 4× 4 1× С учетом (12), (13) можно скомпоновать следующие матричные интерполяционные соотношения: {U}= A U{ L}= A P U R { G}; 3 1× 3 36× 36 1× 3 36× 36 36× 36 1× (14) {εργ}= H {εℵ}y, 6 1× 6 24× 24 1× (15) где { }U T ={v v v1 2 } - компоненты вектора перемещения точки M0 срединной поверхности; 1 3× { UG} - столбец узловых значений компонент вектора перемещения и их производных в глобаль- 36 1× ной системе координат x , t ; PR - матрица, определяющая связь между столбцами {UL} и {UG}. Функционал Лагранжа, необходимый для компоновки матрицы жесткости используемого элемента дискретизации, может быть записан в виде [17] ФL = {σργ} { }εργζ dV - { } { }U P dF 2V 2F или с учетом (9), (14) (16) ФL = 1{UG}T PR T B ST T С T S B dV P U R { 2 G}- (17) 1 36× 36 36× V36 6 6 3× × 3 3× 3 6 6 36× × 36 36× 36 1× 1 T 1 T -1{UG}T PR A P dFT{ } , 2 1 36× 36 36× F36 3× 3 1× где С - матрица упругости; S - матрица перехода от деформаций произвольного слоя к де- 3 3× 3 6× формациям и искривлениям срединной поверхности; B - матрица, компонуемая на основе соот- 6 36× ношений Коши и (14). Минимизацией (16) по {U G} можно получить матрицу жесткости KL и столбец узловых T усилий { f G} элемента дискретизации в формулировке метода перемещений L L 0, (18) 36 36 36 1 36 1 где KL = PR T B ST T С T S B dV P R ; {f G}= PR A T{ }P dF. 36 36× 36 36× V36 6 6 3× × 3 3× 3 6 6 36× × 36 36× 36 1× [1]× F36 3 3 1× × Смешанный функционал для получения матрицы жесткости используемого элемента дискретизации может быть записан в виде ФS = { } {εργζ T σργ}dV - 1 { } {εργζ T σργ}dV - 1 { } { }U T P dF. (19) V 1 3× 3 1× 2V 1 3× 3 1× 2F 1 3× 3 1× Учитывая вышеизложенное, функционал (19) можно преобразовать следующим образом: ФS ={ }εℵ T H T S T С S B dV P R {U G}- y 1 24× V 24 6× 6 3× 3 3× 3 6 6 36× × 36 36× 36 1× -12{ }1 24ε×ℵ T V 24 6H S× T 6 3× T 3 3С× 3 6 6 24S H dV× × {1 24ε×ℵ -}y 12{U1 36× G}T 36 36P×R T F 36 3A P dF× T { }3 1× . (20) y {εℵ}T и по { G}T , можно запиПоследовательно применяя к (20) процедуру минимизации по U y сать следующую систему матричных уравнений: ∂ФS∂ ℵ ≡{ε }Ty G U{ G}- D {εℵ =}y 0; 24 36× 36 1× 24 24× 24 1× ∂ФS∂{UG}T ≡ G T {εℵ -}y {f G}=0, (21) 36 24× 24 1× 36 1× где G = H T S T С S B dV P R ; D = H T S T С S H dV . 24 36× V 24 6× 6 3× 3 3× 3 6 6 36× × 36 36× 24 24× V 24 6× 6 3× 3 3× × ×3 6 6 24 Реализуя подстановку {εℵ =}y D -1 G U{ G} , систему (21) можно привести к виду [ ] [ ]G T D -1 [ ]G {UG}={f G} (22) 36 24 24 24 24 36× × × 36 1× 36 1× или [KS]{UG}={f G}, (23) Процедура компоновки глобальной матрицы жесткости оболочечной конструкции осуществляется стандартным для МКЭ способом [1] из KL или KS посредством матрицы индексов. 3. Результаты и обсуждения Пример расчета 1. Был рассчитан фрагмент цилиндрической оболочки в виде полукольца, загруженного вдоль образующих равномерно распределенной нагрузкой интенсивностью q=1 H/см, шарнирно закрепленного в точках B , B′ (рис.). Были использованы следующие исходные данные: E = 210⋅ 5МПа; ν =0,3; L =1см; толщина стенки h=0,2м. Расчетная схема оболочки И с т о ч н и к: выполнено М.Ю. Клочковым Shell model S o u r c e: made by M.Yu. Klochkov Первоначально был рассчитан фрагмент кругового цилиндра, параметры эллипса поперечного сечения при этом принимались равными a = b = 50 см. Расчеты были выполнены в двух вариантах: в первом варианте использовался МКЭ в форме метода перемещений (16) и (17); во втором варианте был реализован МКЭ в смешанной формулировке (19)…(23). Результаты расчетов сведены в табл. 1, в которой приведены «физические» значения кольцевых напряжений на внутренней σintt , наружной σttout и срединной σttmidl поверхностях оболочки, а также продольной силы N22 в точках A и C оболочки в зависимости от сетки узлов дискретизации. В правой крайней колонке приведены значения аналитического решения. Анализ данных, представленных в табл. 1, показывает, что и в первом и во втором варианте наблюдается устойчивая сходимость вычислений к точному аналитическому решению. В то же время следует отметить, что точность конечно-элементных решений во втором варианте оказалась выше, чем в первом варианте. Пример расчета 2. Во втором примере расчета круговой цилиндр был заменен на эллиптический с соотношением параметров эллипса поперечного сечения a b= 50 см/40 см, равного 1,25. Прочие исходные данные имели те же значения. Результаты повариантных расчетов приведены в табл. 2, структура которой совпадает с табл. 1. Таблица 1 / Table 1 Значения нормальных напряжений и продольной силы в точках приложения нагрузки в круговом цилиндре Values of normal stress and axial force at the points of load application in the circular cylinder Координаты точек x , см; t , рад / Point coordinates x , cm; t , rad σ , Н/см2; N22 , Н / σ , N/cm2; N22 , N Варианты формулировок МКЭ / FEM formulation option Аналитическое решение / Analytical solution Метод перемещений / Displacement method Смешанная формулировка / Mixed formulation Сетка узлов дискретизации / Nodal grid 101×2 151×2 201×2 101×2 151×2 201×2 Точка / Point A: x= 0,00 t =- σintt 6,190 5,470 5,250 5,600 5,250 5,130 - σttout 4,490 4,720 4,830 4,400 4,750 4,870 - σttmidl 5,340 5,100 5,040 5,000 5,000 5,000 5,000 N22 1,0672 1,0191 1,0079 0,9998 0,9999 1,0000 1,0000 Точка / Point C : x= 0,00 t = σintt 6,290 5,500 5,250 5,600 5,250 5,130 - σttout 4,400 4,690 4,830 4,400 4,750 4,870 - σttmidl 5,340 5,100 5,040 5,000 5,000 5,000 5,000 N22 1,0681 1,0192 1,0079 0,9998 0,9999 1,0000 1,0000 И с т о ч н и к: выполнено М.Ю. Клочковым / S o u r c e: made by M.Yu. Klochkov Таблица 2 / Table 2 Значения нормальных напряжений и продольной силы в точках приложения нагрузки в эллиптическом цилиндре при a b=1,25 / Values of normal stresses and axial force at the points of load application in the elliptical cylinder at ab=1.25 Координаты точек x , см; t , рад / Point coordinates x , cm; t , rad σ , Н/см2; N22 , Н / σ , N/cm2; N22 , N Метод перемещений / Displacement method Сетка узлов дискретизации / Nodal grid Аналитическое решение / Analytical solution 101×2 201×2 301×2 351×2 401×2 Точка / Point A: x= 0,00 t =- σintt 6,013 5,578 5,333 5,246 5,133 - σttout -2,489 3,837 3,854 4,345 4,869 - σttmidl 1,735 4,702 4,453 4,793 5,000 5,000 N22 0,347 0,940 0,8905 0,9586 1,000 1,000 Точка / Point C : x= 0,00 t = σintt 6,015 5,492 5,233 5,181 5,196 - σttout -2,484 3,152 5,246 4,841 4,540 - σttmidl 1,739 4,315 5,239 5,010 4,866 5,000 N22 0,349 0,863 1,049 1,002 0,9733 1,000 Окончание табл. 2 / Ending of the Table 2 Координаты точек x , см; t , рад / Point coordinates x , cm; t , rad σ , Н/см2; N22 , Н / σ , N/cm2; N22 , N Смешанная формулировка / Mixed formulation Сетка узлов дискретизации / Nodal grid Аналитическое решение / Analytical solution 101×2 201×2 251×2 251×2 Точка / Point A: x= 0,00 t =- σintt 5,588 5,123 5,123 5,037 - σttout 4,416 4,877 4,877 4,963 - σttmidl 4,999 5,000 5,000 5,000 5,000 N22 0,9997 1,000 1,000 1,000 1,000 Точка / Point C : x= 0,00 t = σintt 5,588 5,123 5,123 5,037 - σttout 4,416 4,877 4,877 4,963 - σttmidl 4,999 5,000 5,000 5,000 5,000 N22 0,9997 1,000 1,000 1,000 1,000 И с т о ч н и к: выполнено М.Ю. Клочковым / S o u r c e: made by M.Yu. Klochkov Анализ табличных данных показывает, что в первом варианте наблюдается медленная сходимость вычислительного процесса. Относительно приемлемые результаты достигаются при сетке дискретизации 351… 401 2× . Кроме того, следует отметить некоторое различие в значениях σtt и N22 в точках A и C оболочки, хотя они должны быть абсолютно одинаковыми, исходя из симметрии расчетной схемы. Во втором варианте расчета можно констатировать быструю сходимость вычислительного процесса и абсолютное совпадение σtt и N22 в точках A и C, что и должно выполняться. Таблица 3 / Table 3 Значения нормальных напряжений и продольной силы в точках приложения нагрузки в эллиптическом цилиндре при a b=2,5 Values of normal stresses and axial force at the points of load application in the elliptical cylinder at ab=2.5 Координаты точек x , см; t , рад / Point coordinates x , cm; t , rad σ , Н/см2; N22 , Н / σ , N/cm2; N22 , N Метод перемещений / Displacement method Аналитическое решение / Analytical solution Сетка узлов дискретизации / Nodal grid 101×2 201×2 301×2 401×2 501×2 Точка / Point A: x= 0,00 t = - σintt -342,0 -22,29 4,098 7,924 9,031 - σttout -880,85 -142,1 -46,52 -17,385 -39,42 - σttmidl -618,3 -83,72 -21,85 -5,048 -15,80 5,000 N22 -123,7 -16,74 -4,369 -1,010 -3,160 1,0000 Точка / Point C: x= 0,00 t = σintt -342,0 -22,23 3,830 7,507 7,062 - σttout -880,84 -141,9 -47,91 -23,50 19,99 - σttmidl -618,3 -83,60 -22,69 -8,386 13,67 5,000 N22 -123,7 -16,72 -4,538 -1,677 2,733 1,0000 Окончание табл. 3 / Ending of the Table 3 Координаты точек x , см; t , рад / Point coordinates x , cm; t , rad σ , Н/см2; N22 , Н σ , N/cm2; N22 , N Смешанная формулировка / Mixed formulation Аналитическое решение / Analytical solution Сетка узлов дискретизации / Nodal grid 101×2 151×2 201×2 251×2 Точка / Point A: x=0,00 t = - σintt 5,500 5,150 5,030 4,973 - σttout 4,512 4,850 4,970 5,025 - σttmidl 4,994 5,000 5,000 5,000 5,000 N22 0,99875 0,9996 0,9999 1,0000 1,0000 Точка / Point C: x=0,00 t = σintt 5,501 5,150 5,030 4,973 - σttout 4,512 4,850 4,970 5,025 - σttmidl 4,994 5,000 5,000 5,000 5,000 N22 0,99875 0,9996 0,9999 1,0000 1,0000 И с т о ч н и к: выполнено М.Ю. Клочковым / S o u r c e: made by M.Yu. Klochkov Пример расчета 3. В третьем примере была рассчитана оболочка с большей кривизной срединной поверхности, чем во втором примере. Соотношение параметров эллипса поперечного сечения было принято равным a b= 50 см/20 см=2,5. Результаты конечно-элементных решений при a b= 2,5 приведены в табл. 3. Как следует из табл. 3, результаты расчетов принципиально различаются по вариантам. В первом варианте численные значения напряжений и продольной силы нельзя признать удовлетворительными, несмотря на значительное сгущение сетки дискретизации. Во втором варианте можно наблюдать устойчивую сходимость численных значений к точному аналитическому решению при сравнительно редкой сетке дискретизации. 4. Заключение Таким образом можно отметить, что при выполнении прочностных расчетов тонкостенных конструкций из оболочек эллипсоидального профиля наиболее эффективным инструментом исследования НДС последних представляется МКЭ в смешанной формулировке. Предложенный вариант функционала (16) позволил получить в результате его минимизации по перемещениям и их производным, а также по деформациям и искривлениям срединной поверхности матрицу жесткости четырехузлового конечного элемента размерностью 36×36, необходимого для дискретизации рассчитываемой оболочечной конструкции. Выполненные исследования позволили сделать следующие выводы: 1. Применение метода подстановки для решения системы матричных уравнений (21) дало возможность избежать увеличения размерности матрицы жесткости конечного элемента при реализации предложенного смешанного варианта МКЭ. 2. При определении НДС тонкостенных конструкций из оболочек эллипсоидального профиля точность вычисления контролируемых параметров прочности предложенным смешанным вариантом МКЭ оказалась существенно выше по сравнению с вариантом расчета, в котором был реализован МКЭ в форме метода перемещений. 3. При расчете оболочки эллипсоидального профиля со значительной кривизной срединной поверхности корректные значения параметров НДС удалось получить только при использовании предложенного варианта смешанного МКЭ. Традиционный вариант МКЭ в форме метода перемещений в этом случае не позволил получить удовлетворительного по точности решения.

About the authors

Mikhail Yu. Klochkov

Volgograd State Technical University

Author for correspondence.
Email: m.klo4koff@yandex.ru
ORCID iD: 0000-0001-6751-4629
SPIN-code: 2767-3955

post-graduate student of the Department of Building Structures, Foundations and Reliability of Structures, Faculty of Construction and Housing and Communal Services

Volgograd, Russia

References

Supplementary files