Применение метода конечных элементов к расчету на прочность круговых цилиндрических оболочек с вырезами Текст научной статьи по специальности «Физика»

Том XVI

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

19 85

№ 3

УДК 539.4

ПРИМЕНЕНИЕ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ К РАСЧЕТУ НА ПРОЧНОСТЬ КРУГОВЫХ ЦИЛИНДРИЧЕСКИХ ОБОЛОЧЕК С ВЫРЕЗАМИ

В. В. Кабанов, Л. П. Железное

Разработан алгоритм определения напряженно-деформированного состояния круговых цилиндрических оболочек и пластин переменной толщины, ослабленных вырезами произвольной формы и величины при произвольной нагрузке. Алгоритм основан на методе конечных элементов в перемещениях с использованием плоских и криволинейных произвольных четырехугольных конечных элементов естественной кривизны. Перемещения точек конечного элемента аппроксимируются степенными полиномами. Система линейных алгебраических уравнений решается методом разложения матрицы на две треугольные. Алгоритм реализован программой для ЭВМ БЭСМ-6. В примерах рассмотрены пластины с круглым отверстием при растяжении, оболочка с круглым и прямоугольным вырезом при растяжении и наддуве. Результаты расчетов сравнивались с результатами других авторов.

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

1. Рассмотрим круговую цилиндрическую оболочку радиуса Я, длиной Ь, толщиной к, нагруженную системой поверхностных распределенных нагрузок <7г(л;, у), локальных сил Рн и моментов Мц, погон-

ных контурных сил Рм(у) и моментов Mki(y) (индексы t=l, 2, 3 соответствуют направлениям осей х, у, z) (рис. 1).

Оболочку разобьем в продольном направлении на т частей, в поперечном— на п частей. Таким образом, оболочка представляется набором тХп криволинейных четырехугольных конечных элементов (см. рис. 1).

Перемещения точек конечного элемента представим в виде U = a^Y] + а2$ + <xa7j + ос4, V — a5Sr¡ + а6£ + а,т] -f а8,

w = а9^3 + a10¡;V + ап?ц + «ИЕ8 + «1B$V + «i«5V + «Л +

-f aí6?2 + a17^Tf¡3 + ZlaW + + «20^ + «217!3 + "iírf + «237! + «24- (1)

Рис. 1 У,

В матричной форме (1) имеет вид

и = Ря, (2)

где Р — матрица связи порядка (3X24), и={и, v, w} — вектор (столбец) перемещений точек элемента et = {aIf a2, . . ., a24} — вектор неизвестных коэффициентов полиномов.

Выразим at через узловые перемещения (?, tj в индексах означают дифференцирование по %, ч\) u = {ui, vb wlt wn, wni, wini,

Up . . ., W^j, Uk, . . ., W^k, Ию ■ • •• Щ-цп}-

Из (1) имеем

a = Вя. (3)

Элементы матрицы В имеют вид

bii = P\jy b2j—p2j, bs j = ps j, ¿>4/= (p3 b5j = (paJ)n, У bü]=-(Ps])^ при S=—a, Tj=— b\

bTj — Pij, bSj=pljt b9J=p3J, bi0J = (p3 bnj = (p3j)rt,

bi2j = (Pv)tn при s = a, ti=—b] ¿is j=Pij, bn}=p2j, bt&J=-p9J, bUJ~(ptJ)u b17j = (p3J)v (4)

¿is j~(j>Bj)in ПРИ \ — а, 1) = b; ¿19 j — Pi j, ¿20 / = p2 j> ¿21 J = Pa я ¿22 j= (Рз j)b ¿23 j = (/?3 b2i j = (P3jhn При 6=i- Л, = j — I, ... , 24,

= при £ = a, -»J — ft, a = 1, ¿> = 1; ]

здесь S, tj в индексах означают дифференцирование.

Отсюда

« = в, и — Р\ и, РХ = РВ-К

(5)

Здесь В — матрица порядка (24X24).

Кинематические и статические соотношения для конечного элемента возьмем в виде

е = (6)

где е={е„ е2, е„ Хи Хг. ХвЬ Т=[Ти Т2, Т„ Ми М2, М$] - векторы деформаций и усилий конечного элемента

( ), о ( ), о О (),()* О О к2 0 -( ),

в, в2 о о о о В2 Вх ООО о о о в3 о о о

0 О О />! 02 О О 0 0 й2 О

о о о о о я3

о о

М ), ( ),

_( 2 ( )ху

Ек

, В2 = чВи

1

= ^ (1 — V) Ви

ЕЬ?

12 (1-^)

1

/)3=-—(1 3 2

(7)

Запишем связь координат х, у с криволинейными координатами т] (см. рис. 1):

х = х(5, у = у(6, т)), (8)

где согласно [7]

■* (5. ч) = /1 (6. ч) + Л (5, ч) + /з (5, ч) + /4 (5, ч) .У (5, Ч)=/1 Ч) У1 +/а (5, Ч) Уг + /з (Е, Ч) Уз +/* (6, ч)у4,

4 у

Л (6, ч) = + О + ч). Л & ч) - у О - &) 0 + ч),

(9)

лг, — узловые координаты.

Производные функций /($, ч) по координатам 5, ч представим формулами:

/' = о/7', Г = нр",

(10)

где

/'-{ДЛЬ /«„/„}.

^ {•/*> /уЬ ^ = \fxxt Уху /уу}>

Х-Ц У-Г]

я=

х* у^

.(П)

здесь ч в индексах означают дифференцирование по ч-94

Из (10) находим

Элементы матриц б-1 и Н~1 имеют вид

gu=Уrlí Д, йа=-Уб/Д, =*-*ч/д» ЙГ22 =- -«е/Д,

'21 '

Й23=-^У£/Дь А,1=лг*/Д1, А82--

Л33 = -*:|/Д1. Д = Уч —Уб, Д1 = Д2-Производные л и у по и| получим из (8) с учетом (9)

хе «= а^, + б^з + с^з + .Уе = «1У1 4- + ^уз + •"■г] — ^2*2 -)- ¿2-^2 -(- Сг-^З

^ = а2У1 + &2.У2 + ^Уз + ^2у4,

где

4 4 4

1

(1 + %

Аг =

1*п( )е + ^0,1 о о

О [й()£ + &(),] Ы>

Ып ( + £22 ( ),,] [^11 ( )5 + ^12 ( )ч] О

О о [-А„( )К-А„( )щ]

О ¿2 [ви ( )£ + ёп ( ),] [- Л„ ( - А,| ( )е,, - Азз ( ]

О 2*2 ( )Е + gn ( )ч] 2 [- Л21 ( - й22 ( )5ч - Л23 ( )„]

(13)

(14)

(15)

(16)

Используя (12), получаем матрицу А как функцию от -ц.

Запишем энергию деформаций и работу внешних сил а для конечного элемента:

*/= ++(17)

5 .$ I

Здесь Я={Я\,Яг>Яъ) — вектор внешней поверхностной нагрузки,

Я к — А 1> Р*Ъ. ,, Мк Л1к 3}, /?г = {Рп , Рп, Рш М1,, М12,

М13\ — векторы контурных и локальных сил, действующие на конечный элемент, иА = {ы, V, да, чи^}, и1={ии vl, тц, ии^и

Введя новые переменные, получим

с1з=:\<1е1[0\\с&с1г1, (11 = \&е\\а]\<П. (18)

Выражение полной потенциальной энергии конечного элемен-

та:

г 1

П, = 1У—| | 7Т(£, 7)) «(5, ц) |(1е1[<?]|

-1-1

11 _ __ _

(19)

-1-1

С учетом (5), (6) имеем

Л.= & (В- ')т Л Р'АЮАР | бе! [0] \й\йгхВ^и —

— Ц* <г Р\6е\[в\\<&ащВ-*Ъ- | Я1Рк\йе1[б]\(йВ-1й-

— Жьи^и^Ки — Ой— <2ки — (?,» ,

(20)

где

0=

I сИ [С] ] (Ка-ц В-\ = Ц Щрк I <1е1 [О] I \

(21)

(Рм) 1} = (Л у). (/7а>2 ] = р2 я (/>*)8 / = />з я 7 = 1, - . .,24, (РкЪ } = [¿Гн (Л > + £12 (/>з Л, (/>*)в } — \ёг\ {Рз М + 2м (Л М (Рк)б у = [МЛ ;)« + А22 (Ра + Л2з (/>8

Узловые перемещения и в системе координат I, -г\ связаны с узловыми перемещениями и в системе координат х, у зависимостью

и

Си,

(22)

где

X, 0 О О

О О О

0 0^0

0 0 0 1

Л,-

1 о о О 1 о О 0 1

Х-цУ^

О О

О О

Xт

(23)

х\, — производные х, у по \ в точке г конечного элемента. Используя (23), из (20) получаем

Пг = &К*и -<?*»- - <?* й, (24)

где

X* — Ст КС, О ос, я;=якс, <)*, = <}£. (25)

Суммируя потенциальные энергии отдельных элементов, получаем потенциальную энергию оболочки. Варьируя потенциальную энергию

ло узловым перемещениям согласно пришшиу возможных перемещений (611 = 0). получим, как и в [7], систему линейных алгебраических уравнений для определения узловых перемещений

Ки'^Я. (26)

Здесь К — матрица упругой жесткости оболочки и'; ф—векторы узловых перемещений и обобщенных узловых сил оболочки.

Матрица К имеет ленточную структуру и получается суммированием элементов матрицы К* с использованием матрицы индексов согласно [8]. Вектор ф получается суммированием векторов (}'= = <?* -г С}1 (2*1- Система (26) решается методом разложения матрицы на две треугольные [9|.

Решив систему (26), найдем вектор и' и по формуле (5) все компоненты напряженно-деформированного состояния оболочки. Полученное решение можно использовать и для пластин, если матрицу А записать для плоского конечного элемента

( ), 0 ( ), 0 0 С)

Аг = 0 ()у (), 0 0 0 (27)

0 0 0 -( )уу -( ),у

2. С целью проверки точности и эффективности разработанного алгоритма были рассмотрены тестовые задачи: длинная пластина н оболочка с круговым отверстием при одноосном растяжении, оболочка с круговым отверстием при внутреннем давлении. Пластина имела ширину 26=10/'. оболочка имела длину Ь — 4,2 м, радиус /?=1,34 м и толщину /г = 0,02 м. Радиус отверстия г = 0,5 м. Сравнение результатов расчета (числитель) с решением в рядах [1] (знаменатель) показано в таблице.

0 1 0 | 11° 23® 34° 46° 55° 67° 78° 90°

кв 2,3 2,1 1,9 1,2 0,6 -......0,5 — 1,4 — 1,8 —2,0

2,1 2,0 1,8 1,1 0,5 -0,4 — 1,0 — 1,7 -1,9

К —0,3 -0,5 -0,7 -0,7 0 1,3 3,6 4,9 5,4

—0,4 -0,6 -0,8 —0,8 0 1,5 3,8 5,2 5,7

Щ 11 12 15 18 17 11 4 -3 _5

13 15 17 20 ' 18 12 4 -3 —5

ч -13 —12 -9 —4 3 9 14 17 18

-13 — 12 -10 _4 3 9 14 17 18

Здесь 0 угол (полярная координата в градусах), /гс = а0с/о) /гс = о^с/з,, = рЩЬ. — коэффициенты концентрации напряжений ос в срединной поверхности оболочки у контура отверстий при растяжении и внутреннем давлении; йв = з0в/о, /г| — о?в[вр — отношения

7— '-Ученые записки Ц.АГИ» Л1 3 97

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

На рис. 2 показаны зависимости коэффициентов концентрации к — = сг1/сг, 1г1С = 0{1а осевых напряжений в срединной поверхности пластины и оболочки по среднему сечению и значение коэффициента концентрации осевых напряжений к^=х>\1а у внутренней поверхности оболочки. Координата х = у/Ь — для пластинки, х = 2у/лК — для оболочки, начало

Рис. 2

координаты помещено у контура отверстия. Пунктиром показано аналитическое решение для пластины [10], кружками — расчет по предлагаемому алгоритму. Сплошные линии соответствуют-числу элементов /?2=13, п= 14, звездочки (к1с) и ромбики (кщ) >1 = 7, л= 14 (рассматривалась одна четвертая часть оболочки и пластины).

Из полученных результатов видно достаточно хорошее соответствие сравниваемых величин.

На рис. 3 и 4 для различных 6 = /?1//г показаны значения /гс, и кв в оболочке длиной 1 = 4,2 м, радиуса ./? = 1,34 м, толщиной /г = 0,02 м, с квадратным вырезом шириной 2г0= 1,0 м и скругленными углами радиуса г = 0,2 м при растяжении. Вырезы имеют подкрепления прямоугольного поперечного сечения толщиной /г 1 н шириной 6 = 0,05, м.

Видно, что максимальная концентрация напряжений в ;срединной поверхности оболочки возникает при 0 = 55°, а изгибных напряжений —

при 6 = 0 и достигает соответственно значений 4,5 и й,1 при 6='1, т. е. когда края не. подкреплены. С увеличением б максимальная концентрация напряжений в срединной поверхности оболочки уменьшается,

Рис. 3

Рис. 4.

а максимальная концентрация изгибных напряжений увеличивается и при б>2 практически .йс зависит от 6.

.,. .Время решения задачи при т— 13 и п= 14 с числом неизвестных М= 1260 и при полуширине ленты матрицы жесткости оболочки В = 96 составило 28 мин машинного времени ЭВМ БЭСМ-6. Решение аналогичной задачи в [6] при той же точности требует 533 конечных элементов с общим числом неизвестных порядка 5500. Таким образом, разработанный алгоритм позволяет эффективно определять напряженно-деформированное состояние оболочек и пластин с вырезами произвольной формы. При этом вид нагрузки не играет существенной роли.

ЛИТЕРАТУРА

1. Г у з ь А. И.'Цилиндрические оболочки, ослабленные отверстиями.— Киев: Наукова Думка, 1977.

2. Utku S., arid М е 1 о s h R. I. Behavior of triangular shell element stiffnes matrices associs ted wiht polyhedral deflection distributions. — A1AA Paper, AIAA 5-th Aerospace Sci Meet., N. Y„ Jan. 1967, N 67—114.

3. С 1 о u g h R. W. and Janson C. P. A finite element approximation for the analusis of thin shells. — Int. J. Solids Struct., 4, 1968.

4. Gallagher R. H., and Yang H. T. Y. Elastic instability predictions for doubly curved shells. — Proc. of the 2 nd Conf. on Matrix Methods in Struct. Mech. AFFDL — TR—68—>150, Air Force Inst, of Technol., Oct. 1968.

5. Bogner F. K., Fox R. D. and S с limit L. A. A Cylindrical shell discrete element.—AIAA J., 5, 1965.

6. Кей С. В., Бейсинджер 3. Е. Расчет тонких оболочек на основе метода конечных элементов. — В сб.: Расчет упругих конструкций с использованием ЭВМ. — Л.: Судостроение, 1974.

7. С е г е р л и н д Л. Применение метода конечных элементов, — М.: Мир, 1979.

8. П о с т н о в В. А., X а р х у р и и И. Я. Метод конечных элементов з расчетах судовых конструкций. — Л.: Судостроение, 1974.

9. Райнт У и л к и н с о и. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра. — М.: Машиностроение, 1976.

10. Савин Г. Н. Концентрация напряжений около отверстии.— М.—Л.: Гостехиздат, 1951.

Рукопись поступила 25/X 1983