Конечно-элементный Анализ устойчивости упругопластической сферической оболочки при всестороннем сжатии Текст научной статьи по специальности «Физика»

Механика

Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 3 (1), с. 158-1(52

УДК 539.3

КОНЕЧНО-ЭЛЕМЕНТНЫЙ АНАЛИЗ УСТОЙЧИВОСТИ УПРУГОПЛАСТИЧЕСКОЙ СФЕРИЧЕСКОЙ ОБОЛОЧКИ ПРИ ВСЕСТОРОННЕМ СЖАТИИ

© 2011 г. А.А. Артемьева х, М.С. Баранова *, А.И. Кибец \ В.И. Романов 2,

А.А. Рябов 2, Д.В. Шошин 1

1 НИИ механики Нижегородского госуниверситета им. Н.И. Лобачевского

2 РФЯЦ-ВНИИЭФ, г. Саров Нижегородской обл.

aranan@mail. rn

Поступила в ридакцию 25.01.2011

В осесимметричной постановке рассматривается потеря устойчивости и закритическое поведение сферической оболочки при всестороннем сжатии. Решение геометрически и физически нелинейной задачи основано на методе конечных элементов и явной конечно-разностной схеме интегрирования по времени типа «крест». Результаты расчета сопоставляются с экспериментальными данными.

Ключивые слова: сферическая оболочка, устойчивость, геометрическая и физическая нелинейность, метод конечных элементов.

Обзор работ, посвященных математическому моделированию потери устойчивости сферических куполов при сжатии, приведен в [1—5]. Сопоставление теоретических значений критических нагрузок с экспериментальными данными показало, что между ними возможно существенное (в 3-4 раза и более) различие. Это связано как с особенностями применяемых математических моделей, так и с условиями проведения экспериментов. В первую очередь сюда относятся начальные отклонения от правильной формы оболочки, несовершенство условий ее закрепления и нагружения. Цель настоящей работы - конечно-элементное решение геометрически и физически нелинейной задачи потери устойчивости и закритического поведения сферической оболочки при ее равномерном всестороннем сжатии, сопоставление результатов расчета с экспериментальными данными и оценка точности применяемой вычислительной модели.

Деформирование оболочки рассматривается в осесимметричной постановке и описывается в переменных Лагранжа с позиций механики сплошных сред в цилиндрической системе координат Х1, X 2 (Х1 - ось вращения, Х2 - радиальная ось). Уравнение движения выводится из баланса виртуальных мощностей [6]:

Ц (°1158П + °225 8 22 + °33 5 8 33 + 2°125812 ) Х

□

хХ2У+ Ц(рЦбЦ + ри2Ьй2)Х2d□ - (1)

□

-1 (Р15й1 + р25й2) X2^ = 0,

где 8.. и Ст.., У, . = 1,3- компоненты тензоров

У У ~ ~

скоростей деформаций и напряжений, й{ - перемещения в системе координат Х, V - параметр симметрии, равный нулю, если рассматривается плоское деформирование, и равный единице, если рассматривается осесимметричное деформирование; р - плотность, р{ - распределенная

нагрузка, □ - область, занимаемая конструкцией, Ор - область действия внешнего давления, точка над символом означает частную производную по времени t. По повторяющимся индексам ведется суммирование. Скорости деформаций определяются в метрике текущего состояния:

¿11 U1,1, ¿22 U2,2 , ¿33 VX2 U2,

12 = 1 (f>1,2 + и 2,1), ®= 1 (£>1,2 - £> 2,1),

t

J UJdt.

(2)

и. . = 8£J./8X., X. = X.\, ° +

‘,J 4 J ’ J J It _°

Уравнения состояния устанавливаются раз-

• V „V •!

дельно для шаровых 8 , О и девиаторных 8..,

У

О. составляющих скоростей деформаций и напряжений. Зависимость шаровых компонент скоростей деформаций и напряжений имеет вид

ог = -3а: 8г ,

¿V _ (¿11 + 822 +^33)

3

где К - коэффициент объемного сжатия.

(3)

p

При упругопластическом деформировании девиаторные составляющие скорости деформации 8. раскладываются на пластические 8'р и

У У

упругие 8. компоненты:

У

8' = 8'е +8'р 8' р +8' р +8' р = 0

ьи у у ’ ьп ^ь22 ь33 -°-

Девиаторные составляющие тензора напряжений вычисляются с помощью соотношений теории течения с кинематическим и изотропным упрочнением:

D а.. = 2Gs'e, S.. = а., — р., р. = 2ев ■

її і 11^11 11 >і і <Ь 11

р.. = f p.dt, =IS.., S ..S.. < -а2, (4)

П/ J 'У ’ . У’ У У 3 T ’

Л

ОТ = ОТ <Х 128 X 128 = Л2 8'р8'р , X = | 128^

’ 3 0

Здесь Dj - производная по Яуманну, G -модуль сдвига, от - предел текучести, Sjj и р.

- компоненты тензоров активных напряжений и микронапряжений. Система уравнений (1)-(4) дополняется начальными условиями и кинематическими граничными условиями.

Решение определяющей системы уравнений (1)-(4) основано на методе конечных элементов и явной конечно-разностной схеме интегрирования по времени типа «крест» [6-8]. Расчетная область покрывается лагранжевой сеткой из 4-узловых конечных элементов. В узлах сетки

определяются перемещения й, скорости й и ускорения й в общей системе координат X = [Х1 X2]Т , используемой для стыковки конечных элементов (КЭ). В каждом элементе вводится локальный прямоугольный базис |х} = [х1 х2 ]Т, отслеживающий его вращение

как жесткого целого пошаговым пересчетом направляющих косинусов осей х.

Конечный элемент, в общем случае искаженный, с помощью билинейного изопарамет-рического преобразования отображается на квадрат -1 < ^ < 1 (У = 1,2):

X =Х N. (^ |2) X ,(У = 1,2),

к=1

Nl^l, ^) = (1 -^(1 Ч2)/4 =

= (1 -^-|2 + Щ/4,

N 2^!, |2) = (1 + ^)(1 —1-)/4 = = (1 + Si — ^ — Щ/4, N3(^1,12) = (1 + ^)(1 + іу/4 = = (1 + |, +^2 +^1^2>/4, N 4(^1,^2> = (1 — У(1 + ^2>/4 = = (1 — ^1 + ^2 — ^1^2) / 4

где хк - координаты узлов в базисе j x}. Компоненты скорости узлов КЭ проектируются в локальный базис и аппроксимируются внутри элемента с помощью функций формы Nk из (5):

и, =1 ukN &, ^). (6)

к=1

Группируя коэффициенты при одинаковых степенях ^1^2, получим

x i = a3 + а&1 + a^2 + ^1^2 , (7)

ui = d3 + d1^1 + d1 ^2 + d1^1^2. (8)

Скорости деформаций e аппроксимируются в виде суммы безмоментных и моментных составляющих:

811 = 811 + 811^2 , 822 = 822 + S22^1 ,

812 = 812 ’ 833 = 833 + 833^1 + 833^2 , (9)

где s3k = const - значения компонент скорости

деформаций в центре КЭ, 81кк = 5s kkj52, х =

=const - значения компонент градиента скорости деформаций в центре КЭ. Полные деформации аппроксимируются линейными функциями аналогично (9). Пластические и упругие компоненты деформаций определяются из уравнений состояния и могут иметь в пределах конечного элемента более сложную зависимость от пространственных координат. В соответствии с этим для уточнения численного решения значения пластических деформаций и напряжений определяются в выбранных фиксированных точках конечного элемента. Для улучшения сходимости численного решения пластические деформации в элементе корректируются из условия а и, = 3, i = 1,2 [9]. При выполнении

численного интегрирования в (1) для определения статически эквивалентных узловых сил применяются квадратурные формулы [13]. В каждом конечном элементе мощность виртуальной работы в (1) выражается через матрицу масс, узловые ускорения, силы и виртуальные скорости перемещения [6]. Заменяя в (1) интегрирование по области Q суммированием по элементам, получим дискретный аналог уравнений движения

(10)

(5)

где ті - массы узлов КЭ-сетки; и/ , F¡J - компоненты ускорений узлов и результирующих узловых сил в общей системе координат X. Для интегрирования определяющей системы уравнений по времени воспользуемся явной конечноразностной схемой типа «крест». Величина временных шагов выбирается исходя из условия устойчивости Куранта.

АI = 0

х. N N. "V. _ ДГ = ( ✓ _ «и* (.5 мс ✓ 9

Ч _ * ✓ у«.ф

ц Д? = 1 2 мс ***

** ’’ ф*

Д?=1 б мс

Рис. 1

Рис. 2

Изложенный выше 4-узловой конечный элемент реализован в рамках вычислительной системы «Динамика-3» [11]. На его основе решена осесимметричная задача упругопластического деформирования медной сферы (модуль упругости Е = 68.6 ГПа, коэффициент Пуассона ц = 0.376, плотность р = 8.95 г/см3, предел текучести от = 0.108 ГПа, модуль линейного упрочнения g = 0.582 ГПа) постоянной толщины ^1уЪ0 = 23.4, где R0 - начальный радиус, ^ -начальная толщина) при линейном увеличении сжимающего давления, равномерно распределенного по внешней поверхности. Скорость нарастания внешнего давления выбиралась из условия квазистатического нагружения, исключающего динамические эффекты. Для численного исследования выпучивания и закритического поведения сферической оболочки применялась сетка п х п2 = 1 х 40, где п1, п2 - число конечных элементов по толщине и меридиану оболочки

соответственно. Для инициирования потери устойчивости оболочки в подобласти 0° < ф < 18° (ф - угол, отсчитанный от оси вращения оболочки) задавался ее начальный прогиб, линейно распределенный по углу поворота с максимальным значением w0 = 0.01Л0 в полюсе.

Результаты численного решения задачи представлены на рис. 1, 2. На рис. 1 изображены формы верхней половины образующей оболочки для различных временных слоев, отсчитываемых от начала потери устойчивости р(г0) = ркр, Дг = г - г0. На рис. 2 сравниваются

остаточные конфигурации оболочки, полученные в расчете (пунктирная линия) и в эксперименте [12] (сплошная линия).

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

женном состоянии. Изменение окружных и меридиональных напряжений по толщине не превышает 8%. Выпучивание сферы начинается при р = ррр~ 0.0092 ГПа. На начальном этапе закритического поведения оболочки вмятина носит локальный характер (ф < 4°), что позволяет применить теорию пологих оболочек для определения критического значения внешнего давления р;р. В [1] приведена формула, полученная из аналитического решения упругой задачи в геометрически нелинейной постановке: рар = 0.2Е(Н / Д)2/(1 -ц2)1/2. (11)

Для рассматриваемой оболочки р^р = =0.027 ГПа, что превышает в 2.9 раза ракр , найденное из численного решения задачи. Но, как отмечено в [1], формула (11) не учитывает начальное несовершенство формы оболочки, которое снижает значение критической нагрузки. Кроме того, согласно результатам конечноэлементного расчета, к моменту потери устойчивости оболочки интенсивность напряжений достигает предела текучести, и ее дальнейшее деформирование происходит уже в упругопластической области, а формула (11) выведена в предположении упругого деформирования материала. Следует отметить, что критическое напряжение сткр, соответствующее потере устойчивости сферической оболочки, связано с критическим давлением рр следующей формулой [1]:

ст = р 0.5Л /Ъ. (12)

кр * кр 4 '

Подставляя в (12) р^, найденное из конечно-элементного решения рассматриваемой задачи, получим сткр и стТ, что подтверждается

данными численного решения.

Изгиб оболочки после потери устойчивости сопровождается возрастанием скорости смещения ее полюса. Вокруг него образуется локальная вмятина, радиус и глубина которой со временем увеличиваются. На начальном этапе (Дг <0.8 мс) узлы конечно-элементной сетки оболочки смещаются внутрь исходной сферы: для любого узла оболочки ее текущий радиус г (ф) < Д0. Когда глубина вмятины в полюсе

достигает приблизительно % радиуса оболочки R0 (Дг = 1.2 мс), под действием меридиональных напряжений ее сегмент за точкой перегиба выходит за пределы начальной конфигурации (см. рис. 1: г(ф) > Я0 при 60° < ф < 90°), что соответствует экспериментальным данным [12]. Анализ напряженно-деформированного состояния выявил значительное искривление траекто-

рии деформаций в пространстве Ильюшина [13] в области перегиба образующей. Это позволяет утверждать, что изгиб сферического купола происходит в условиях сложного нагружения. Процесс упругопластического выпучивания оболочки до перехода в новую равновесную форму продолжается примерно 2 мс. При этом полюсная точка разгоняется до максимальной скорости V и 97 м/с. Из рис. 2 видно, что в процессе деформирования произошло значительное изменение формы оболочки. Расчетное wр и

экспериментальное '№э значения перемещения

полюса (W = w / Ъ) в остаточном положении равны 28.5 и 26.8 соответственно. Их расхождение 5 = ^ - wэ) х 100% / Wр составляет 6%.

Хорошее качественное и количественное согласование расчетных и экспериментальных данных по остаточной форме указывает на то, что вычислительная система «Динамика-3» и реализованная в нем конечно-элементная методика позволяют с достаточной степенью точности описывать потерю устойчивости и закритиче-ское поведение сферической оболочки в условиях сложного нагружения и больших формоизменений.

Работа выполнена в рамках реализации ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы, гранта Президента РФ для поддержки ведущих научных школ РФ НШ- 4807.2010.8, а также при поддержке РФФИ (проекты №11-08-00557-а, 11-08-97023-

р_поволжье_а ).

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

1. Вольмир А.С. Устойчивость деформируемых систем. М.: Наука, 1967. 984 с.

2. Коноплев Ю.Г., Саченков А.В. Исследование прочности и устойчивости пологих сферических оболочек под действием локальных нагрузок // Ис-след. по теор. пластин и оболочек. Казань: Изд-во Казанского ун-та, 1967. Вып. 5. С. 161-188.

3. Ворович И.И., Минакова Н.И. Проблема устойчивости и численные методы в теории сферических оболочек // Итоги науки и техники. Механика твердых деформируемых тел. М.: ВИНИТИ, 1973. Т.

7. С. 5-86.

4. Григолюк Э.И., Мамай В.И. Механика деформирования сферических оболочек. М.: Изд-во МГУ, 1983. 114 с.

5. Гудрамович В.С. Устойчивость упругопластических оболочек. Киев: Наукова думка, 1987. 216 с.

6. Баженов В.Г., Прокопенко М.Б. Численное решение осесимметричных нелиненых нестационарных задач динамики составных упругопластических конструкций // Прикладные проблемы прочности и пластичности. Численное моделирование физикомеханических процессов: Межвуз. сб. Нижегородский ун-т, 1991. С. 55-63.

7. Баженов В.Г., Кибец А.И., Цветкова И.Н. Численное моделирование нестационарных процессов ударного взаимодействия деформируемых элементов конструкций // Проблемы машиностроения и надежности машин. 1995. № 2. С. 20-26.

8. Артемьева А.А., Баженов В.Г., Кибец А.И. и др. Верификация конечно-элементного решения трехмерных нестационарных задач упругопластического деформирования, устойчивости и закритиче-ского поведения оболочек // Вычислительная механика сплошных сред. 2010. Т. 3. № 2. С. 5-14.

9. Метод конечных элементов в механике твердых тел / Под ред. А.С. Сахарова и И. Альтенбаха. Киев: Вища школа - Лейпциг: ФЕБ Фахбухферлаг, 1982. 480 с.

10. Демидович Б.Д., Марон И.А. Основы вычислительной математики. М.: Наука, 1970.

11. Программный продукт «Пакет прикладных программ для решения трехмерных задач нестационарного деформирования конструкций, включающих массивные тела и оболочки, «Динамика-3» (ШШ «Динамика 3»). Сертификат соответствия Госстандарта России № РОСС RU.ME20.H00338/2000.

12. Рябов А.А., Романов В.И., Зефиров С.В. Численное исследование упругопластического выпучивания сферической оболочки // Прикладные проблемы прочности и пластичности. Численное моделирование физико-механических процессов: Межвуз. сб. М.: Товарищество научных изданий КМК, 1999. С. 125-128.

13. Ильюшин А.А. Пластичность. Основы общей математической теории. М.: Изд-во АН СССР, 1963.

FINITE-ELEMENT ANALYSIS OF ELASTOPLASTIC SPHERICAL SHELL STABILITY UNDER ALL-AROUND COMPRESSION

A.A. Artemyeva, M.S. Baranova, A.I. Kibets, V.I. Romanov,

A.A. Ryabov, D. V. Shoshin

In the axisymmetric formulation, the buckling failure and postbuckling behavior of an elastoplastic spherical shell under uniform all-around compression are considered. The solution of the geometrically and physically nonlinear problem is based on the finite element method and the explicit finite-difference time-domain scheme. The computation results are compared with experimental data.

Keywords: spherical shell, stability, geometric and physical nonlinearity, finite element method.