CC BY
18
КОНЕЧНО-ЭЛЕМЕНТНЫЙ АНАЛИЗ НДС ОБОЛОЧЕК ВРАЩЕНИЯ С УЧЕТОМ ДЕФОРМАЦИЙ ПОПЕРЕЧНОГО СДВИГА1
Ю.В. КЛОЧКОВ, доктор техн. наук, профессор А.П. НИКОЛАЕВ, доктор техн. наук, профессор ТР. ИЩАНОВ, аспирант
Волгоградский государственный аграрный университет 400002, г. Волгоград, пр. Университетский, 26, e-mail: ishchanov. volgau@yandex. ru
В данной работе на основе четырехугольного элемента дискретизации излагается алгоритм конечно-элементного расчета оболочек вращения с учетом деформаций поперечного сдвига при различных вариантах отсчета угла наклона нормали в процессе деформирования.
КЛЮЧЕВЫЕ СЛОВА: скалярная аппроксимация, конечный элемент, круговой цилиндр, поперечный сдвиг.
При расчете оболочек вращения наиболее часто используется теория тонких оболочек, основанная на гипотезах Кирхгофа-Лява [1-3]. Однако, в ряде случаев (например, при расчете короткопролетных конструкций), пренебрежение деформациями сдвига не является вполне корректным.
В таких ситуациях, как правило, используются теории оболочек типа Тимошенко [4-7]. Решение систем дифференциальных уравнений, описывающих процесс деформирования оболочек вращения с учетом деформации поперечного сдвига аналитическими способами весьма затруднительно. Поэтому в настоящее время используются численные методы расчета, как правило, метод конечных элементов (МКЭ) [8-13].
1. Геометрия оболочки
Срединная поверхность оболочки вращения может быть задана радиус-вектором:
ft0 = xi + r • sin(i)/+ г • cos(i) /с, (1.1)
где х - осевая координата; г = г(х) - радиус вращения; t - угловой параметр, отсчитываемый от оси OZ против хода часовой стрелки.
Дифференцированием (1.1) по дуге меридиана и дуге окружности s2 = rdt можно получить касательные орты локального базиса:
= JR° = хд? + Гд • sin(i)/+ Гд • cos(i) к;
go = = cos(t)/- sin(t) /с, (1.2)
где нижние индексы 1, 2 после запятой обозначают операцию дифференцирования по криволинейным координатам s-l и s2 соответственно.
Орт нормали к срединной поверхности определяется векторным произведением
п0 = х е° = — Гд? + Хд • sin(t)/+ хд • cos(t) к, (1.3)
где
1
i + ы2
1 Работа выполнена при поддержке гранта РФФИ № 15-41-02346 р_поволжье_а.
48
Орты локального базиса (1.2) и (1.3) и их производные по координатам и я2 могут быть представлены матричным выражением [14]:
(е0) = К]{Г}; = Ы(е0}; (1.4)
где
(ё°)т = {ё°1ё°2г10); (1)т = {Т/*}.
Вектор перемещения точки срединной поверхности оболочки вращения и его производные по глобальным координатам я2, при учете (1.4), могут быть представлены компонентами, отнесенными к локальному базису данной точки:
V = У1*?! + Р2в2 + УП0; = 1 + + £<хп°;
= + + 1ацП0, (1.5)
где а,р последовательно принимают значения 1, 2; V1, г>2,г> - тангенциальные и нормальная компоненты вектора перемещения; Ьа, Ь^р, ^р, - много-
члены, содержащие компоненты вектора перемещения, их первые и вторые (для Ьар) производные по глобальным координатам з2. Положение точки, отстоящей от срединной поверхности на расстоянии £ в исходном и деформированном состояниях, определяется соответствующими радиус-векторами:
до<; = до + ^0; я? = до<; + £ (1.6)
Входящий в (1.6) вектор перемещения точки, отстоящей от срединной поверхности на расстоянии может быть определен следующим образом:
у = г? + ^(у х п0), (1.7)
где
у = -у2е£ + у1«?? + Л„гг0
- вектор углов поворота нормали [7].
Соотношение (1.7) можно рассматривать как вариант, в котором вектор углов поворота нормали отсчитывается от исходного положения нормали. Данный вариант соответствует подходу, описанному в [5,6].
Если поворот нормали отсчитывать от ее деформированного состояния, то формулу (1.7) следует записывать в следующем виде:
У = г? + ^(г?п х п0) + х п0), (1.8)
где г?п - вектор разности нормалей в деформированном и исходном состояниях [7,14]:
гЯ = Й-п0 = -рп2е? + рп1е% + рпп0. (1.9)
Здесь п - орт нормали в деформированном состоянии
г1 = а1х а2/1а1х а21, (110)
где
аа = Яа = + г?) а.
Ковариантные компоненты тензора деформаций определяются соотношением механики сплошной среды [15]:
4з = (ё*3-Ваз)П, (111)
где gaß, ёаз, Saß, ёаз - ковариантные компоненты метрического тензора в деформированном и исходном состояниях, определяемые соответствующими скалярными произведениями
gaß=ga^gß'; ga3=ga^g3; gaß=ga^gß; ga3 = • (112)
Входящие в (1.12) векторы базиса ga,g3,ga,g3 могут быть найдены дифференцированием соответствующих радиус-векторов
= =
Й = ^ = (1.13)
Соотношения (1.11) могут быть представлены в виде суммы:
£aß = £aß + (Haß;
£аЗ = £аЗ + Í*a3, (114)
где saß; xaß - деформации и искривления срединной поверхности оболочки
вращения; £^3, £а3 - деформация сдвига в точках с радиус-векторами R и Ж;
£зз, £зз - линейные деформации вдоль нормали в тех же точках.
2. Физические соотношения тонких оболочек
Контравариантные компоненты тензора напряжений в произвольном слое оболочки, отстоящем от срединной поверхности на расстоянии (, определяются через ковариантные компоненты тензора деформаций соотношениями механики сплошной среды [15]:
атп = ÄI1(£)gmn + 2 MgmrgnP£rp, (2.1)
где верхние и нижние индексы т, п, у, р последовательно принимают значения 1, 2, 3; А,^ - коэффициенты Ляме; gmn - контравариантные компоненты метрического тензора;
11(е) = ётпе1п - первый инвариант тензора деформаций.
Соотношение (2.1) может быть представлено в матричном виде
[о™}=[С]{еЦ, (2.2)
6X1 6X6 6Х/
где
\Т
{arnn}T = (alla12a13a22a23a33}; {^J = {^^^f^^}.
1X6 1X6
Принимая во внимание общепринятую в теории тонких оболочек гипотезу а33 = 0 [1-3], из (2.1) можно получить следующую зависимость:
£33 = / (£11'£12'£13'£22'£2З)- (2 3)
Используя зависимость (2.3), можно уменьшить размерность матричного выражения (2.2):
{а^}=[С]{£Ц (2.4)
5X1 5х5^5х7
где [С] - матрица упругости.
5X5
3. Четырехугольный конечный элемент
В качестве элемента дискретизации выбирается четырехугольный фрагмент срединной поверхности оболочки вращения с узлами I,], к, I, отображае-
мый для удобства численного интегрирования на квадрат в локальной системе координат —1 < < 1.
Столбец узловых неизвестных конечного элемента в глобальной и локальной системах координат выбирается в виде:
=}тн}тштшт\ (3.1)
1х44 V 1х12 1х12 1х12 1х4 1х4 )
{^Г = ) {<}Т№}Т№}Т{у}}Т{У1}Т [, (3.2)
1х44 V 1х12 1х12 1х12 1х4 1х4
где
{Яу}Т = {я1 -У^д -9,19,2 -9,2};
1х12
= - ^Чс - Ч^Ч - 9/77};
1х12
{Уу}Т = {у1у}укУ11
1х4
Здесь под дт (т = г, у, к, I) понимается компонента вектора перемеще-
ния Рат или Рт.
Для вычисления компонент вектора перемещения д и компонент векто-
ра углов поворота нормали уа используются интерполяционные выражения следующего вида:
ч = {<РГ ка = {^)т{куа}т, (3.3)
1х12 12х1 1х4 4х1
где = {^1^2 - ^12) - матрица строка, содержащая произведения по-
1х12
линомов Эрмита третьей степени;
{'Ф)Т = - матрица-строка,
1х4
содержащая билинейные функции локальных координат.
Дифференцированием (3.3) можно получить производные компонент вектора перемещения и компонент углов поворота нормали:
Ч,а = + {ч^
= + (3.4)
Для получения матрицы жесткости и столбца узловых усилий четырехугольного конечного элемента можно воспользоваться функционалом Лагран-жа:
I {аУР)М = I {и)т{Р)(1Р, (3.5)
где
(£Ур) {£112£122£13£222£2з} ; (3 6)
= {а11а12а13а22а23};
{У)т = {у1^2^}, {Р)т = {Р1Р2Р} - столбцы компонент вектора перемещения и внешней нагрузки.
С учетом (1.14), (2.4), (3.3) функционал (3.5) может быть преобразован к виду
(3.7)
{УуГ}>*г[ [В]г [Г]г [С] [Г] [В] ^[Рк]{иуГ} =
1x44 44X44-/ 44X10 10x5 5x5 5x10 10x44 44x44 44x1
= {иТу}Т[РК]Т [ №
1X44 44x44 У 44x3 3x1
где [Г] - матрица перехода от столбца (3.6) к столбцу
5x10
{£а
= {£112£122£13£222£23^112х122х13х222х23};
1x10
[В] - матрица дифференциальных и алгебраических операторов, необходимая
10x44
для перехода от столбца {еар\ к ст
олбцу {Щ} - (3.2); [РК] - матрица перехода
44x44
от столбца {иу} к столбцу {^у}, определяемая с использованием соотношений
= + (3.8)
Входящая в (3.7) матрица [Л] имеет следующую структуру:
\шт {о}г {0}т {0}г {о}7]
1Т12 1x12 1x12 1x4 1x4
Г {0}г {0}г
1x12 1x4 1x4 • \Т ГП\Т ГП\Т
[А] =
3x44
{0}г {ф}т {0}т {0}т {0}
К J Чч/ЧО -1ч/ Л. 1 VI
1x12 1x12
{0}г {0}г {#>}г {0}т {0}
1-1x12 1x12 1x12 lx4
(3.9)
В результате минимизации (3.7) по {^у} можно получить следующее матричное выражение
[кГ]{Щ}= [рг], (3.10)
44x44 44x1
где
[КГ] = [РК]Т Г [В]т [Г]т [С] [Г] [В] (IV [Рк] - матрица жесткости;
44x44 44x44 44x10 10x5 5x55x1010x44 44x44
Ю = [^г [ № {РЫР
44x1 44x44 J 44x3 3x1
- столбец узловых усилий четырехугольного конечного элемента в глобальной системе координат.
Рис. 1 52
Пример расчета 1. В качестве примера была рассчитана цилиндрическая оболочка, загруженная вдоль образующей распределенной нагрузкой интенсивности q и имеющая на диаметрально противоположной образующей шарнирные опоры, препятствующие вертикальному смещению (рис. 1).
Вследствие наличия плоскостей симметрии оболочка моделировалась одной лентой конечных элементов, ориентированной в кольцевом направлении. Были приняты следующие исходные данные:
Расчеты были выполнены в двух вариантах:
в первом варианте при формировании матрицы жесткости КЭ использовалось соотношение, соответствующее отсчету угла поворота нормали от ее исходного положения (1.7);
во втором варианте отсчет угла поворота нормали осуществлялся от ее деформированного положения - (1.8) .. (1.10).
Результаты повариантного расчета представлены в виде диаграммы (рис. 2), на которой приведены значения физических напряжений в точке приложения нагрузки на внутренней (верхняя часть диаграммы) и внешней (нижняя часть диаграммы) поверхностях цилиндра в зависимости от числа элементов дискретизации .
Как видно из диаграмм, сходимость вычислительного процесса во втором варианте существенно лучше, чем в первом, а численные значения напряжений близки к значению , вычисленному по формуле сопротивления ма-
териалов [16] для задачи расчета кольца с двумя сосредоточенными силами.
250 т
м
" 200 о
зН 150 5 100
&1 и
5> 0 я £1
05 -50
5 -100
я -150 ГО
-200
50
-250
-17
-18
2,4
0,3
8,5
180,2
-18
9,9
-18
-1
5,1
90
145
-186,8
-190
9,5
-18
193
7,6
90
1 вариант расчета
2 вариант расчета
1 вариант расчета
2 вариант расчета
2
5
4
9
9
7
1
Рис. 2
При сопоставлении результатов повариантного расчета (таблица 1) следует также отметить наличие в первом варианте расчета «скачка» в значениях напряжений в узлах смежных элементов дискретизации (Рис. 1). Величина данного «скачка» уменьшается при увеличении числа элементов дискретизации. Во втором варианте расчета вышеупомянутый «скачок» в значениях практически не наблюдается.
Таблица 1
пэ Я/см*\ 25 49 97 145 193
1 вариант °22 172,4 -170,3 182,1 -180,2 187,0 -185,1 188,7 -186,8 189,5 -187,6
133,1 -131,9 162,3 -160,8 177,1 -175,4 182,0 -180,3 189,5 -187,6
н 153,4 -151,3 172,3 -170,6 182,1 -180,3 185,4 -183,5 187,0 -185,1
& И ^2 152,8 -150,7 172,3 -170,6 182,1 -180,3 185,4 -183,5 187,0 -185,1
Пример расчета 2. Был рассчитан жестко защемленный по торцам цилиндр, нагруженный внутренним давлением интенсивности q (рис. 3). Были приняты следующие исходные данные:
Я = 1,0 м; I = 1,0 м; Е = 2 • 105 МПа; V = 0,3; д = 5 МПа.
Вследствие наличия осевой симметрии оболочка моделировалась одной лентой КЭ, ориентированной вдоль образующей.
Расчеты, как и в примере 1, выполнялись в двух вариантах. Результаты по-вариантного расчета представлены в виде диаграммы (рис. 4), на которой показаны значения физических напряжений а11 в жесткой заделке на внутренней (верхняя часть диаграммы) и внешней (нижняя часть диаграммы) поверхностях цилиндра в зависимости от числа элементов дискретизации пэ.
Анализ диаграммы показывает существенно лучшую сходимость вычислительного процесса во втором варианте расчета по сравнению с первым вариантом.
Для достижения аналогичного уровня точности в первом варианте требуется на порядок большее число элементов дискретизации, чем во втором варианте.
Ь
зЯ Я Я
£ К ft
600
500
400
300
372
200
Я Я
Я 500
0
к я я
й -100 я
Я -200 го
-300
-400
,2
,7
349,5-35
48
7 455
75 ,1460
,4-340,4 354
2 вариант расчета
1 вариант расчета
2 вариант расчета
1 вариант расчета
3
1
1
2
4
3
6
4
8
0
9
6
1
2
1
8
Рис. 4
Анализ результатов вычисленных напряжений, представленных на обеих диаграммах, позволяет сделать вывод о предпочтительности второго варианта расчета, соответствующего отсчету угла поворота нормали от ее деформированного состояния, при котором можно получать удовлетворительные по точности значения напряжений при относительно редкой сетке элементов дискретизации.
В то же время следует отметить, что и первый вариант компоновки матрицы жесткости четырехугольного конечного элемента, при котором отсчет угла поворота нормали осуществлялся от ее исходного состояния, позволяет получать приемлемые результаты при весьма существенном сгущении сетки дискретизации рассчитываемой оболочки.
Л и т е р а т у р а
1. НовожиловВ.В. Теория тонких оболочек. - Л.: Судостроение, 1962. - 431с.
2. Голованов А.И., Тюленева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. - М.: Физматлит, 2006. - 391 с.
3. Скопинский В.Н. Напряжения в пересекающихся оболочках. - М.: Физматлит, 2008. - 399 с.
4. Тимошенко С.П., Войновский-Кригер С. Пластинки и оболочки. - М.: Наука, 1966. - 636 с.
5. Рикардс Р. Б. Метод конечных элементов в теории оболочек и пластин. - Рига: Зинатне, 1988. - 248 с.
6. ВольмирА.С. Гибкие пластинки и оболочки. - М.: Гос. изд-во техн.-теорет. лит., 1956. - 419 с.
7. Черных К.Ф. Линейная теория оболочек. Ч. 2. Некоторые вопросы теории. -Ленинград: Изд-во ЛГУ, 1964. - 394 с.
8. Bathe K.J., Zhang H., Ji S. Finite element analysis of fluid flows fully coupled with structural interactions // Computers & Structures., 1999. - T. 72, № 1-3 p. 1-16.
9. Косицын С.Б., Ч.С. Линь. Численный анализ напряженно- деформированных состояний пересекающихся цилиндрических оболочек обделок тоннелей, взаимодействующих с окружающим массивом грунта, с учетом последовательности их возведения // IJCCSE. - 2015. - Т. 11. - № 2. - C. 101 - 106.
10. Dodla S., Bertram A., Kruger M. Finite element simulation of lamellar copper-silver composites // Comput. Mater. Sci. - 2015. - Vol. 101. - P. 29 - 38.
11. Mozolevski I., Prudhomme S. Goal-oriented error estimation based on equilibrated-flux reconstruction for finite element approximations of elliptic problems // Comput. Meth. Appl. Mech. and Eng. - 2015. - Vol. 288. - P. 127 - 145.
12. Kang Z., Bui T.Q., Hirose S., Nguyen D.D., Saitoh T. An extended consecutive-interpolation quadrilateral element (xcq4) applied to linear elastic fracture mechanics // Acta Mech. - 2015. - Vol. 226. - № 12. - P. 3991 - 4015.
13. Liu Y., Zhang H., Zheng Y., Zhang S., Chen B. A nonlinear finite element model for the stress analysis of soft solids with a growing mass // Int. J. Solids and Struct. - 2014. - Vol. 51. - № 17. - P. 2964 - 2978.
14. Николаев А.П., Клочков Ю.В., Киселев А.П., Гуреева Н.А. Векторная интерполяция полей перемещений в конечно-элементных расчетах оболочек: Монография. -Волгоград: ФГБОУ ВПО Волгоградский ГАУ, 2012. - 264 с.
15. Седов Л.И. Механика сплошной среды. - М.: Наука, 1976. - Том.1. - 536 с.
16. Беляев Н.М. Сопротивление материалов. - М.: Наука, 1976. - 608 с.
References
1. Novozhilov, V.V. (1962). The Theory of Thin Shells, L .: Sudostroeniye, 431p.
2. Golovanov, A.I., Tuleneva, O.N., Shigabutdinov, A.F. (2006). The Finite Element Method In Statics and Dynamics of Thin-Walled Structures, Moscow: Fizmatlit, 391 p.
3. Skopinsky, V.N. (2008). Stresses in the Intersecting Shells, M .: Fizmatlit, 399 p.
4. Timoshenko, S.P., Voynovskiy - Krieger, S. (1966 ). Plates and shells, M .: Nauka, 636 p.
5. Rickards, R.B. (1988). The Finite Element Method in the Theory of Shells and Plates, Riga: Zinatne, 248 p.
6. Volmir, A.S. ( 1956J. Flexible Plates and Shells, Moscow: GITTL, 419 p.
7. Chernykh, K.F. ( 1964). The Linear Theory of Shells, Part 2: Some Questions of the Theory, Leningrad: Leningrad State University, 394 p.
8. Bathe, K.J., Zhang, H., Ji, S. (1999). Finite element analysis of fluid flows fully coupled with structural interactions, Computers & Structures, Vol. 72, No 1-3, p. 1-16.
9. Kositsyn, S.B., T.X. Lin (2015). Numerical analysis of stress - strain state of intersecting cylindrical shells lining tunnels, interacting with the surrounding array of soil, with taking into account the sequence of their construction, IJCCSE, Vol. 11, No 2, p. 101 - 106.
10. Dodla, S., Bertram, A., Kruger, M. (2015). Finite element simulation of lamellar copper-silver composites, Comput. Mater. Sci., Vol. 101, p. 29 - 38.
11. Mozolevski, I., Prudhomme, S. (2015). Goal-oriented error estimation based on equilibrated-flux reconstruction for finite element approximations of elliptic problems, Comput. Meth. Appl. Mech. and Eng., Vol. 288, p. 127 - 145.
12. Kang, Z., Bui, T.Q., Hirose, S., Nguyen, D.D., Saitoh, T. (2015). An extended consecutive-interpolation quadrilateral element (xcq4) applied to linear elastic fracture mechanics, Acta Mech., Vol. 226, № 12, p. 3991 - 4015.
13. Liu, Y., Zhang, H., Zheng, Y., Zhang, S., Chen, B. (2014). A nonlinear finite element model for the stress analysis of soft solids with a growing mass, Int. J. Solids and Struct., Vol. 51, No 17, p. 2964 - 2978.
14. Nikolaev, A.P., Klotchkov, Y.U., Kiselev, A.P., Gureeva, N.A. (2012 ). Vector interpolation displacement fields in finite-element calculations of shells: a monograph, Volgograd: Volgograd VPO GAU , 264 p.
15. Sedov, L.I. ( 1976). Continuum Mechanics, Moscow.: Nauka, Vol. 1, 536 p.
16. Belyaev, N.M. ( 1976). Strength of materials, Moscow: Nauka, 608 p.
FINITE ELEMENT ANALYSIS 0F STRESS-STRAIN STATE OF SHELLS OF REVOLUTION WITH TAKING INTO ACCOUNT THE STRAIN OF TRANSVERSAL SHEARING
Yu.V. Klochkov, A.P. Nikolaev, T.R. Ischanov
In this paper, based on the Foursquare element sampling algorithm is presented finite element calculation of shells of revolution considering transverse shear strains in different variants of the reference tilt angle of the normal process of deformation .
KEYWORDS: the scalar approximation, finite element, circular cylinder , transverse shear.
