CC BY
9
Численные методы расчета конструкций
КОНЕЧНО-ЭЛЕМЕНТНЫЙ АНАЛИЗ ОСЕСИММЕТРИЧНО НАГРУЖЕННЫХ ОБОЛОЧЕК ВРАЩЕНИЯ С ВЕТВЯЩИМСЯ МЕРИДИАНОМ ПРИ УПРУГО-ПЛАСТИЧЕСКОМ ДЕФОРМИРОВАНИИ
Ю.В. КЛОЧКОВ, доктор техн. наук, профессор А.П. НИКОЛАЕВ, доктор техн. наук, профессор А.Ш. ДЖАБРАИЛОВ, канд. техн. наук, доцент Волгоградский государственный аграрный университет 400002, г. Волгоград, пр. Университетский, 26.
В настоящей работе излагается алгоритм расчета осесимметрично нагруженных оболочек вращения с ветвящимся меридианом на основе теории малых упругопла-стических деформаций.
КЛЮЧЕВЫЕ СЛОВА: оболочка вращения, матрица жесткости, конечный элемент, ветвящийся меридиан.
Оболочки вращения сложной геометрии, в том числе с ветвящимся меридианом, находят широкое применение в строительной, авиакосмической, нефтяной и других отраслях современного машиностроения [1, 2, 3]. В зонах ветвления меридиана возникают концентрации напряжений, вследствие чего, именно в этих областях материал оболочки переходит в пластическую стадию работы. Поэтому при расчете подобного рода конструкций необходимо обязательно учитывать физическую нелинейность применяемого материала. [4, 5].
1. Конечный элемент и интерполяция перемещений. В качестве элемента дискретизации выбирается фрагмент меридиана с узлами / и /, выделенный двумя плоскостями, перпендикулярными оси вращения, который для удобства численного интегрирования отображается на прямоугольный отрезок с локальной координатой -1 < / < 1. Связь между глобальной £ (длина дуги меридиана) и локальной координатами осуществляется выражением:
^ = + И^/, (1)
2 2
где и - координаты узлов конечного элемента в глобальной системе координат.
При использовании векторной интерполяции перемещений [6], в сочетании с шаговой процедурой нагружения, матрицы-строки векторных узловых неизвестных в локальной и глобальной системах координат выбираются в виде:
{\Уу ^ = \аУ' АУ]АУ/ АУ/}, \АУГ } = {\У' АУ]АУ^ АУ,{ }, (2) 1х4 1х4
На основании (1) можно сформировать матричное соотношение между
к; } и ЙТ }:
Ау }=м{АУ/ }. (3)
4x1 4x4 4x1
Входящие в структуру \аУТ} приращения векторов перемещений узлов
конечного элемента и приращения их производных можно представить в локальном узловом базисе:
АУт= Аита0т + Аwma0т; АУ,? = I^ а0т + ¡^ а0т, (4)
где Аит и Аwm - приращения меридиональной и нормальной компонент векто-
1т т 1т т
ра перемещения узлов конечного элемента (КЭ); ¡,1 , ¡,1 , /,ц , /,ц- многочлены, содержащие приращения Аит, Аwm, а также их первые производные по дуге меридиана верхний индекс т обозначает узлы КЭ. На основании (4) столбец {аУу } можно представить матричным произведением:
кг }=ик } (5)
4х1 4х8 8x1
где {пу У = А1 Аwi Аи]' Аw j¡,S ¡,\ ¡1 ¡,{ }
В структуру матрицы [А] входят узловые векторы локального базиса а10т, а0т, которые могут быть представлены в локальном базисе внутренней точки КЭ:
а0т = dmlal0 + dm2a0; а0т = ¿Ца? + ¿2т2а0. (6)
Приращение вектора перемещения внутренней точки КЭ на (j +1) -м шаге нагружения интерполируются через столбец (2):
АУ = (}Т ^У;}, (7)
1x4 4x1
где {р }Т = (р2.....р4 } - матрица строка, содержащая полиномы Эрмита 3-й сте-
lx 4
пени. Входящие в (7) приращение вектора перемещения внутренней точки КЭ и его производные по глобальной координате £ запишутся в виде
АУ = Аиа0 + Аwa0; АУ, 5 = ¡1а10 + ¡1а 0; АУ, ю = ¡^50 + ¡11а 0, (8)
где ¡1 = Аи,5 -kАw ; ¡1 = Аuk + Аw,5 ¡^ = Аи,^ -k(Аuk + 2Аw,5 ) - Аwk,5 ,
¡11 = Аw,^ +k(2Аи,5-Аwk) + А^,5 .
С учетом (3), (5) интерполяционное выражение (7) запишется следующим образом:
АУ = (}Т ШМ ={р}Т Шк}. (9)
lx4 4х4 4х8 8 Х1 1x4 4x8 8 х8 8 Х1
Входящую в выражение (9) матрицу [А] с учетом (6) можно представить матричной суммой:
А] = а]0 [А1 ]+ а0 [А ]. (10)
4x8 4x8 4x8 Принимая во внимание (8) и (10), соотношение (9) запишем следующим образом:
Аиа0 + Аwa0 = {р }Т (а0 [А1 ] + а0 [А2 ])[G]{пу } (11)
Путем приравнивания соответствующих выражений при ортах а10 и а0, стоящих в левой и правой частях (11), можно записать следующие интерполяционные зависимости для приращений компонент вектора перемещения на (j +1 )-м шаге нагружения:
Аи = {у}Т А ][о]{пу }; А^ = {у}Т [Л2 %у }. (12)
Дифференцируя (7) по дуге меридиана 5, можно получить производную приращения вектора перемещения:
ау ,,=у л} ^ {у } . (13)
С учетом (4) и (5)...(10), соотношение (13) можно записать в виде:
« + /1« 0 = У, / } ^ А ] + а 0 [Л2 ])° ]{гу } . (14)
Из (14) можно получить интерполяционные зависимости следующего вида:
/,1=У „ ¡д5 [Л11° Ьу } /,1=У / № ^} (15)
Аналогично выводятся соотношения для /)ц и /,ц:
/,¡1 4„, } (§)2 [Л1 Ю]ру } /,,1 = У„,/ } (§)2 [Л2 ЮЬУ } (16)
Входящий в структуру (15), (16) столбец {пу } может быть выражен через стандартный столбец узловых варьируемых параметров:
{пу }=[г]{АиТ }, (17)
где {АиТ }= |{диуГ} {^Т } |; } = {у Ад,^ Ад,I}
Здесь под Ад понимается приращение меридиональной Аи или нормальной Ам> компонент вектора перемещения на (/ + /)-м шаге нагружения.
2. Матрица жесткости и столбец узловых нагрузок на (г + у) -м шаге нагружения. Равенство работ внешних и внутренних сил на(/ + 1)-м шаге нагру-жения можно записать в следующем виде:
т
|Цв} М+^К =|К}(М+М)^, (18)
V F
где
т 4 £ПА£22 } - приращения деформаций в произвольном слое об°
лочки; } {аст0 Т = {аст11 Аст22 } - напряжения и их приращения на (/ +1)-м шаге нагружения; {Аи }Т = {АиАм>} -приращения компонент вектора перемещения на (/ + 1)-м шаге нагружения; {Р},{АР}- столбец внешней поверхностной нагрузки и её приращение на (/ +1 )-м шаге нагружения.
Входящие в структуру функционала (18) приращения напряжений на (/ +1 )-м шаге нагружения {аст°^} могут быть выражены через приращения деформаций {А^ }:
^Мо^} (19)
Входящая в (19) матрица пластичности [СП] на шаге нагружения формируется на основе деформационной теории пластичности [7] путем дифференцирования компонент тензора напряжений по компонентам тензора деформаций [8].
Приращения деформаций в произвольном слое оболочки могут быть выражены через приращения деформаций срединной поверхности:
52
kß} =[г ]{Aa},
''aß 2 x1
(20)
2x4 4x1
где [Г ] =
2x4
1 0 4 0
0 1 0 4
Входящий в (20) столбец приращений деформаций в точке срединной поверхности на основании соотношений теории тонких оболочек [9] может быть
выражен через столбец узловых варьируемых параметров {а} :
{Аеар Цв]^ } (21)
Столбец {Аи } с учетом (12), (17) может быть представлен матричным произведением:
{Аи } (22)
{(У А ЫТ ] {р }Т [а2 ][ы 1т ][]
С учетом (19)...(22) функционал (18) запишется в следующем виде:
{< } [Т ]Т | [5]Т [г Т [сп ][г ^У [Т ]{< }=
где
[N ] =
y )
(23)
можно записать
J[N]T {AP }dF + J[N J p}dF -J[b]T [Г]Т [G]dV
v F F V y
В результате минимизации функционала (23) по {aU^} следующее матричное соотношение:
[K ][AUyr ]={f } + {Я }, (24)
где [K] = [T]T J [B]T [Г]T [Сп \Г][B]dV ; {f } = J[Nf {AP }dF - матрица жестко-
V F
сти и столбец узловых нагрузок КЭ на (j + 1)-м шаге нагружения;
( \
lTfnW ГГпиГгЛТ
- поправка Ньютона-Рафсона.
M = IW PW [г]Т Gdv
\F V ,
3. Условия сопряжения в узле ветвления меридиана. При выводе условий сопряжения n оболочек в узле ветвления меридиана столбец узловых неизвестных одной из них принимается за основной. Столбцы узловых варьируемых параметров остальных (n - 1) оболочек должны быть выражены через столбец узловых неизвестных основной оболочки, исходя из кинематических условий сопряжения [10]. Одним из таких условий является инвариантность приращений векторов перемещений на (j +1) -м шаге нагружения:
ЛУ(1) = ЛУ(2) = ... = Л% = ... = Л V(n), (25)
где AV(i) = Au(if(i) + Aw(i)a(i) ;
Vçî ) -орт i -ой оболочки, касательный к мери-
0
диану и направленный к узлу ветвления; а^ - орт нормали.
Базисные орты и а0) можно выразить через базисные векторы декартовой системы координат:
Н-) = М(г).} где = ^5(0)} {}Г = {,к} (26)
Принимая, например, первую оболочку за основную можно, с использованием (26) сформировать матричную зависимость между базисными ортами - -й и основной оболочек:
И(г) = Ь М(1), где Ьг] = И(г) И-1). (27)
Из (25) с учетом (27) можно получить выражения для приращений компонент вектора перемещения - -й оболочки через приращения компонент основной оболочки:
Аи (.) = /1(Аи (1), А^(1)); ) = /2 (Аи(1), а^(1)). (28)
При жестком соединении п оболочек в качестве условия сопряжения можно использовать предположение о том, что углы поворота нормалей в узле ветвления в процессе деформирования остаются равными для всех оболочек:
А^(1) а0 = 50 = 50 = = ^к5о (28)
д^) (1) asV(2) (2) ...... дsV(i) (') ...... дsv{n) <п) ' ^
Из (28) можно получить выражения для приращения производной компоненты вектора перемещения - -й оболочки через узловые неизвестные основной оболочки: Ам>,Sv(i) = /3 (А^,Sv(l), Ам(ц) (29)
Компоненты Аи,S^)(г = 1..п), всех сопрягаемых оболочек остаются свободно варьируемыми. В результате столбец узловых неизвестных основной оболочки в узле сопряжения примет вид:
{АиУ }(1) = S(1), S(1), Аu,S(2), Аu,ф),...,Аu,*(п)} (30)
На основании (28)...(30) формируются матричные соотношения:
К" },) = [Р ](.) № }(1). (31)
Матрицы жесткостей и столбцы узловых нагрузок элементов дискретизации (п - 1) оболочек, непосредственно примыкающих к узлу ветвления меридиана умножаются на входящую в (34) матрицу р ]:
Р Г К ^1 Р Г({/}+{*}) (32)
Пример расчета. Была решена задача по определению напряженно- деформированного состояния (НДС) оболочечной конструкции, состоящей из цилиндра и примыкающих к нему двух конусов, загруженной внутренним давлением интенсивности q (Рис 1.). Были приняты следующие исходные данные:
модуль упругости Е = 7,5 -105МПа, коэффициент Пуассона V = 0,32., Яц = 90см, t = 1см, q = 0,5МПа , 11 = 80см, 12 = 100см, 13 = 60см, а = 450, / = 30°.
В расчете использовалась прямая упрочнений (рис. 2), которая была получена с помощью формулы:
= к(£ - £0 ^^ (33)
где к = 14444,4 МПа; £0 = 0,002 ст0 = 200МПа.
Результаты расчетов представлены в таблице 1, в которой приведены численные значения меридиональных (стм) и кольцевых (стк ).напряжений на внешней и внутренней поверхностях цилиндра в опорном сечении I, в узле ветвления меридиана II, а также в сечениях конусов III и IV, отстоящих от узла
ветвления на 1^2 толщины оболочки. Каждая сопрягаемая оболочка при этом разбивалась на 50 конечных элементов.
Таблица 1
Сечение Напряжения, МПа Число шагов нагружения
90 110 130
I аМ 3,97 3,97 3,97
а"м 3,95 3,95 3,95
аВ 45,0 45,0 45,0
а 45,0 45,0 45,0
II аМ 208,801 208,801 208,801
а"м -222,060 -222,059 -222,059
аВ 35,477 35,474 35,471
а"к -135,461 -135,469 -135,474
III аМ 476,26 476,264 476,969
аМ -471,62 -471,624 -472,323
а'к 201,64 201,640 201,990
а"к -222,77 -222,775 -223,092
IV аМ 708,05 708,057 709,148
аМ -720,2 -720,203 -721,291
аВ 309,31 309,31 309,828
аВ -343,77 -343,770 -344,275
Анализ табличного материала показывает устойчивую сходимость вычислительного процесса с увеличением числа шагов нагружения. Максимальных значений напряжения достигают, как и следовало ожидать, в зонах, близких к узлу ветвления меридиана. Расчетная схема позволяет вычислить меридиональное напряжение на срединной поверхности цилиндра в опорном сечении, исходя из условия равновесия конструкции.
Давление, действующее на нижний конус, вызывает появление растягивающих напряжений (а р ), а на нижний конус - сжимающих (ас):
я(R2 - r,2) q n(r2 - R2) q
а -w х = 21,46МПа ас = —-ц— ^ = 17,5МПа .
2яИц t c 2яRц t
цц
Результирующее меридиональное напряжение на срединной поверхности цилиндра в опорном сечении равно:
ам = ар + ас = 21,46 -17,5 = 3,96МПа.
Кольцевые напряжения в опорном сечении могут быть вычислены по формуле ак = qR4jt = 0,5МПа • 90см/ 1см = 45МПа./ Из таблицы видно, что в опорном сечении I, значения ам и ок, полученные в результате конечно-элементного расчета совпадают с ом и ок, вычисленным исходя из условия равновесия.
На основании вышеизложенного можно сделать вывод о достаточной точности разработанного алгоритма конечно-элементного расчета осесимметрично нагруженной оболочки вращения с ветвящимся меридианом при упругопласти-ческом деформировании.
Л и т е р а т у р а
1. Zhang Junbo, Li Xikui (2011). A mixed finite element and mesh-free method using linear complementarity theory for gradient plasticity, Comput. Mech, 47, №2, p. 171-185.
2. Galishin A.Z., Shevchenko Yu. N.(2011). Determining the axisymmetric elastoplastic state of thin shells with allowance for the third invariant of the stress deviator, Int. Appl. Mech, 46, №8, p. 868-876.
3. Poh L. H., Peerlings R.H. J., Geers M.C.D., Swaddiwudhipong S. (2011). An implicit tensorial gradient plasticity model-formulation and comparison with a scalar gradient model, Int. Solids and Struct., 48, №8, p. 2595-2604.
4. Shang B.P., Du. H. (2010). Application and finite element simulation of plastic blanking deformation, Comput and Theor. Nanosci., 5, №8, p. 1631-1635.
5. Dzierzanowski Grzegorz (2012). Stress energy minimization as a tool in the material layout desingn of shallow shells, Int. J. Solids and Struct, 49, № 11-12, p. 1343-1354.
6. Dzhabrailov A. Sh., Klochkov Yu. V., Marchenko S.S., Nikolaev A. P. (2007).The finite element approximation of vector fields in curvilinear coordinates, Russian Aeronautics, Vol. 50, № 2, p. 115-120.
7. Малинин Н.Н. (1975). Прикладная теория пластичности и ползучести. - М: Машиностроение. - 400 c.
8. Николаев А.П., Клочков Ю.В., Киселев А.П., Гуреева Н.А. (2009). Расчет оболочек на основе МКЭ в двумерной постановке. - Волгоград «ИПК» Нива.- 432 с.
9. Новожилов В.В. (1962). Теория тонких оболочек. - Л: Судпромиздат. - 432 с.
10. Dzhabrailov A. Sh., Klochkov Yu. V., Nikolaev A. P. The finite element analysis of shells of revolution with a branching meridian, Russian Aeronautics, Vol. 52, № 1. p. 22-29.
A FINITE ELEMENT ANALYSIS OF AXISYMMETRIC LOADED SHELLS OF REVOLUTION WITH A BRANCHING MERIDIAN UNDER ELASTIC-PLASTIC DEFORMING
Yu.V. Klochkov, A.P. Nikolaev, A.Sh. Dzhabrailov
Volgograd State Agricultural University, Volgograd
In this work, the calculation algorithm of axisymmetric loaded shells of revolution with branching meridian on the basis of the theory of small elastic-plastic deformations is studied.
KEYWORDS: shell of revolution, rigidity matrix, finite element, branching meridian.
