CC BY
7
Численные методы расчета конструкций
РАЗРАБОТКА И ОЦЕНКА ЧИСЛЕННОГО АЛГОРИТМА РАСЧЕТА НА УСТОЙЧИВОСТЬ НЕЛИНЕЙНО ДЕФОРМИРУЕМОЙ ПОЛОГОЙ ЦИЛИНДРИЧЕСКОЙ ОБОЛОЧКИ
С.А. ИВАНОВ, младший научный сотрудник ЦНИИСК им. В. А. Кучеренко,
109428Москва, ул. 2-я Институтская, д. 6; ivanov@tsniisk.ru
В статье исследовано поведение пологой цилиндрической оболочки из упругопла-стического материала при нагружении равномерно распределенной вертикальной нагрузки. Приведены основы методики решения задачи с применением вариационно- разностного подхода и теории малых упругопластических деформаций. Получены кривые равновесных состояний оболочки, кривая зависимости верхней критической нагрузки от кривизны оболочки в безразмерном виде, приведено сравнение величин нагрузок с результатами расчета по методу конечных элементов.
КЛЮЧЕВЫЕ СЛОВА: пологая цилиндрическая оболочка, устойчивость, геометрическая нелинейность, физическая нелинейность, упругопластический материал, вариационно-разностный метод, кривая равновесных состояний, критическая нагрузка.
I. Введение
В данной статье представлены результаты численного исследования поведения пологой цилиндрической оболочки на прямоугольном плане из материала, обладающего нелинейной зависимостью между напряжениями и деформациями, при действии равномерно распределенной нагрузки. Для решения задачи применен метод шагового нагружения в сочетании с вариационно- разностным подходом (ВРМ). Кинематическая модель построена с учетом деформаций поперечного сдвига и может быть использована для расчета как тонких оболочек, так и оболочек средней толщины. Разработанный алгоритм реализован на ЭВМ в виде программы на языке ФОРТРАН-90, программа протестирована на задаче об изгибе пластинки из упруго-пластического материала.
II. Постановка задачи
Задача расчета нелинейно деформируемой оболочки формулируется как вариационная и сводится к нахождению функций перемещений, дающих минимальное значение полной потенциальной энергии системы. Выражение для полной потенциальной энергии оболочки, занимающей область □, имеет вид:
^ п п
где & = (епе22еикпк22киеие23) - вектор, компонентами которого являются составляющие тензора деформаций; N = (АГП И22 И12МпМ22М12 QlЪ Q2Ъ )т - вектор, компонентами которого являются усилия, определяемые по соотношениям теории малых упруго-пластических деформаций; с] = (£/, т1 т2 с/2 )Т - вектор внешней нагрузки; V = (и V вх 02 м>)т - вектор, компонентами которого являются функции перемещений; и,у,в1,в2,м>- соответственно функции тангенциальных перемещений, углов поворота поперечных сечений и прогиба.
Компоненты тензора деформаций для пологой оболочки записываются с учетом деформаций поперечного сдвига [1]:
ди м> 1 (дм/ \ дv 1
еи =--1-—ь— — ; е77=-+ —
11 дх Я 2{дх) ду 2
2 л. ,
<3и>
¿¿V ди д\\! дм)
ду) дх ду дх ду'
дви
дв.
2 .
дв, дв, дм! „ ^ ^ ¿ту дх ду
(2)
¿7Х
где Л - радиус кривизны срединной поверхности оболочки вдоль направляющей. Усилия N11, Ы22, ..., 623 , действующие в оболочке и входящие в выражение (1), определяются по известным формулам [2], согласно теории малых упруго-пластических деформаций:
^22 = В 21 (а1' а2 )е11+^22 (а1 ■> а2 )£22+^21 (а1 ■> а2 )К\\Г^22 (а1' а2 )К22 з Ы12 = Въъ(а1,а2)£12+С зз(а,1,с);2)лг12;
Ми = Сп(а1,а2)£п+С12(а1,а2)Е22+Вп(а1,а2)кп+В12(а1,а2)к:22;
а2)х:11+£)22(а1,ог2)Аг22; (3)
Ми = С33(а1,а2)е12+033(а1,а2)к12; 013 = А,3(а15а2)£13; б23 = В33(аиа2)£23,
ы? ст 1
Д11(а1,а2) = Д22(а1,а2)= — ---¿¿г;
-а/2
А/2
, (7 1/
2?12(а1;а2) = 2?21(а1;а2) = —--с1г;
-Ы2£, 1-у
где
А/2 а 2
Сп(а1,а2) = С22(а1,а2)= / —-
-к 12 I
А/2 А/2
-V СГ
йг;
(• С • И
С12(а15а2) = С21(а15а2)= ] ——
-А/2 1 ~
'12
А/2
-V2
А/2
-¿¿г;
(4)
ст 1 сг г
В33(а1,а2)= | ——-С33(а15а2)= | — —--
-А/2 ^ 2(1 + У) -А/2 2(1 + 1/)
йг;
ст- г2
А/2<т,
Д1(а1,а2)= [ —--¿¡г; В^а^а^ = В21(аг,а2) = [ —
-А/2^1 "V2 -А/2^1"^2
йг ;
Въъ(сс1,ос2)= |
А/2 (7,. г2
-¿¿г.
-А/2^- 2(1 + у)
В формулах (4) о, - интенсивность напряжений; е, - интенсивность деформаций. При активной деформации интенсивность напряжений о{ является вполне определенной функцией интенсивности деформаций е{ независимо от вида напряженного состояния, то есть о, = ог(ег). Коэффициент Пуассона V изменяется в процессе деформирования по закону:
у = (5)
2 2 егЬ0
где Е0 и у0 - соответственно модуль упругости и коэффициент Пуассона в начальном состоянии. В рамках принятой модели интенсивность деформаций определяется по формуле:
1
е< =■
1-у'
1-У + У
>121 + ^2 4у-у2 -1^11^22 +
3(1 -У)2 /2 2 , 2
^12 + £13 + 23 4
4
1/2
Теория малых упругопластических деформаций, строго говоря, справедлива лишь при простом нагружении, однако и в случае сложного нагружения, достаточно близкого к простому, указанная теория дает вполне приемлемые по точности результаты близкие к экспериментальным данным. Все формулиров-
ки, приводимые для деформационной теории пластичности, справедливы и в случае нелинейно-упругого материала.
III. Методика решения
Задача решалась вариационно-разностным методом в сочетании с шаговым нагружением. Условия минимума функционала (1) имеют вид
dIJ/dzi =F(z,p) = 0 , (6)
где z - вектор узловых перемещений, p - параметр нагрузки. В качестве ведущего параметра принимается длина дуги кривой равновесных состояний As. Тогда дополнительно к основной системе n уравнений (6) вводится n +1-е уравнение вида F0(z,p,s) = 0.
Основная система и вспомогательное уравнение решаются двумя различными методами: Ньютона-Рафсона и самокорректирующимся методом Рунге-Кутта. При итерационном подходе на каждом шаге m, которому соответствует значение параметра продолжения sm, искомые функции z(sm) и p(sm) находятся путем последовательных приближений к точному решению с использованием итерационных формул и вспомогательного уравнения:
R(zk(sm))Azk+1 =QA рк+1 + Qpk(sm)-VW(zk(sm));
z*+i =zt+rtAzt+1, (7) + (8)
где k - номер итерации; m - номер шага по ведущему параметру; zk - вектор узловых перемещений на итерации с номером k; zk+1 - значение искомой функции на итерации с номером k + 1; Azk+1 - вектор приращения узловых перемещений на итерации с номером k + 1; R(zk (sm)) - матрица Гессе на итерации k; VIV(/,-. (sm)) - градиент функции W(/,-. (sm)), представляющей собой потенциальную энергию деформации; Q - нормированный вектор внешней нагрузки; 5zк+1 = 5zk+Azk+1; Ърк+1 =Ърк +Арк+1.
Формулы (7) и (8) описывают итерационный процесс нахождения решения на сфере с центром в точке С ^т ^Р ^т ^ и радиусом As- (схема Крисфилда).
Компоненты вектора градиента и элементы матрицы Гессе вычисляются по формулам на основе явных выражений:
dz, JJ ^ '
d2W(z) „ ¿e) тъ*^ - 11 -+8 Df
>) (9)
, - ¿п.
дг, дг, дг,
3 1 j J
В формулах (9) присутствует упругопластическая матрица -О(е), связывающая усилия и деформации на основе соотношений (3). Элементы матрицы В(е) зависят от интенсивности деформаций и вычисляются на основе реальной диаграммы деформирования материала.
В качестве тестовой задачи исследована квадратная пластинка со сторонами длиной 1 м и толщиной 1 см, свободно опертая по контуру. Для этой задачи разработанный алгоритм дает хорошую сходимость. Значения нагрузки и перемещений близки к представленным в работе [3], где они были получены методом конечных разностей, однако без учета деформаций поперечного сдвига.
Исследована цилиндрическая оболочка на квадратном плане (рис.1) с размерами в плане а = Ь = 1 м, толщиной к = 1 см. Граничные условия соответствуют шарнирно-подвижному опиранию по всем сторонам (опирание оболочки на вертикальные диафрагмы), т. е. приняты следующие соотношения:
при х = 0 и х — а и = 91=>у = 0; при у = 0 и у = Ь V = 02 = = 0.
Оболочка нагружена равномерно распределенной вертикальной нагрузкой, нагружение только активное. В качестве зависимости о,- (е,) принята диаграмма деформирования низкоуглеродистой стали при одноосном растяжении-сжатии, с пределом текучести от = 2100 кг/см2, начальным модулем упругости Е0 = 2Д-106 кг/см2. Диаграмма, аппроксимированная естественными кубическими сплайнами, приведена на рис. 2. Начальное значение коэффициента Пуассона Vo = 0,3, затем на каждой итерации он пересчитывается по формуле (5). Оболочка покрывается сеткой 8^8 элементов, в силу симметрии задачи расчет производится для четверти оболочки.
Рис. 2. Диаграмма работы материала (зависимость о, (е,)).
IV. Результаты расчета
На рис. 3 и 4 приведены для сравнения кривые равновесных состояний оболочки, полученные различными методами, при толщине И = 1 см и различных значениях главной кривизны к1 = 1/Л. Кривые без значков на рисунке соответствуют к\ = 0,2 м-1, отмеченные квадратами - к1 = 0,3 м-1, треугольниками - к1 = 0,4 м-1, кругами - к1 = 0,5 м-1. Для каждого значения кривизны получены две кривых: одна по расчету с учетом геометрической и физической нелинейности (сплошная на рис. 3, 4), другая с учетом только геометрической нелинейности (штриховая на рис. 3, 4). Кривые на рис. 3 получены с использованием самокорректирующегося метода Рунге-Кутта, на рис. 4 - метода Ньютона-Рафсона. При заданных величинах параметров вычислительного процесса для оболочки с кривизной к1 = 0,5 м-1 значение верхней критической нагрузки, полученное методом Рунге-Кутта, выше значения, полученного методом Ньютона-Рафсона.
Данное различие может являться следствием накопления погрешностей при шаговом вычислительном процессе, реализующем процедуру Рунге-Кутта. На рис. 5 приведены кривые зависимости безразмерной верхней критической
нагрузки д*сг =(дсг/Е0)(а/И)4от безразмерной кривизны оболочки к[ — кх а2/И.
Рис. 3. Кривые равновесных состояний оболочки. Самокорректирующийся метод Рунге-Кутта
0 0.5 1 1.5 2
Рис. 4. Кривые равновесных состояний оболочки. Метод Ньютона-Рафсона
Обе кривые получены с применением метода Ньютона-Рафсона, в качестве ведущего параметра принята длина дуги кривой равновесных состояний по схеме Крисфилда, с шагами приращения ведущего параметра Ду= 0,10 и Ду = 0,15. Штриховой линией показана кривая, полученная при расчетах в физически линейной постановке (линейно-упругий материал с Е0 = 2,1-10б кг/см2, v0 = 0,3), сплошной - при расчетах в физически нелинейной постановке (упругопласти-ческий материал с диаграммой деформирования, приведенной на рис. 2).
10 20 30 40 50 60
Рис. 5. Зависимость безразмерной верхней критической нагрузки от безразмерной
кривизны оболочки
В табл. 1 приведены результаты расчета цилиндрической оболочки с к1 = 0,4 м-1 и к = 1 см без учета нелинейной работы материала по разработанной методике и их сопоставление с результатами, полученными при расчете методом конечных элементов в ПК Лира 9.4.
Таблица 1
Величина ВРМ, метод Рунге-Кутта, сетка 8x8 ВРМ, метод Ньютона-Рафсона, сетка 8x8 МКЭ, сетка 8x8
Нагрузка, соответствующая первой предельной точке, кг/см2 4,154 4,011 3,9
Прогиб центра оболочки, соответствующий первой предельной точке, см 1,3 1,4 0,8078
В табл. 2 приведены результаты расчета цилиндрической оболочки с к1 = 0,4 м-1 и к = 1 см с учетом нелинейной работы материала по разработанной методике и их сопоставление с результатами, полученными при расчете методом конечных элементов в ПК Лира 9.4. Таблица 2
Величина ВРМ, метод Рунге-Кутта, сетка 8x8 ВРМ, метод Ньютона-Рафсона, сетка 8x8 МКЭ, сетка 8x8
Нагрузка, соответствующая первой предельной точке, кг/см2 3,121 2,986 3,2
Прогиб центра оболочки, соответствующий первой предельной точке, см 0,800 0,800 0,8069
V. Заключение
Вариационно-разностный алгоритм решения задачи устойчивости пологой цилиндрической оболочки позволяет получить картину поведения оболочки из нелинейно-упругого или упругопластического материала и оценить устойчивость оболочки по величине верхней критической нагрузки. Результаты весьма близки к значениям, полученным при расчете методом конечных элементов при той же густоте сетки, как для линейно-упругого материала, так и для упруго-пластического. Значения безразмерной верхней критической нагрузки qcr для оболочки из упругопластического материалов меньше, чем для оболочки из линейно-упругого материала при том же значении безразмерной кривизны к*, и при увеличении кривизны разность этих нагрузок возрастает.
Л и т е р а т у р а
1. Милейковский И.Е., Трушин С.И. Расчет тонкостенных конструкций. - М.: Стройиздат, 1989. - 200 с.
2. Трушин С. И. Метод конечных элементов. Теория и задачи: Учебное пособие. -М.: Издательство АСВ, 2008. - 256 с.
3. Стрельбицкая А.И., Колгадин В.А., Матошко С.И. Изгиб прямоугольных пластин за пределом упругости. - Киев: Наукова Думка, 1971. - 243 с.
WORKING OUT AND ESTIMATION OF NUMERICAL ALGORITHM OF STABILITY ANALYSIS OF THE NONLINEAR DEFORMABLE SHALLOW CYLINDRICAL SHELL
Ivanov S.A.
In this article the behavior of the shallow elastic-plastic cylindrical shell under the uniformly distributed load is analyzed. The basics of the problem solving by applying the finite difference energy method and incremental technique are presented.
The equilibrium curves and the dimensionless correlation of the upper critical load and the shell's curvature were obtained as the result. The upper critical load values were compared with the values that were obtained by the finite element method.
KEY WORDS: shallow cylindrical shell, stability, geometrical nonlinearity, physical nonlinearity, elastic-plastic material, finite difference method, equilibrium curve, critical load.
