FINITE ELEMENT ANALYSIS OF THE SHEAR FLEXIBLE CYLINDRIC SHELL

详细

The finite element of a shell with low value of transverse shear stiffness is considered. As the basic kinematic variables transverse shear strains are in eveidence. The results of numerical analysis showing urgency of the model are presented.

全文:

Композитные конструкции, обладающие высокими удельными механическими характеристиками, находят широкое применение в производстве авиационной и ракетно-космической техники. Ряд особенностей композитных материалов резко выделяет их из ряда традиционных конструкционных материалов, что вынуждает разрабатывать новые расчетные модели, в которых адекватным образом учитываются эти особенности. Главная среди них - низкая сдвиговая жесткость по отношению к трансверсальным напряжениям. Учет данного обстоятельства приводит к повышению порядка разрешающих уравнений за счет введения в рассмотрение углов трансверсального сдвига. В настоящей работе рассматривается новый четырехугольный конечный элемент оболочки, среди основных узловых переменных которого присутствуют углы трансверсального сдвига. Разработанная модель основывается на математической теории пластин Рейсснера-Мидлина. Вывод основополагающих уравнений и матриц теории метода конечных элементов выполняется вариационным методом, который предполагает использование в качестве исходного выражения полной потенциальной энергии для элемента оболочки. Основные преобразования при этом касаются выражения потенциальной энергии деформации, которое в геометрически нелинейной постановке для сдвиговой модели имеет вид + Nap(Sap +8ра + ЮаЮр) + МаХа + МрХр + (1) + Мав (Хар + Хра ) + Qa Va + QpVp } dа dР , где N, Q, M - внутренние погонные усилия; е, х -обобщенные деформации начальной поверхности, ю -углы поворота бесконечно малого элемента оболочки. Рассмотрим цилиндрическую оболочку радиусом R, нагруженную давлением p. Криволинейную систему координат свяжем с начальной поверхностью, где ось x направлена вдоль образующей цилиндра, ось y -по окружности начальной поверхности, ось z -по нормали к ней. Запишем функционал Лагранжа для задачи статического расчета в геометрически линейной постановке. Выражение потенциальной энергии деформации (1) в данном случае примет следующий вид: + My Х y + Mxy Х xy + Qx V x + Qy V y ) dxdy, где внутренние погонные усилия связаны с обобщенными деформациями начальной поверхности с помощью следующих физических соотношений: Nx = B11Sx + B12Sy + C11xx + C12Xy , Ny = B21Sx + B22Sy + C21Xx + C22Xy , Nxy = B33Sxy + C33Xxy , (3) Mx = C11S x + C12S y + D11x x + D12X y My = C21Sx + C22S y + D21Xx + D22Xy , Mxy = C33Sxy + D33Xxy , Qx = KxVx , Qy = Ky Vy , где B, C, D и K - жесткостные параметры: В определяют мембранную жесткость оболочки, С - так называемая смешанная жесткость, D - изгибная жесткость, К - жесткость оболочки при трансверсальном сдвиге. Обобщенные деформации определяются следующими геометрическими соотношениями: 1 dv d2 w V2 du dv + D2 + B3 dy R dy dy2 dy dx 2 d w d\Vx ^V y 1 dv —— + —- +---2dy dx R dx dxdy du dv dy dx + 2C33 \ — + — Я 2 V d w ( d\Vx cty y 1 dv , —LJL + —L +---2dy dx R dx dxdy + Kx V2 + Ky V2} dxdy. du dx ’ dx dv w = dy + R 50 y du dv 8 xy = ay+dx d0x d0y + D 33 Х x =■ dy dy dx dw Здесь a и b определяют границы элемента оболочки вдоль оси x, а c и d - вдоль оси y. Потенциал внешних сил для цилиндрической оболочки, нагруженной давлением p, определяется следующим интегралом: dw dx = V x В выражении (4) фигурируют следующие кинематические параметры: 0x и 0y - углы наклона сечений, измеряемые в плоскостях xz и yz соответственно, u и v -перемещения точек начальной поверхности вдоль осей x и y соответственно, w - прогибы точек начальной поверхности, yx и yy - углы сдвига или осреднен-ные по высоте сечения деформации трансверсального сдвига, вычисляемые по следующим формулам [1]: h-s л h-s b d П = J!pwdxdy . (8) Выражение полной потенциальной энергии цилиндрической оболочки имеет вид E = U + П . (9) Здесь потенциальная энергия деформации задается интегралом (7), а потенциал внешних сил - интегралом (8). Функционал (9) позволяет получить разрешающие уравнения теории МКЭ. Будем рассматривать модель с сеткой прямоугольных конечных элементов с внутренней нумерацией узлов, осуществляемой при обходе контура против часовой стрелки. Вектор узловых кинематических параметров 6е имеет следующий вид: 6е = {61 62 63 б4>г, (10) где 6; (i = 1, 2, 3, 4) - вектор кинематических параметров i-го узла: -i II— о h i . dz, eyZ dZ , (5) V x = V y = где h - толщина оболочки, (h - s) и (-s) - координаты по высоте верхней и нижней поверхности оболочки соответственно. С учетом соотношений (4) для углов наклона сечений (0x и 0y) представим выражения для кривизн в виде d 2 w dx2 dv ^2 dv dvx. dx y 1 dv d w dy R dy dy2 Х x =■ Х y =■ dV y 1 dv , d2 w T y +- — - 2- R dx dxdy (6) + dy dx Х xy 6; = * \Ц [f . (11) Подставим физические (3) и геометрические (4), (6) соотношения в подынтегральное выражение (2). В результате получим bd ( /^\2 Кинематические параметры можно выразить через их узловые значения следующим образом: [ dw V [ dw V w = pUw1 + pi2\ -r\ + Лэ\ -r-I + p18w2 + ~,2 \ d w 2 JJ{b.. (f dv ^_ dx dx2 du + 2C12- ( du ( dv U = dx (dy dx a c [ dw V (dw V [ dw' + p19\ -T-\ + p1,10\ -T~\ + p1,15 w3 + ^1,16 I -T~\ + + D11 + 2B1 dx dy dx dx dx2 dx (dy R J (dy y 1 dv , dw V [ dw V ( dw + p1,17 \ \ + p1,22w4 + p1,23 \ \ + p1,24 \ 2 d w dv w | du — + — I + 2C12 — dy RJ dx + B2 dy dx dy dy R dy dy2 (7) dw [ dw V [ dw V — = p21w1 + p— ^ \ + p23 \ -Г- \ + p28w2 + dx (dx J (dy 2 d w , dv w + 12 [~dy + ~R dx dx dw V ( dw V [ dw + p29 \ \ + p2,10 \ \ + p2,15 + p2,16 d2 w )( 3v y 2 d w ^V x _ dx dv 1 dv dx dy dx 12 dy R dy dy2 dx' dw V [ dw V (dw + p2,17 \ \ + p2,22w4 + p2,23 \ \ + p2,24 \ 2 d w 1 dv y + 2C dy dx 22 dy R dy dy2 dy R dw [ dw V [ dw V — = Pэ\R\ + p32 \ — \ + p33 \ — \ + p38w2 + dy (dx J (dy [ dw V [ dw V [ dw V + p39 \ ~~T \ + p3,10 \ "T” \ + p3,15w3 + p3,16 \ ~r~ \ + (12) ( dx L ( dy L ( dx J„ 2xy 3x2 y 3x2 2x 2y2 y3 y p29 = ^-—+—;---; p2i0 = ——+^—7+— ab a b a a , ab ab a 6x2 y 6xy y 3 y2 2 y3 p215 =—+-f--—+^-г —^; , a b a b ab ab b a 2x2 x3 x a3b 2 3 xx 2 ab ab a a _ xy 3x y 3xy 2x y 2xy p1,15 = T+ ' +" " a3b b3 a a3b b a ab ab2 [ dw V [ dw V [ dw V + p3,17 \ T~ \ + p3,22w4 + p3,23 \ T~ \ + p3,24 \ ~T \ ; (dy J3 (dx J4 (dy J4 Vx = p44 Vx1 + p4,11Vx2 + p4,18Vx3 + p4,25Vx4 ; Vy = p55Vy1 + p5,12Vy2 + p5,19Vy3 + p5,26Vy4 ; u = p66u1 + p6,13u2 + p6,20u3 + p6,27u4 ; v = p77v1 + p7,14v2 + p7,21v3 + p7,28v4 , где ptj - компоненты матрицы преобразования P размерностью 7x28 со следующими ненулевыми компонентами: , xy 3x2 y 3xy2 2 x3 y p11 = 1 —— + —t^- + xy 3x2y 3xy2 2x3y 2xy3 3x2 2x3 p18 =-T + —u---TT + ^T~ + TT~ + ~---T; ab a b ab a b b a a a 2 3 3 2 ,,2 3 x y x y x x 2xy xy xy p19 =4-—tt+т—; p1,10 =—^-+^+--r- 2 3 2 3 x y x y xy xy a,16=—г+~£; p1,17=-u+^2; ab a b ab ab xy 3x2 y 3xy2 2 x3 y 2 xy3 3y2 2y3 p1,22 = ~---;---- + —^^^^ +---• -.2 3 3 <■, 2 2x y x y xy x 2x p12 =—-—T- —-+—7--+x; ab a b b a a = 2 xy2 xy3 2 y2 + xy y3 + p13 =—г---и T~+--TT + y; ab ab b a b 2xy3 - 3^ - ^ 2^ 2y3 ; b3a a2 b2 a3 b3 ’ ab a 2b ab2 ab a b ab ab a b ab ab a b a = 2x2y + x3y + xy ; = 2y2 y3 y _ p1,23 = , + 2, + , ; p23 = ab a2b b 2xy 3x2 y y2 y3 p2,16 = 7”+ T~; p2,17 = 7+ 72; ab a b ab ab 6 x2 y 6 xy y 3 y2 2 y3 p., =--^------• 2,22 a3b a 2b ab ab2 b3a 4 xy 3x2 y y y2 y3 p223 =—+—; p2 24 ---=^r; , ab a b b , ab ab 6xy2 x 3x2 2x3 6y 6y2 6 xy p31 = b3a ab + a2b a3b b2 + b3 + ab2 ; p32 = 2 xy 3xy p3,16 = T + ~YT; p3,17 = T+ TT'; ab a b ab ab = 6 xy2 x 3x2 2x3 6y 6y2 6 xy ^ъ’22 b}a ab a2b a}b b2 b ab2 ’ 2 x2 x3 x. ^3 23 =---+--i--+--; , ab abb ab a b b = 1+ixy - 3 xy2 + 3/ - x - ; p33 ab ab2 b2 a b ’ 6xy2 x 3x2 2x3 6xy p38 —+~r TT72; b a ab a b a b ab p39 = , 2 ; p3,10 = ' , ' 2 ' ab a b ab ab a 6 xy2 x 3 x2 2x3 6xy p3 15 = —3----+—2---3—+—2; 3,15 b a ab a b a3b ab2 = 2 xy 3xy2 + 3 y2 2 y = 1 + xy x y p3,24 = , 2 + 2 ,; p44 = 1 + , , ; ab ab b b ab a b xy x xy xy y p4,11 = T^_ ; p4,18 =~T ; p4,25 = 7 + T ; ab a ab ab b , xy x y xy x p55 = 1 T T ; p5,12 = T^_ ; ab a b ab a 4 xy 3xy x 2 3 xx 2 3 2 3 p ^_^_+ у_. (13) p1,24 _ . , 2 l , 2 ; ab ab b b 6x2 y 6xy y 3 y2 2 y3 6x 6x2 p21=-+ab- ab+ab2- ^" a2+~; 4 xy 3 x2 y 3x2 4x y p— =1 +—,---TT+~2----T; ab a b a a b 6 x2 y 6 xy y 3 y2 2 y3 6 x 6x2 p28 = ^+ * "^T + Тз + ~--Г; a b a b ab ab baa a _xy,y. _ i , xy x y p5,19 = , ; p5,26 = , ^ , ; p66 = 1 + , ab ab b ab a b xy x xy xy y p6,13 = T+ ; p6,20 T ; p6,27 = T^T ; ab a ab ab b xy x y xy x p77 = 1 T T; p7,14 = T^_ ; ab a b ab a _ x^ . _ xy y p7,21 = , ; p7,28 = , ^ , . ab ab b (15) (16) 2EH t 1 -ц 11 - B22 - Ен (H3 - h3) 12(1 -Ц2н) , Он (H3 - h3) (2i) a2b ; 14; ab 14 ab a2b ab2 ab 2 Подставляя разложения (12) в выражения потенциальной энергии деформации (7) и потенциала внешних сил (8) и минимизируя их сумму (9) по компонентам вектора узловых кинематических параметров (10), (11), получим систему алгебраических уравнений равновесия Ke 6е = Fe, (14) где Ke - матрица жесткости конечного элемента цилиндрической оболочки; Fe - вектор эквивалентных узловых сил. Вектор эквивалентных узловых сил удобно представить в блочном виде: Здесь p - величина давления. В общем случае матрица жесткости конечного элемента цилиндрической оболочки Ke, размерность которой 28x28, имеет плотное заполнение, т. е. в ней немного нулевых элементов. Для ее вывода разработана специальная программа [2]. Вывод разрешающих уравнений при реализации вариационного подхода МКЭ для оболочек предполагает получение выражения полной потенциальной энергии для всего ансамбля конечных элементов модели: F -{ /2 /3 0 0 0 0}Т; F2 - // /9 h10 0 0 0 0}Т; F3 - {./15 f16 f\7 0 0 0 0}Т; F4 - {j~22 /23 /24 0 0 0 0}Т; ab - ■; h - a 2b p- 24 ; /3- ab 2 = p— 24 ■; / ab )-• 4 /9 =- P /ю = P /15 = P 24 a 2b fn - p 24 ; Fe = (17) 24 ab I-1 где подвекторы F], F2, F3, F4, состоящие из 7 компонентов каждый (в случае нагрузки равномерным давлением) имеют вид где N - число элементов в модели; U - потенциальная энергия деформации i-го элемента; П - потенциал внешних сил i-го элемента. Лб = - p—; /17 = - p—; f22 = p—; Минимизация (17) по узловым неизвестным каждого узла системы приведет к глобальным уравнениям равновесия К А = Fs, (18) где К - глобальная матрица жесткости пластины; Fe - глобальный вектор эквивалентных узловых сил всей пластины; А - глобальный вектор узловых неизвестных: А = { Й! 62 ... б; ... 6N}r. (19) Следует иметь в виду, что при решении СЛАУ (18) необходимо учесть граничные условия. Одной из самых распространенных конструкций, в расчетах которых необходимо учитывать трансвер-сальный сдвиг, является трехслойная оболочка с податливым заполнителем. Будем полагать, что жесткость несущих слоев существенно выше жесткости промежуточного слоя. В этом случае можно допустить, что параметры жесткости оболочки B, Си D обеспечены несущими слоями. Будем полагать материалы несущих слоев и заполнителя изотропными. При совпадении начальной поверхности со срединной и при симметричной структуре слоистого пакета смешанные жесткости (С) равны нулю. Остальные параметры жесткости вычисляются по формулам где H - полная толщина пакета; h - толщина слоя заполнителя; t - толщина каждого из несущих слоев; Eн - приведенный модуль упругости материала несущих слоев; G„ и G3 - модули сдвига материалов несущих слоев и заполнителя соответственно. Как известно, в последнее время получила распространение конечно-элементная модель податливых при трансверсальном сдвиге оболочек, в которой в качестве основных узловых переменных фигурируют перемещения точек начальной поверхности и углы поворота сечения (0Х и 0y). Эта модель с пятью кинематическими параметрами в узле обладает определенным недостатком, который называется эффектом сдвигового запирания (ЭСЗ). Покажем, что наша КЭ модель с семью кинематическими параметрами в каждом узле лишена данного недостатка. Для этого проведем следующий численный эксперимент. Выполним две серии расчетов для оболочек различной толщины. В первой серии будет задействована рассматриваемая в настоящей работе модель с 7 узловыми кинематическими параметрами. ^ - ^ - 1Т+Т G,, + G, D12 - D21 - Н-н D11, D33 - B33 - 2GJ, D11 - D22 - B12 - B21 - Н- н B11, 12(1 -ц2) (22) (24) 2 I я 2 А2 гл д w +D дх2 1 дv д2 w А R ду ду2 дх2 дv w ду + R ди ду + ду дх + B2 (26) dxdy. + B3 +D3 Ah 12 12 Djj - D22 - ATih A33h 12 12 12 D12 - D21 - (27) Размерность матрицы жесткости элемента сократится до 20. Вектор эквивалентных узловых сил сохранит блочный вид (15), но подвекторы (16) будут иметь размерность 5 с теми же значениями ненулевых компонентов. Рассмотрим цилиндрическую панель, с защемленными торцевыми (искривленными) кромками и свободными прямыми (см. рисунок), нагруженную равномерным давлением на внешней поверхности. Радиус кривизны панели - 5 м. Размеры в плане 1 м на 1 м. Будем полагать оболочку изготовленной из однонаправ- Во второй серии будем использовать конечно-элементную модель, основывающуюся на классической теории оболочек Кирхгофа. В этой классической модели в каждом узле присутствует 5 кинематических переменных: перемещения и и v, прогиб w и два изгибных угла поворота нормали (dw/dx и dw/dy). Таким образом, вектор узловых неизвестных 4-узлового элемента оболочки в данном случае будет иметь вид (10) со следующим подвектором узловых параметров в каждом узле: Процедура вывода разрешающих уравнений равновесия для одного конечного элемента точно такая же, что описана выше для модели, учитывающей трансверсальный сдвиг. Только во всех выражениях необходимо положить уХ = уу = 0 и, следовательно, 0x = -dw/dx и 0У = -dw/dy. Выражение энергии деформации при этом будет иметь вид U - -2 jjjB1.fl' 1 ду - 2_ R дх дхду д2w А2 (23) dw дУ „ _ ди I dv w . + 12 "дХ [ду+ ~r 1+ 12 ду R ду ду2 22 1 ду д w ' + D2 2 д w у A22 - Е sin4 ф + Е2 cos4 ф + + 2 (е.ц12 + 2G12 j sin2 ф cos2 ф; A12 - A21 - E1M-12 + + (e. + E2 - 2 (е.ц12 + 2G12) sin2 ф cos2 ф; A33 - (E. + E2 - 2 Е.ц12) sin2 фcos2 ф + G12 cos2 2ф, где E. и E2 - приведенные модули упругости, вычисляемые по формулам E1(2) - E1(2)/(1 -М-12М21) , (25) а коэффициент Пуассона ц21 определяется из условия симметрии упругих постоянных: М-12 Е1 ленного композита с углами укладки ф = ±45°. Зададим следующие механические свойства (углепластик): - модуль упругости вдоль волокон Е1 = 180 ГПа; - модуль упругости поперек волокон Е2 = 6,2 ГПа; - модуль сдвига G12 = 5 ГПа; - коэффициент Пуассона д12 = 0,007; - плотность р = 1500 кг/м3. Упругие параметры пластины, которую будем считать однослойной и условно однородной, вычислим по формулам A11 - Е cos4 ф + Е2 sin4 ф + + 2 (е.ц12 + 2G12) sin2 ф cos2 ф; Модули трансверсального сдвига примем равными G12. Коэффициенты жесткости пластины вычислим по формулам B11 - A11 h ; B22 - A22h ; B21 - A21 h ; B33 - A33h ; D33 - Н-21 - " A„ h3 Ah E2 Kx - Ку - G.2h . Композитная панель и картина деформирования Таблица 1 Прогибы однородной композитной панели Толщина оболочки, h, мм Прогиб в центральном узле, мм Расхождение, % Давление, Па Классическая модель Модель, учитывающая сдвиг 50 4,593 5,070 9,4 1000000 40 8,579 9,144 6,18 1000000 30 18,566 19,219 3,40 1000000 20 4,993 5,055 1,23 100000 10 1,7817 1,7811 0,034 10000 5 2,8689 2,8619 0,24 10000 1 -0,522263 -0,522286 0,0044 100 Таблица 2 Прогибы трехслойной панели Толщина слоя заполнителя, h, мм Прогиб в центральном узле, мм Расхождение, % Давление, Па Классическая модель Модель, учитывающая сдвиг 50 4,718 10,717 55,98 100000 40 7,302 14,673 50,24 100000 30 6,369 11,119 42,72 50000 20 2,778 4,117 32,52 10000 10 4,990 6,101 18,21 5000 5 1,610 1,765 8,78 500 3 0,675 0,708 4,66 100 1 2,02 2,045 1,22 100 0,5 1,4358 1,4444 0,595 50 0,1 0,39015 0,39065 0,128 10 Определим прогибы в центральном узле, вычисляемые по классической (Кирхгофа-Лява) и сдвиговой моделям для различных значений толщины панели. Результаты расчетов представим таблично (табл. 1). Анализ этих данных позволяет заключить, что по мере увеличения толщины оболочки классическое решение дает существенно заниженный прогиб. С другой стороны, по мере уменьшения толщины пакета решение по сдвиговой модели приближается к классическому. Еще более рельефно отмеченные эффекты проявляются в расчетах по классической и сдвиговой моделям, выполненным для трехслойной панели. Длина панели (по образующей - 1 м, ширина (в окружном направлении) - 0,2 м. Радиус кривизны - 5 м. Толщина стальных несущих слоев - 1 мм, толщина слоя заполнителя варьируется. Коэффициенты жесткости вычисляются по формулам (21). С помощью специально разработанной программы [3] выполнена серия статических расчетов панели при нагрузке давлением. При этом оцениваются прогибы в центральном узле, полученные для трехслойной оболочки с различной толщиной слоя заполнителя. Параллельно выполняются расчеты по двух КЭ моделям: сдвиговой и без учета сдвига (Кирхгофа-Лява). Результаты представлены таблично (табл. 2). Данные таблицы показывают, что по мере уменьшения толщины слоя заполнителя результаты решения с использованием сдвиговой модели приближаются к результатам, полученным на основе классического решения, что свидетельствует об отсутствии ЭСЗ в разработанной КЭ модели. Эти же результаты лишний раз подчеркивают необходимость применения сдвиговой модели для расчета конструкций, весьма податливых при трансверсальном сдвиге. Разработана конечно-элементная модель цилиндрической оболочки, податливой при трансверсальном сдвиге. Разрешающие уравнения теории МКЭ получены вариационным способом. Вектор узловых неизвестных содержит семь независимых кинематических параметров, включая осредненные по высоте сечения углы трансверсального сдвига. В результате численного эксперимента на примерах расчетов монолитной композитной и трехслойной оболочек показано, что в КЭ модели с полным набором узловых кинематических параметров отсутствует эффект сдвигового запирания.