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

Полный текст

Введение1 Конструкции из тонких оболочек находят самое широкое применение в строительстве и архитектуре [1; 2], машиностроении, авиастроении, химической, нефтяной и газовой промышленностях и т.д. При проектировании и реконструкции такого рода объектов в настоящее время используют чис- 1 Исследование выполнено при финансовой поддержке РФФИ и Администрации Волгоградской области в рамках научного проекта № 18-41-340007 р_а. © Клочков Ю.В., Николаев А.П., Вахнина О.В., 2019 This work is licensed under a Creative Commons Attribution 4.0 International License ленные методы анализа их напряженно-деформированного состояния (НДС) с применением высокопроизводительной вычислительной техники [3- 7]. Одним из наиболее распространенных численных методов анализа НДС тонкостенных конструкций является метод конечных элементов (МКЭ) в различных формулировках [8-18]. Несмотря на значительное количество публикаций, посвященных данной проблематике, по-прежнему актуальной является задача совершенствования конечно элементных алгоритмов в плане решения проблем совместности используемых конечных элементов, повышения точности численных решений и других важных аспектов по данному направлению. ЧИСЛЕННЫЕ МЕТОДЫ РАСЧЕТА КОНСТРУКЦИЙ 117 Klochkov Yu.V., Nikolaev A.P., Vakhnina O.V. Structural Mechanics of Engineering Constructions and Buildings, 2019, 15(2), 117-126 1. Геометрические соотношения где = a a a ) a ( 0 0 0 0 2 11 22 12 - детерминант метриче- Срединная поверхность может быть задана радиус-вектором ского тензора, компоненты которого определяются скалярными произведениями r r r r r r R0 = x (θ1 ,θ2 )i + y (θ1 ,θ2 ) j + z (θ1 ,θ2 )k , (1) a0 = a0 × a0 = 1+ (r )2 ; 11 1 1 , x 0 r0 r0 где (θ1 , θ2 ) - криволинейные координаты поверхa12 = a1 × a2 = r, x × r,θ ; ности оболочки. 0 r0 r0 2 2 Выражение (1) может быть конкретизировано для оболочек различного типа. Например, для трехосного эллипсоида оно выглядит следующим образом a22 = a2 × a2 = (r,θ ) + r . Здесь под r понимается функция r ( x,θ) (8) или r r r r r (θ) , определяемая по (3) или (4). R0 = xi + r ( x,θ)sinθ j + r ( x,θ)cosθ k , (2) Производные локального базиса точки M0 где θ - угловая координата, отсчитываемая против хода часовой стрелки от вертикальной оси в поперечном сечении эллипсоида плоскостью, пермогут быть получены по деривационным формулам [19] r0 0ρ r0 0 r0 r0 0ρ r0 пендикулярной оси Ox . Входящая в (2) функция r ( x,θ) имеет вид aα,β = Гαβ aρ + bαβ a ; a,α = -bα aρ , (9) r ( x,θ) = (1-( x / a)2 )bc / c2 sin2 θ + b2 cos2 θ, (3) где 0 Г αβ 0ρ - символы Кристоффеля второго рода; 0ρ bαβ и bα - ковариантные и смешанные компогде a, b и c - параметры трехосного эллипсоида ненты тензора кривизны. Здесь и ниже греческие индексы последовательно принимают значения 1, 2. при записи его уравнения в каноническом виде В процессе деформирования оболочечной кон- (x2 a2 + y2 b2 + z2 c2 = 1). струкции точка M0 и отстоящая от нее на Если в формуле (3) первый сомножитель расстоянии ζ точка M 0ζ займут новые почислителя принять равным единице, то можно получить следующее выражение ложения M и рами M ζ , определяемые радиус-вектоr r r r r r r (θ) = bc / c2 sin2 θ + b2 cos2 θ. (4) r ρ r0 R = R0 + v; r0 Rζ = R + ζ a, (10) Заменяя в (2) функцию r ( x,θ) формулой (4), где v = v aρ + va 0 - вектор перемещения точможно получить выражение для радиус-вектора эллиптического цилиндра: r r r r ки M . Входящий во вторую формулу (10) орт нормали a точки M определяется векторным произведением R0 = xi + r (θ)sinθ j + r (θ)cosθ k .¨ (5) Дифференцированием (2) или (5) по x и θ можно получить касательные векторы локальноr r r a = a1 ´ a2 / a , r r 2 (11) го базиса в произвольной точке M0 поверхности где aα = R,α ; a = a11a22 - (a12 ) - детерминант оболочки: r r r r 0 0 0 0 метрического тензора деформированного состояния, который может быть представлен в виде a1 = R,x ; a2 = R,θ . (6) 0 ( ρ ) Орт нормали в точке M 0 определяется выраa = a 1+ 2ερ , (12) жением ε ρ где ρ - смешанные компоненты тензора дефорr 0 r 0 r 0 0 маций оболочки в точке срединной поверхности. a = a1 ´ a2 / a , (7) 118 NUMERICAL METHODS OF ANALYSIS OF STRUCTURES Клочков Ю.В., Николаев А.П., Вахнина О.В. Строительная механика инженерных конструкций и сооружений. 2019. Т. 15. № 2. С. 117-126 Дифференцированием (11) можно получить ì ü производные орта нормали к деформированной { L }T í{ 1L }T { 2 L }T { L }T ý; Uy = vy vy vy (17) поверхности оболочки 1´27 î 1´9 1´9 1´9 þ r r r r r r r G T ì 1G T 2G T G T ü a,α = (a1,α ´ a2 + a1 ´ a2,α ) a +(a1 ´ a2 )(1 ) a , (13) ,α {U y } = í{vy } {vy } {vy } ý, (18) 1´27 î 1´9 1´9 1´9 þ где 1 a » (1- ερ ) a0 . L Т i j k i j k i j k ρ Деформации в точке M ζ определяются где {qy } = {q q q q,ξ q,ξ q,ξ q,η q,η q,η }; G Т i j k i j k i j k известным соотношением механики сплошной среды [20]: {qy } ={q q q q,x q,x q,x q,θ q,θ q,θ }. Здесь под q понимается компонента вектора ζ ( 0 ) 1 2 εαβ = gαβ - gαβ / 2, (14) перемещения v , v или v. где r r r r Компоненты вектора перемещения точки внутренней области конечного элемента выражаются через их узловые значения с помощью интерпоg = g × g = Rς × Rς ; αβ α β ,α ,β ляционных зависимостей вида [7-18] r r r r r r g0 = g × g0 = (R0 + ζa0 ) ×(R0 + ζa0 ) . { }T { L }T αβ α β ,α ,β q = φ qy , (19) Соотношения (14) могут быть представлены суммой T где {φ} 1´9 9´1 - матрица-строка функций формы, со- αβ αβ αβ εζ = ε + ζÀ , r r r r (15) держащая двумерные полиномы третьей степени. Рассматриваемый треугольный конечный элемент является совместным по компонентам вектора где ε = (a 0 × v + v × a 0 ) 2 ; перемещения, но несовместным по их производ- αβ α ,β ,α β ным. Если вычислить производные нормальной ( r0 r r r0 r r0 r r r0 0ρ r0 компоненты вектора перемещения вдоль нормалей Àαβ = aα × a,β + v,α × aβ + aα × aβ + a,α × v,β + aα ×bβ × aρ к сторонам конечного элемента в точках 1, 2, 3, со- 0ρ r0 r0 ) ответствующих серединам сторон, то в значениях +bα × aρ × aβ 2. этих производных в смежных элементах I и II будет наблюдаться различие (рис. 1). 2. Треугольный конечный элемент В качестве конечного элемента был выбран фрагмент срединной поверхности тонкой оболочки треугольной формы с узлами i, j, k, расположенными в его вершинах. Для реализации процедуры численного интегрирования по площади конечного элемента треугольный фрагмент срединной поверхности отображается на прямоугольный треугольник с локальной системой координат 0 £ ξ, η £ 1. Глобальные координаты x и θ точки внутренней области конечного элемента выражаются через глобальные координаты узлов зависимостями x = (1-ξ - η) xi + ξ x j + η xk ; θ = (1- ξ - η)θi + ξθ j + ηθk . (16) Рис. 1. Векторы нормалей в смежных элементах [Figure 1. Vectors of normals in adjacent elements] Столбец узловых варьируемых параметров конечного элемента в локальной ξ,η и глобальной x, θ системах координат был выбран в следующем виде: m m¢ ¶v ¶v r ¹ r , (20) ¶nm ¶nm¢ ЧИСЛЕННЫЕ МЕТОДЫ РАСЧЕТА КОНСТРУКЦИЙ 119 Klochkov Yu.V., Nikolaev A.P., Vakhnina O.V. Structural Mechanics of Engineering Constructions and Buildings, 2019, 15(2), 117-126 где индексы m и m¢ принимают значения до- ¶v ¶v ¶v 2 полнительных узлов 1, 2, 3, и 1¢, 2¢, 3¢ (рис. 1) соответственно. r ¶n2 = cos γ + cos δ = ¶x ¶θ что при Отмеченное различие обусловлено тем, r = é( T {φ,ξ } × ξ, x + T {φ,η} o η, x )cos γ + вычислении, например, производных ¶v2 ¶n2 и êë ¶v2¢ ¶n2¢ в точках 2 и 2¢ смежных элементов I и + ({φ,ξ } × ξ,θ + {φ,η} η,θ )cos δúû {vy }, (24) r II используются значения нормальной компонен- T T ù L r ты вектора перемещения в узлах i и i¢ , находягде γ и δ - углы между нормалью n2 и векторащихся за пределами общей границы j - k ных элементов I и II. смежми локального базиса в данной точке (рис. 1). Для смежного элемента II соотношения (23) Решение отмеченной проблемы совместности треугольных конечных элементов предлагается осуществить за счет использования множителей Лагранжа. Тогда выражение (20) может быть заи (24) берутся с противоположным знаком. Если рассматривать треугольный элемент I в отдельности, то для него можно записать равенство писано в следующем виде: ¶ r 1 2 λ ¶v1 + λ n1 ¶ r 3 ¶v2 + λ n2 ¶ r ¶v3 = 0. n3 (25) æ ¶v ¶v ö λ m - m¢ = 0, (21) m ç r r ÷ ç nm nm ÷ С учетом (22)-(24) соотношение (25) можно è ¶ ¶ ¢ ø представить в матричном виде где λm - множитель Лагранжа в дополнительном é ù ì¶v1 ü r узле m . 1 T ï¶n ï ï ï {d1} [PR ] ê 1´27 27´27 ú Входящие в (21) производные нормальной ï¶v ï ê ú G G компоненты вектора перемещения вдоль норма- {λ}T r2 ={λ}T {d }T [P ] {U } ={λ}T [D]{U } = 0, (26) í ý 2 R y y r 1´3 ï¶n2 ï 1´3 ê 1´27 27´27 ú 27´1 1´3 3´27 27´1 ê T ú лей к серединам сторон nm могут быть выраже- ¶v3 {d } [P ] ны через столбцы (17) и (18) узловых неизвест- ï r êë úû ï ï î¶n þï 3 1´27 R 27´27 ных треугольного конечного элемента: 3 3´27 3´1 T ¶vm = {d }Т {U L } = {d }Т [P ]{U G }, (22) где {λ} = {λ1 λ2 λ3}. r m y m R y ¶nm 1´27 27´1 1´27 27´27 27´1 Функционал Лагранжа для треугольного элемента с учетом дополнительного условия (26) загде [PR ] } y - матрица перехода от столбца {U L к пишется в виде T T = { ς } { αβ} V +{ } [D]{U G }y столбцу {U G}. ФL ò V εαβ o d λ y Структура входящих в (22) матриц-строк {dm } зависит от ориентации треугольных конечных элементов на срединной поверхности обо- -ò {U}T {P}dF, F T T (27) лочки. Если сетку дискретизации сориентировать где { ς } { ς ς ς } { αβ } { 11 22 12} εαβ = ε11 ε22 2ε12 ; σ = σ σ σ - вдоль линий главных кривизн (рис. 1), то для матрицы-строки деформаций и напряжений в точке треугольного элемента I будут справедливы зависимости M ζ ; {U}T = {v1 v2 v} - матрица-строка компо- ¶v1 ¶v ( = - {φ }T ×ξ +{φ }T ×η ){vL }; нент вектора перемещения точки M ;{P} - стол- ¶ r =n1 ¶θ ,ξ ,θ ,η ,θ y бец внешней поверхностной нагрузки. Столбец {σαβ } на основании закона Гука [20] ¶v ( ¶v3 = - {φ }T ×ξ +{φ }T ×η ){vL }. (23) может быть представлен матричным произведением ¶ r =n3 ¶x ,ξ ,x ,η ,x y { αβ } [ ]{ ς } Для узла 2 на стороне j - k треугольного ко- σ 3´1 = C 3´3 εαβ , 3´1 (28) нечного элемента можно записать следующее выражение: где [C] - матрица упругости. 3´3 120 NUMERICAL METHODS OF ANALYSIS OF STRUCTURES Клочков Ю.В., Николаев А.П., Вахнина О.В. Строительная механика инженерных конструкций и сооружений. 2019. Т. 15. № 2. С. 117-126 Деформации в произвольном слое оболочки, отстоящем на расстоянии ζ от срединной поверхности, с учетом (15) и (19) могут быть представлены в виде произведения матриц: 3. Пример расчета Была рассчитана эллипсоидальная оболочка, поверхность которой описывается радиус-вектором (5). Эллиптический цилиндр загружен в середине дву- { ς } [ ]{ } εαβ = Г εαβ = мя диаметрально противоположными сосредоточенными силами Q (рис. 2). Приняты следующие [ ][ ]{ L } [ ][ ][ ]{ G }. исходные данные: Q = 453, 6 Н; длина образую- = Г B Uy = Г B PR Uy (29) 4 С учетом (28) и (29) функционал (27) преобщей L = 26см; модуль упругости 7, 38 ´10 МПа; разуется к виду коэффициент Пуассона n = 0,3125 ; толщина Ф ={UG}T [P ]T ò[B]T [Г]T [C][Г][B]dV [P ]{UG}+ стенки t = 0,24 см; параметры эллиптического L y R R y V поперечного сечения b = 12,58 см, c = 11,43 см. y y R +{λ}Т [D]{U G }-{UG }T [P ]T ò [ A]T {P}dF , (30) В силу наличия плоскостей симметрии рассчитана 1/8 часть эллиптического цилиндра. é{φ}T {0}T F {0}T ù ê ú где [ A] = ê{0}T {φ}T {0}T ú - матрица функ- 3´27 ê ú ê{0}T {0}T {φ}T ú ë û ций формы. y Выполняя минимизацию (30) по {U G }T и {λ}T , получим следующую систему уравнений: ì ¶ФL ï G T = [К ]{U G }+ [D]Т {λ}-{F} = 0; ï¶{U y } í 27´27 y 27´1 27´3 3´1 27´1 (31) Рис. 2. Расчетная схема эллиптического цилиндра ï ¶ФL ï T î¶{λ} y = [D]{U G } = 0, 3´27 27´1 [Figure 2. Elliptical cylinder design] Расчеты проводились в двух вариантах: R где [К ] = [P ]T ò [B]T [Г]T [C][Г][B]dV ; 27´27 V R {F} = [P ]T ò [ A]T {P}dF . в первом варианте использовался треугольный конечный элемент без множителей Лагранжа с матрицей жесткости [К ] порядка 27 ´ 27 ; во втором варианте применялся совместный элемент с мат- 27´1 F Систему (31) можно записать в более комрицей жесткости [Кλ ] размерностью 30 ´ 30 . пактной форме: λ} = {Fλ }, G [Кλ ]{U y (32) В качестве контрольного варианта использован четырехугольный конечный элемент также с девятью степенями свободы в узле (17), (18) с по- 30´30 30´1 30´1 рядком матрицы жесткости 36 ´ 36 [21]. Резуль- é[К ] [D]T ù таты повариантных расчетов представлены в таблице, в которой приведены численные значения где [К ] = ê27´27 27´3 ú; нормальных напряжений σxx и σθθ на внутрен- λ 30´30 ê ú ê[D] [0] ú ней σin и наружной σout поверхностях оболочки ë 3´27 3´3 û в точке N с координатами x = L 2; θ = π 2 , G T ì G T T ü Т ì ü {U y λ} = í{U y } {λ} ý; {Fλ } = í{F}{00 0}ý. а также величины прогиба v под сосредоточенной 1´30 î 1´27 1´3 þ 1´30 î 1´27 1´3 þ силой Q в зависимости от густоты сетки дискретизации. ЧИСЛЕННЫЕ МЕТОДЫ РАСЧЕТА КОНСТРУКЦИЙ 121 Klochkov Yu.V., Nikolaev A.P., Vakhnina O.V. Structural Mechanics of Engineering Constructions and Buildings, 2019, 15(2), 117-126 Значения напряжений и прогиба [Table. Stress and deflection values] Таблица Напряжения (МПа) в точке N с координатами x = L 2 , θ = π 2 [Stress (MPa) in N point with coordinates x = L 2, θ = π 2 ] Вариант расчета [Variant of calculation] Четырехугольный КЭ [Quadrilateral CE] [21] I II Сетка дискретизации [Sampling grid] 6×6 7×7 8×8 6×6 7×7 8×8 7×7 σin xx 8,67 12,12 14,04 -11,59 -10,87 -10,36 -15,22 σout xx -13,23 -15,67 -17,06 10,47 10,56 10,45 12,02 σin θθ -21,15 -17,07 -14,55 -47,01 -46,18 -45,52 -46,95 σout θθ 17,31 15,02 13,40 37,05 38,16 38,47 40,68 Прогиб под силой Q [Deflection under force Q] v × 10-2 м -0,2997 -0,2985 -0,2968 -0,2651 -0,2662 -0,2666 -0,2771 Как видно из таблицы, численные значения нормальных напряжений кардинально различаются между собой в зависимости от варианта расчета. Значения прогиба под сосредоточенными силами в первом варианте расчета оказались завышенными примерно на 12 % по сравнению со xx Так, в первом варианте расчета σin на внутренней вторым вариантом расчета. поверхности оказались растягивающими, а на наружной - сжатыми. В действительности же, если проанализировать деформированное состояние обо- В крайней правой колонке таблицы представлены численные значения нормальных напряжений и прогиба, полученные при использовании лочки в сечении, перпендикулярном оси Ox, прочетырехугольного конечного элемента 36 ´ 36 при ходящем через точки приложения сил Q (рис. 3), можно отметить, что внутренняя поверхность эллиптического цилиндра в точке N сжата, а наружная - растянута, что и наблюдается во втором варианте расчета. сетке дискретизации 7 ´ 7 . Сопоставляя результаты повариантных расчетов со значениями крайней правой колонки, можно отметить следующее. Напряжения σxx , полученные при использовании четырехугольного элемента 36 ´ 36 , имеют такие же знаки, что и напряжения второго варианта расчета. Напряжения σθθ , которые существенно больше, чем σxx во втором варианте расчета, практически совпали или достаточно близки по своим значениям, полученным при применении четырехугольного элемента. Значения σxx в первом варианте можно признать неприемлемыми из-за несоответствия картине деформирования оболочки, а значения σθθ оказались в 4 раза заниженными по сравнению со вторым и контрольным вариантами расчета. Заключение Рис. 3. Деформация цилиндра [Figure 3. Cylinder deformation] На основании выполненного анализа НДС тонкой оболочки в виде эллиптического цилиндра, за- 122 NUMERICAL METHODS OF ANALYSIS OF STRUCTURES Клочков Ю.В., Николаев А.П., Вахнина О.В. Строительная механика инженерных конструкций и сооружений. 2019. Т. 15. № 2. С. 117-126 груженного системой двух сосредоточенных сил, можно заключить, что использование несовместных треугольных элементов приводит к существенным погрешностям расчета, вплоть до неприемлемых. Для корректного анализа НДС тонких оболочек необходимо использовать совместный треугольный элемент, матрица жесткости которого формируется в соответствии с алгоритмом, изложенным в настоящей статье.

Об авторах

Юрий Васильевич Клочков

Волгоградский государственный аграрный университет

Автор, ответственный за переписку.
Email: klotchkov@bk.ru
SPIN-код: 9436-3693

доктор технических наук, профессор, заведующий кафедрой высшей математики, Волгоградский государственный аграрный университет. Опубликовал 165 научных статей, 4 монографии, 4 наименования учебно-методической литературы

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Анатолий Петрович Николаев

Волгоградский государственный аграрный университет

Email: anpetr40@yandex.ru
SPIN-код: 2653-5484

доктор технических наук, профессор, профессор кафедры прикладной геодезии, природообустройства и водопользования, Волгоградский государственный аграрный университет. Опубликовал 149 научных статей, 6 монографий, 5 наименований учебно-методической литературы

Российская Федерация, 400002, Волгоград, Университетский пр., 26

Ольга Владимировна Вахнина

Волгоградский государственный аграрный университет

Email: ovahnina@bk.ru
SPIN-код: 3593-0159

кандидат технических наук, доцент кафедры высшей математики, Волгоградский государственный аграрный университет. Опубликовала 47 научных статей, 1 монографию, 8 наименований учебно-методической литературы

Российская Федерация, 400002, Волгоград, Университетский пр., 26

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