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

УДК 539.3

А. Д. Матвеев, А. Н. Гришанов

Одно- и двухсеточные криволинейные элементы трехмерных цилиндрических панелей и оболочек*

A. D. Matveev, A. N. Grishanov

Single and Double Grid Curvilinear Elements of ThreeDimensional Cylindrical Panels and Shells

Изложены процедуры построения однородных односеточных и композитных двухсеточных криволинейных конечных элементов для расчета трехмерных упругих цилиндрических панелей и оболочек. Матрицы жесткости и узловые усилия криволинейных конечных элементов определяются в локальных декартовых системах координат, а системы уравнений метода конечных элементов для дискретных моделей оболочек и панелей — в глобальных декартовых системах координат. Связь между локальными и глобальными системами координат осуществляется с помощью матриц вращений, которые определяются только для векторов узловых перемещений криволинейных конечных элементов. При построении основных соотношений криволинейных конечных элементов для аппроксимации перемещений используются известные полиномы и уравнения трехмерной задачи теории упругости, записанные в локальных декартовых системах координат данных конечных элементов. Достоинства предлагаемых элементов заключаются в следующем.

Одно- и двухсеточные криволинейные элементы описывают трехмерное напряженное состояние цилиндрических оболочек и панелей, учитывают их сложный характер закрепления и нагружения. Двухсеточные элементы учитывают неоднородную и микронеоднородную структуру панелей и оболочек и порождают дискретные модели малой размерности. Ключевые слова композиты, упругость, цилиндрические оболочки, панели, одно- и двухсеточные криволинейные конечные элементы.

БОТ 10.14258/1гуа8и(2014)1.1-19

Construction of homogeneous single-grid and composite double-grid curvilinear finite elements for calculating parameters of three-dimensional cylindrical panels and shells is presented. In this paper, stiffness matrix and nodal forces of curvilinear finite elements are determined in local Cartesian reference systems. The set of equations of finite element method for discrete models of shells and panels utilizes global Cartesian reference system. Relation between local and global coordinate systems is carried out using the rotation matrix defined only for vectors of curvilinear finite elements nodal displacements. Relations for finite element characteristics of cuboid for further displacement approximation are based on known polynomials and equations of threedimensional elasticity problem in local Cartesian reference system.

The proposed single-grid and double-grid curvilinear elements provide three-dimensional stress state description while taking into account complexity of attaching and loading. Double grid curvilinear elements consider the heterogeneous and micro-inhomogeneous structure of panels and shells and serve as the basis for discrete models of small dimension.

Key words: composites, elasticity, cylindrical shells, panels,

single — grid and multi — grid curvilinear finite elements.

Введение. В данной работе для расчета однородных и композитных цилиндрических панелей и оболочек по методу конечных элементов с применением мелких разбиений предложены криволинейные одно- и двухсеточные конечные элементы 1, 2 и 3-го порядков. Двухсеточные конечные элементы формы прямоугольного параллелепипеда для анализа деформирования упругих тел неоднородной струк-

туры изучены в работах [1, с. 3; 2, с. 161]. Здесь рассматриваются следующие конечные элементы. Криволинейный односеточный конечный элемент V 1-го порядка размерами кех хкеу хкег показан на рисунке 1, где ае — угол раствора; О1х1 у1г1 — локальная декартовая система координат и г1О1 у1 — плоскость симметрии конечного элемента V ; cd — ось цилиндрической оболочки; Щ Щ — радиусы

* Работа выполнена при финансовой поддержке РФФИ (проект 11-01-00053).

нижней и верхней поверхностей конечных элементов V; № — толщина; Н‘ — длина конечного элемента V; Н‘х = аеЩ , узлы отмечены точками. Прямоугольники размерами № хк‘у есть боковые грани, криволинейные прямоугольники со сторонами Не, аеКё, аеКё — торцевые грани конечного элемента V. Форма конечного элемента V есть прямая призма высотой Н‘ , основанием которой является криволинейный прямоугольник. Однородные односеточные криволинейные элементы 2-го и 3-го порядков и композитные двухсеточ-

г1|

Рис. 1. Конечный элемент Уе 1-го порядка

ные конечные элементы 1, 2 и 3-го порядков имеют такую же геометрическую форму, как и конечного элемента V (рис. 1). Криволинейный односеточный конечный элемент 2-го порядка имеет 20 узлов (рис. 2а), конечный элемент 3-го порядка — 32 узла (рис. 26), узлы отмечены точками. При мелком разбиении углы раствора предлагаемых криволинейных конечных элементов малы (угол ае мал; рис. 1), поэтому формы криволинейных конечных элементов мало отличаются от форм прямоугольных параллелепипедов [3, с. 136; 4, с. 212].

б

Рис. 2. Конечные элементы 2-го и 3-го порядков

В связи с этим для аппроксимации перемещений криволинейных конечных элементов используем известные полиномы 1, 2 и 3-го порядков (которые применяются при построении основных соотношений

элементов формы прямоугольного параллелепипеда) и уравнения трехмерной задачи теории упругости, записанные в локальных декартовых системах координат конечных элементов (рис. 1-3).

Достоинства предлагаемых конечных элементов состоят в следующем. Криволинейные одно- и двухсеточные конечные элементы учитывают сложный характер крепления и нагружения панелей и оболочек, описывают трехмерное напряженное состояние, сколь угодно точно (при мелких разбиениях) представляют упругое поведение панелей и оболочек в случае действия локальных нагрузок. Двухсеточные конечные элементы учитывают неоднородную и микронеодно-родную структуру панелей и оболочек, образуют дискретные модели, размерности которых на несколько порядков меньше размерностей базовых моделей, напряжения могут быть определены в любом компоненте неоднородной структуры оболочки и панели. Реализация метода конечных элементов для двухсеточных дискретных моделей оболочек и панелей требует меньше памяти ЭВМ и временных затрат, чем для базовых моделей.

1. Односеточные криволинейные элементы. Рассмотрим односеточный однородный криволинейный конечный элемент V е 1-го порядка со сторонами кех хкеу хк (см. рис. 1). Перемещения, деформации и напряжения элемента V е удовлетворяют соотношениям Коши и закону Гука [3, с. 300], которые отвечают трехмерной задачи теории упругости [5, с. 40], т. е. во всей области для конечного элемента V е реализуется трехмерное напряженное состояние. Для аппроксимирующих функций перемещений ые, у е, w е элемента Уе используем полином 1-го порядка, записанный в локальной декартовой системе координат 01х1 у1г1 (рис. 1),

«е , V , ^ = а1 + а2 Х1 + а3 У1 + а4 *1 + а5 Х1У1 + аб *1Х1 + (1) + а7 *1У1 + а8 Х1У 1 ^

Обозначим: «' - перемещение узла г конечного элемента V в направлении оси 01 х1 ; х'1 , у'1, г'1 - координаты г -го узла конечного элемента V ; г = 1,...,8; номера узлов отмечены на рисунке 1. Коэффициенты а{ в (1) находим из условия ые (х1 , у1,г[) = «', г = 1,...,8, которое представим [С]а = и , где а = {а1,...,а8}т, и = {«1,...,«1 }г ; [С] — матрица размерности 8х8 ; Т — транспонирование. Отсюда следует: а = [С]-1 и . Используя найденный вектор а в (1), определяем функцию «е. Аналогично строим аппроксимирующие функции перемещений Уе , w е конечного элемента V е, которые запишем

8 8 8

«е = ^ N«1, у е = ^ N^1, м’е = ^ N^1, (2)

1=1 1=1 1=1

где «1 , V , wl1 — перемещения узла г ; Nj = Nj (х1, у1, г1) — функция формы г -го узла конечного элемента V ; г = 1,...,8.

Вектор функций перемещений ие = {«е , У е , }т ко-

нечного элемента V е представим

и = N К, (3)

где Ь\ — вектор неизвестных конечного элемента V е , имеющий в локальной декартовой системе координат 01 х1 у1г1 структуру

1 {«;,...,«8,у^1,...,у8,^!,...^^} , (4)

^е ] — матрица функций формы конечного элемента V е

Полную потенциальную энергию конечного элемента V е запишем в матричной форме

IV =1 Г ет а йV -Г UTFdV -Г ита йБ , (5)

е 20 veee Jve е е и Бе е е где ее = Ц, еЛ, е*, 7х,у,, 7уЛ, !хл }Т ;

а е = {ах , % , аг1 ,Тху , Ту1г1 , Тх1г1 }Т — векторы функций

деформаций и напряжений [3, с. 300], найденных в локальной декартовой системе координат 01х1 у1г1 конечного элемента V е; Б , це — векторы объемных и поверхностных сил конечного элемента V е; V е, Бе — область и поверхность конечного элемента V е.

Соотношения Коши и закон Гука для конечного элемента V е имеют вид

е е = В 6 , ае = [Ве ]е е, (6)

где [Ве ], [Бе ] — матрицы деформаций и модулей упругости конечного элемента V е.

Подставляя (3), (6) в (5), получаем

V61) =1 [ 6)т[В ]тБЩ6^-

2 ° е

б)N]тРеdV-£ (6])г ]т qeйБ . Из выполне-

ния дWe (ёге)/ 36] = 0 следует [К1 ] 6] = Ре‘, где [К1] = £ [Ве ]т[Бе ][Ве ]йУ , ре = X N]тре-/5 N]тqeйБ , здесь [К]] — матрица жесткости и Р1 — вектор узловых сил конечного элемента V е (которые определяются в локальной декартовой системе координат 01 х1у1^1 конечного элемента V е).

Процедуры построения односеточных криволинейных конечных элементов 2, 3-го порядков, имеющих такую же форму, как и конечный элемент Ve (рис. 1), аналогичны выше описанной.

2. Композитные криволинейные двухсеточные элементы. Рассмотрим трехмерные композитные криволинейные двухсеточные конечные элементы 1, 2, 3-го порядков двух типов.

2.1. Изложим процедуру построения двухсеточного конечного элемента 3-го порядка 1-го типа, который обозначим через Va (рис. 3), где 0хуг — локальная декартовая система координат двухсеточного конечного элемента Va; г0у — плоскость симметрии; ка — толщина; кау — длина двухсеточного конечного элемента Va. Считаем, что между компонентами неоднородной структуры двухсеточного конечного

элемента Уа связи идеальны. Функции перемещений, напряжений и деформаций компонентов двухсеточного конечного элемента Уа удовлетворяют закону Гука и соотношениям Коши, которые отвечают трехмерной задаче теории упругости, т. е. во всей области двухсеточный конечный элемент Уа реализуется трехмерное напряженное состояние. Область двухсеточного конечного элемента Уа представляем базовым разбиением Яа, которое состоит из однородных односеточных криволинейных элементов Уе 1-го порядка (рис. 1). Базовое разбиение Яа учитывает неоднородную структуру двухсеточного конечного элемента Уа и порождает мелкую трехмерную узловую сетку На. На мелкой сетке На определяем крупную сетку На. Узлы крупной сетки На на рисунке 3 отмечены точками — 32 узла. На сетке На строим аппроксимирующие функции иа, Уа, ша двухсеточного конечного элемента Уа, используя полином 3-го порядка (записанный в локальной декартовой системе координат Охуг; рис. 3)

иа, Уа, ша = а1 + а2 х + а3 у + а4 г + а5 ху + а6 гх + а7 гу + +а8 хуг + а9 х2 + а10 у2 + а11г2 + а12 х2 у + а13 х2 г + а14 у2 х + +а15 у 2г + а16 г2 х + а17 г2 у + а18 х3 + а19 у3 + а20г3 + а21х3 у + +а22 х3 г + а23 у3 х + а24 у3 г + а25г3 х + а26 г3 у + а27 х2 уг + +а28 у2 хг + а29 г2 ху + а30 х3 уг + а31 у3 хг + а32 г3 ху.

По процедуре, аналогичной процедуре п. 1, на сетке На для двухсеточного конечного элемента опре-

деляем аппроксимирующие функции перемещений и , V , ш вида

а > а* а.

= VN.u., v =VN.v., w =VN,1

/ j і і’ a / j і і’ a / j і

(7)

где и, V., ші — перемещения і -го узла крупной сетки На; N. = (х, у, 2) — функция формы і -го узла

сетки Н ; і = 1,...,32 .

Вектор ёа неизвестных крупной сетки На в локальной декартовой системе координат Охуг имеет структуру

Sa = {ul,...,u32, vl,...,v32,wl,...,w32} .

(8)

Пусть ось О1 у1 конечного элемента (располо-

женного в базовой модели двухсеточного конечного элемента Уа; рис. 1) параллельна оси Оу (рис. 3) и между осями О1 х1 и Ох угол равен ве, вектор Ре, де неизвестных конечных элементов V, в декартовой системе координат Охуг имеет структуру

с Г е ее ее е ->Т

6е = {u1,..., и8 ,v1,...,v8 ^..^ } ,

где uie , vie

(9)

конечного эле-

ш. — перемещения узла мента V , і = 1,...,8.

Между векторами ^ , Ре, имеем связь й,1 = [Те ]ёе

[3, с. 237], где [Те ] — матрица вращений [4, с. 64], размерности 24 х 24, в силу (4), (9) имеет вид

[Te ] =

[El] [E0] [E2]

[E0] [E] [E0]

-[E2] [E0] [E1]

(і0)

здесь [E0 ],..., [E2] — матрицы размерности 8x8; [E0] — нулевая и [E] — единичная матрицы; [EJ = cos в [E], [E2] = sin в [E].

Используя равенство 81 = [Г ]ёе, получим соотношения [3, с. 237]

[Ke ] = [Te ]T [K‘][Te ] , Pe = [Te ]T Pe1

(ii)

где [Ке ] — матрица жесткости и Ре — вектор узловых сил конечного элемента Уе, определяемые в локальной декартовой системе координат Охуг двухсеточного конечного элемента Уа (рис. 3).

Полную потенциальную энергию I а для базового разбиения Яа запишем

,1

Па =Е (2 ST [Ke ]Se - ST Pe ),

e=1 2

(і2)

где т — общее число конечного элемента Уе базового разбиения Яа.

С помощью аппроксимаций (7) вектор Ре, д е выражаем через вектор 6

6е = А ]6а, (13)

где [а; ] — прямоугольная матрица, размерности 24x96, е = 1,...,т .

Используя (13) в (12), из условия дПа (6а )/д6а = 0 получаем [Ка ] 6а = ¥а, где

т т

[Ка ] = £[ А! ]Т [Ке ][А; ], Ба =£[А; ]Т Ре , (14)

е=1 е=1

[Ка ] — матрица жесткости; Ба — вектор узловых сил двухсеточного конечного элемента Уа 1-го типа.

2.2. Кратко рассмотрим процедуру построения двухсеточного конечного элемента 3-го порядка 2-го типа, который обозначим через . Двухсеточный конечный элемент V имеет размеры, форму и неоднородную структуру, как и двухсеточный конечный элемент Уа (рис. 3). В данной процедуре используем мелкую На и крупную На сетки, базовое разбиение Яа и аппроксимирующие функции перемещений иа, va, ша двухсеточного конечного элемента Уа (см. п. 2.1), т. е. 6а — вектор узловых неизвестных двухсеточного конечного элемента V. На базовом разбиении Яа двухсеточного конечного элемента V с помощью метода конденсации [4, с. 248] строим суперэлемент . Систему уравнений для разбиения Яа представим

[ЛД +[ A12 ]Ss = B1 , [AJSl +[A22]Ss = B2

(і5)

где 61, 65 - векторы перемещений соответственно внутренних и граничных узлов сетки На разбиения Яа суперэлемента 05; [А11], [А22] ([А12], [А21]) — квадратные (прямоугольные) матрицы; В1, В2 — векторы узловых сил.

Из (15) для суперэлемента получаем систему уравнений [К5 ] 65 = Р , где [К5 ] — матрица жесткости; Р5, 65 — векторы узловых сил и неизвестных суперэлемента С ; [К] = [А22]-[А21][АП]-1[А12], Р5 = В2 — [А21 ][А11 ]-1 В1. Полную потенциальную энергию П5 суперэлемента запишем в виде

=1

=1

=1

математика и механика

П =1ST [K ]S - STP .

(і6)

Используя (7), между векторами 6,, 6а установим связь

6 = [ А ]6а , (17)

где [А, ] — прямоугольная матрица; 6а — вектор узловых неизвестных двухсеточного конечного элемента У,.

Подставляя (17) в (16), из условия дП, / д6а = 0 получаем [К, ] qa = Р , где

[Ks ] = [As ]T [Ks ][As ], Fs = [As ]T Ps

(і8)

сеточного конечного элемента Vs — 2Не. Модуль Юнга связующего материала равен 1, волокон — 10, коэффициент Пуассона — 0,3. На внешней поверхности панели в узлах і,}, к заданы силы qz = 0,01, где к = 19 , і = 46 + 9(а-1), } = 10 + 9(в-1), а = 1,...,5 , в = 1,...,9 . Ось О2у0 (рис. 5) параллельна оси Оу (см. рис. 3), угол между осями О2х0 и Ох равен в,, вектор ёа узловых неизвестных двухсеточного конечного элемента V, в глобальной декартовой системе координат О2 х0 у0 20 имеет вид

Л0 — Г„0 „0 ,,0 ,,0 „,0 „,0 \T

Sa {U1 ,..., U32 , V1 ,..., V32 , wl ,... W32 }

(і9)

[К, ] — матрица жесткости и Fs — вектор узловых сил двухсеточного конечного элемента V, 2-го типа.

Отметим, что матрицы жесткости и векторы узловых сил двухсеточного конечного элемента Va ( V ) определяются в локальной декартовой системе координат Оху2 (рис. 3). Процедуры построения криволинейных двухсеточных конечных элементов 1-го и 2-го порядков двух типов, имеющих форму, как двухсеточный конечный элемент Va (V,), аналогичны выше описанным процедурам.

3. Результаты расчетов. Рассмотрим в глобальной декартовой системе координат О2 х0 у0 20 расчет композитной прямоугольной цилиндрической панели V0, угол раствора которой а0 = п /2 ; О2х0у0 — плоскость симметрии.

где м°, V0, — перемещения узла і двухсеточного

конечного элемента V,; і = 1,...,32.

С помощью матрицы вращений [Т ], размерности 96 х 96 , между векторами ёа, ё0 установим связь ёа = [Т ]ё1. Матрица [Т ] в силу (8), (19) имеет вид

[Ts ] =

[Els ] [E0 ] [E2;

[E0 ] [Es ] [E0:

-[E2 ] [E0 ] [Els;

(20)

где [E0] ,..., [ ] — матрицы размерности 32x32 ;

[E0 ] — нулевая и [Es ] — единичная матрицы; [E; ] = cos Д [Es ], [E2 ] = sin fis [Es ].

Используя равенство ёа = [Г ]ё0, построим соотношения [3, с. 237]

[Ks0] = [Ts]T[Ks][Ts], Fs0 = [Ts]TFs

(2i)

Рис. 5. Панель У0

При у0 = 0 панель закреплена. Базовая дискретная модель панели состоит из однородных изотропных конечного элемента Уе (рис. 1), которая порождает криволинейную сетку г,},к размерами 82x82x19 (рис. 5). Для конечного элемента Уе имеем: ае = а0 / 81 , к‘х = аеЩ , кеу = Ь /81 = 0,5 ,

к = к/18 = 0,278, где Ь = 40,5 , к = 5, Ь — длина;

к — толщина панели. Внешний и внутренний радиусы панели равны 29 и 24. Двухсеточная дискретная модель панели состоит из двухсеточного конечного элемента У, 3-го порядка 2-го типа размерами

9к‘ х9ке х9к (рис. 3), к = 9к‘ , ка = 9ке .

х у г х г г ? у у

Двухсеточный конечный элемент У, армирован волокнами (с криволинейным сечением со сторонами ке , кге ), направленными по оси Оу , сечения волокон закрашены (рис. 4). Расстояние между волокнами по дуговой координате равно 2к , по толщине двух-

где [К,0] — матрица жесткости; Б,0 — вектор узловых сил двухсеточного конечного элемента У,, которые определяются в глобальной декартовой системе координат О2 х0 у0 г0 панели У0.

Для двухсеточного конечного элемента У, находим угол в, и по формуле (20) вычисляем матрицу [Т ]; , = 1,...,N ; N — общее число двухсеточного конечного элемента У,. Используя (21), определяем матрицы жесткости [К,0] и векторы узловых сил Б,0 двухсеточного конечного элемента У,, , = 1,...,N, с помощью которых строим систему уравнений метода конечных элементов для дискретной модели панели У0 .

Анализ результатов показывает, что максимальные перемещения двухсеточной шк = 3,694 и базовой

= 4,002 моделей отличаются на 7,7%. Максимальные эквивалентные напряжения двухсеточной ак = 0,424 и базовой а0 = 0,431 моделей отличаются на 1,6 % . Число неизвестных базовой модели равно 378594, двухсеточной — 4968. Ширина ленты системы уравнений метода конечных элементов базовой модели равна 9474, двухсеточной — 1206. Двухсеточная модель занимает в 600 раз меньше памяти ЭВМ, чем базовая. Время реализации метода конечных элементов для двухсеточной модели в 27 раз меньше, чем для базовой модели.

Расчеты показывают, что для однородной изотропной панели У0 (рис. 5) погрешность для макси-

мальных перемещений (по оси О2 г0) равна 3,0%, для максимальных эквивалентных напряжений —

11,6%, т. е. двухсеточный конечный элемент У, мож-

но использовать для анализа деформирования однородных изотропных цилиндрических оболочек и панелей.

Библиографический список

1. Матвеев А.Д. Некоторые подходы проектирования упругих многосеточных конечных элементов // Деп. в ВИНИТИ № 2990-В00 / Институт вычислительного моделирования СО РАН. — Красноярск, 2000.

2. Матвеев А.Д. Многосеточное моделирование композитов нерегулярной структуры с малым коэффициентом наполнения // ПМТФ. — № 3. — 2004.

3. Зенкевич О. Метод конечных элементов в технике. — М., 1975.

4. Норри Д., де Фриз Ж. Введение в метод конечных элементов. — М., 1981.

5. Самуль В.И. Основы теории упругости и пластичности. — М., 1982.