CC BY
35
____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА
Том 156, кн. 4 Физико-математические науки
2014
УДК 539.3
МНОГОСЕТОЧНЫЕ КРИВОЛИНЕЙНЫЕ ЭЛЕМЕНТЫ В ТРЕХМЕРНОМ АНАЛИЗЕ ЦИЛИНДРИЧЕСКИХ КОМПОЗИТНЫХ ПАНЕЛЕЙ И ОБОЛОЧЕК С ПОЛОСТЯМИ И ОТВЕРСТИЯМИ
А.Д. Матвеев, А.Н. Гришанов
Аннотация
Представлены процедуры построения криволинейных двухсеточных конечных элементов (ДвКЭ) и сложных многосеточных конечных элементов (МнКЭ) для расчета трехмерных упругих статически нагруженных композитных цилиндрических панелей, оболочек с полостями и отверстиями сложной формы. Базовые дискретные модели трехмерных композитных панелей, оболочек, учитывающие их неоднородную (микронеоднородную) структуру и сложную форму, имеют очень высокую размерность. ДвКЭ и сложные МнКЭ проектируются на основе базовых дискретных моделей панелей и оболочек с применением интерполяционных полиномов и уравнений трехмерной задачи теории упругости, записанных в локальных декартовых системах координат данных элементов. Предлагаемые элементы описывают трехмерное напряженное состояние в композитных панелях и оболочках, учитывают их структуры и порождают дискретные модели малой размерности. Расчеты панелей волокнистой структуры с полостями показывают, что максимальные перемещения и эквивалентные напряжения базовых и двухсеточных (многосеточных) дискретных моделей панелей отличаются на 1^12 %. Реализация метода конечных элементов для двух- и многосеточных дискретных моделей панелей требует в 103 + 104 раз меньше объема памяти ЭВМ и в 102 + 103 раз меньше временных затрат, чем для базовых.
Ключевые слова: композиты, упругость, цилиндрические оболочки и панели, сложные многосеточные и двухсеточные криволинейные конечные элементы.
Введение
Как известно, общий недостаток теорий деформирования упругих композитных цилиндрических панелей и оболочек заключается в том, что в их основе лежат гипотезы (определяющие законы перемещений и напряжений), которые порождают неустранимую погрешность в приближенных решениях. В частности, существующие теории не учитывают сложный характер закреплений, например частичное закрепление по толщине толстых панелей и оболочек, не всегда достаточно точно описывают деформирование композитных панелей и оболочек с полостями и отверстиями сложной формы при действии локальных нагружений. С помощью этих теорий затруднительно исследовать напряженное состояние трехмерных композитных панелей и оболочек, которые подкреплены ребрами жесткости, имеют переменную толщину, нерегулярную структуру, малый коэффициент наполнения.
В настоящей работе представлены процедуры построения в локальных декартовых системах координат криволинейных ДвКЭ и сложных МнКЭ, которые используются для расчета по методу конечных элементов (МКЭ) линейно упругих трехмерных композитных цилиндрических панелей и оболочек с полостями и отверстиями сложной формы при статическом нагружении. Сложные МнКЭ и ДвКЭ формы прямоугольного параллелепипеда, применяемые для анализа упругих тел
47
48
А.Д. МАТВЕЕВ, А.Н. ГРИШАНОВ
неоднородной структуры, построены в работах [1—3]. Процедуры построения трехмерных криволинейных композитных ДвКЭ в локальных декартовых системах координат с применением известных интерполяционных полиномов 1-го, 2-го и
3-го порядков изложены в [4, 5]. В отличие от ДвКЭ, рассмотренных в работах [4, 5], здесь рассмотрены процедуры построения криволинейных композитных ДвКЭ, имеющих полости и отверстия сложной формы. Полости и отверстия применяются в панелях и оболочках, например, для уменьшения веса, размещения электронного оборудования, различного рода механизмов и т. д. Как известно, базовые дискретные модели трехмерных композитных панелей и оболочек, которые учитывают их сложную форму, неоднородные и микронеоднородные структуры, имеют очень высокую размерность. Предлагаемые криволинейные элементы проектируются на основе базовых (мелких) разбиений на конечные элементы (КЭ) панелей и оболочек. Для построения композитного криволинейного ДвКЭ, который имеет полости (отверстия), применяем две вложенные криволинейные сетки: мелкую и крупную. Мелкая сетка порождена базовым разбиением ДвКЭ, которое состоит из криволинейных однородных КЭ 1-го порядка, учитывает его неоднородную (микронеоднородную) структуру и сложную форму полостей и отверстий. На мелкой сетке определяем крупную сетку; узлы крупной сетки ДвКЭ расположены вне областей полостей и отверстий. Представлены две процедуры построения трехмерных криволинейных композитных ДвКЭ.
Согласно первой процедуре на базовом разбиении (то есть на мелкой сетке) строим функционал Wa полной потенциальной энергии ДвКЭ Va в матричной форме [6, 7]. С помощью функций перемещений ДвКЭ Va, построенных по МКЭ на крупной сетке, узловые перемещения мелкой сетки выражаем через узловые перемещения крупной. Минимизируя функционал Wa по узловым перемещениям крупной сетки, получаем формулы для вычисления матрицы жесткости и вектора узловых сил ДвКЭ Va .
Суть второй процедуры заключается в следующем. На базовом разбиении ДвКЭ Vb с помощью метода конденсации [6] строим криволинейный суперэлемент. Вершины суперэлемента совпадают с узлами крупной сетки ДвКЭ Vb. Функционал Wb полной потенциальной энергии суперэлемента представляем в матричной форме. C помощью функций перемещений ДвКЭ Vb, построенных по МКЭ на крупной сетке, в функционале Wb узловые перемещения суперэлемента выражаем через узловые перемещения крупной сетки. Из минимизации функционала Wb по узловым перемещениям крупной сетки получаем формулы для вычисления матрицы жесткости и вектора узловых сил ДвКЭ Vb.
Предложена процедура построения криволинейных сложных МнКЭ, которые проектируются с применением композитных криволинейных ДвКЭ. Кратко изложена процедура построения однородных трехмерных криволинейных односеточных КЭ 1-го, 2-го и 3-го порядков в локальных декартовых системах координат. При построении ДвКЭ и сложных МнКЭ используем криволинейные однородные КЭ 1-го порядка. На рис. 1 представлен односеточный однородный криволинейный КЭ Ve 1-го порядка с характерными размерами hex х hzy х hez, где az - угол раствора КЭ Ve, O\x\y\z\ - локальная декартова система координат, z\O\y\ -плоскость симметрии, cd - ось цилиндрической панели, оболочки, hez - толщина, hy - длина КЭ Vz; Rf, Rf - радиусы нижней и верхней поверхностей КЭ Vz, узлы отмечены точками (8 узлов). Прямоугольники размерами hZ, х hzy есть боковые грани, криволинейные прямоугольники - торцевые грани КЭ Vz . Форма КЭ Vz - прямая призма высотой hzy. Отметим, что для учета неоднородной (микронеоднородной) структуры, сложных форм полостей и отверстий необходимо использовать мелкие разбиения панелей и оболочек. Так как при мелком разбиении угол
МНОГОСЕТОЧНЫЕ КРИВОЛИНЕЙНЫЕ ЭЛЕМЕНТЫ.. .
49
а
7,
б
Рис. 1. КЭ Ve 1-го порядка
Рис. 2. КЭ Ve 2-го и 3-го порядков
раствора ае мал (рис. 1), то форма КЭ Ve мало отличается от формы прямоугольного параллелепипеда.
В связи с этим при построении по МКЭ функций перемещений для однородных криволинейных КЭ Ve 1-го, 2-го и 3-го порядков используем соответственно известные интерполяционные полиномы 1-го, 2-го и 3-го порядков (которые используют при построении функций перемещений для КЭ 1-го, 2-го и 3-го порядков формы прямоугольного параллелепипеда [6, 7]) и уравнения трехмерной задачи теории упругости, записанные в локальных декартовых системах координат Oyx\y\z\ данных КЭ (рис. 1, 2). Таким образом, в КЭ Ve реализуется трехмерное напряженное деформированное состояние. Поскольку при построении основных соотношений для криволинейных ДвКЭ и сложных МнКЭ используем КЭ Ve 1-го порядка (рис. 1), то ДвКЭ и сложные МнКЭ также описывают трехмерное напряженное состояние в композитных панелях и оболочках. Матрицы жесткости и векторы узловых сил криволинейных ДвКЭ и сложных МнКЭ определяем в локальных декартовых системах координат, а системы уравнений МКЭ для дискретных моделей оболочек и панелей - в глобальных декартовых системах координат. Связь между локальными и глобальными системами координат осуществляем с помощью матриц вращений [6], которые определяем только для узловых перемещений ДвКЭ и сложных МнКЭ.
1. Односеточные однородные криволинейные конечные элементы
Кратко рассмотрим процедуру построения в локальной декартовой системе координат O\X\y\Z\ матрицы жесткости и вектора узловых сил для односеточного однородного криволинейного КЭ Ve 1-го порядка, рис. 1 (который используется при построении мелких разбиений панелей, оболочек). При аппроксимации функций перемещений ue, ve, we КЭ Ve применяем полиномы 1-го порядка (записанные
50
А.Д. МАТВЕЕВ, А.Н. ГРИШАНОВ
в декартовой системе координат O\x\y\z\) вида
ue, ve, we = ai + a2'x\ + азу1 + a4Zi + a^xiyi + agzixi + a7Ziyi + asx\y\z\, (1)
где ai = const, i = 1,..., 8.
Используя (1), с помощью МКЭ [6, 7] строим функции перемещений ue, ve, we для КЭ Ve, которые запишем в виде
^2 Niui , ve = ^2 NiVi , We = ^2 NiWl ,
i=1
i=1
i=1
где ul, v1 *, wl * - перемещения i-го узла; Ni - функция формы i-го узла КЭ Ve, i =1,..., 8.
Вектор-функцию перемещений Ue = {ue, ve, we}T КЭ Ve представим в матричной форме
Ue = [Ne ]Sl (2)
где [Ne] - матрица функций формы КЭ Ve, S1 - вектор узловых перемещений КЭ Ve .
Вектор Si, отвечающий декартовой системе координат Oixiyizi, представим в форме
(3)
Sl = {u1,...,u8,vl,...,v8,wl,...,ws}T,
где индекс T означает транспонирование.
Для полной потенциальной энергии КЭ Ve составим выражение
We = ~ eTVe dV - UT Fe dV - UTqe dS,
(4)
Ve
Ve
u
e
где ee {exi ,eyi, ezi T^xiyi T^xizi ,1yi zi } , Ve {vxi, Vyi, Vzi ,Txiyi ,TxiZi ,TyiZi }
вектор-функции деформаций и напряжений (отвечающие трехмерной задаче теории упругости [7]); Fe, qe - векторы объемных и поверхностных сил КЭ Ve; Ve, Se - область и поверхность КЭ Ve соответственно.
Соотношения Коши и закон Гука для КЭ Ve имеют вид [7]
ee = [Be] Si, Ve = [De][Be] S\, (5)
где [Be], [De] - матрицы деформаций и модулей упругости КЭ Ve.
Подставляя (2), (5) в функционал (4) и минимизируя его, получаем
[Kl ] = J [Be]T [De][Be] dV, Pe = J [Ne]T Fe dv + J [Ne]T qe dS, (6)
Ve Ve Se
где [Kl] - матрица жесткости и P l - вектор узловых сил КЭ Ve.
Процедуры построения однородных односеточных криволинейных КЭ Ve 2-го и 3-го порядков (рис. 2), имеющих такие же геометрические формы, как и КЭ Ve 1-го порядка (рис. 1), аналогичны вышеописанной.
2. Композитные криволинейные ДвКЭ с полостями и отверстиями
Рассмотрим две процедуры построения в локальных декартовых системах ко-
ординат основных соотношений (матриц жесткости и узловых усилий) для трех-
мерных криволинейных композитных ДвКЭ с полостями и отверстиями сложной
формы.
МНОГОСЕТОЧНЫЕ КРИВОЛИНЕЙНЫЕ ЭЛЕМЕНТЫ.. .
51
Рис. 3. ДвКЭ vn , vn
Рис. 4. Сечение ДвКЭ V£ , V^ (суперэлемента Ga)
2.1. Первая процедура построения ДвКЭ. Основные положения первой процедуры (не теряя общности суждений) рассмотрим на примере построения ДвКЭ УП, имеющего полость сложной формы. На рис. 3 показан ДвКЭ УП в локальной декартовой системе координат Oxyz c характерными размерами ЬП х bn х ЬП, где zOy - плоскость симметрии, ЬП - толщина, by - длина ДвКЭ УП, n = !,..., No, No - общее число ДвКЭ УП в дискретной модели панели (оболочки).
Считаем, что ДвКЭ Vya армирован непрерывными волокнами, направленными вдоль оси Oy. Между компонентами неоднородной структуры ДвКЭ У^ связи идеальны. Компоненты есть изотропные однородные тела. Функции перемещений, напряжений и деформаций компонентов ДвКЭ УП удовлетворяют закону Гука и соотношениям Коши, которые отвечают трехмерной задаче теории упругости. Область ДвКЭ У£ представляем базовым разбиением Ra, которое состоит из однородных КЭ Уе 1-го порядка (рис. 1) с характерными размерами hex х hzy х hez, пусть hy, hez, ae = const, e = 1,...,M, M - общее число КЭ Ve базового разбиения Ra. Имеем by = 9hzx, by = 9hy, ЬП = 9hzz. В центре области ДвКЭ У^
расположена полость сложной формы с характерными размерами 3h| х 3hy х 3hZ, то есть полость имеет форму криволинейного параллелепипеда.
Базовое разбиение Ra учитывает неоднородную (микронеоднородную) структуру ДвКЭ УПа, сложную форму полости и порождает мелкую криволинейную трехмерную сетку ha. Поскольку в КЭ Vz (рис. 1) реализуется трехмерное напряженно-деформированное состояние (см. п. 1), то в области ДвКЭ Уу также реализуется трехмерное напряженное состояние. Отметим, что базовые разбиения Ra двухсеточных КЭ УПа (которые представляют область панели, оболочки) образуют базовую дискретную модель панели, оболочки. На рис. 4 показано поперечное сечение (при y = 4hy) ДвКЭ Vya , которое представлено мелкой сеткой базового разбиения, сечения волокон закрашены. На мелкой сетке ha определяем крупную криволинейную трехмерную узловую сетку Ha . Считаем, что угол раствора ДвКЭ Уу мал (что характерно при использовании мелких разбиений Ra или когда ДвКЭ Vya представляют пологую панель, оболочку). В этом случае форма ДвКЭ Vya (рис. 3) мало отличается от формы прямоугольного параллелепипеда. В связи с этим при построении по МКЭ на крупной сетке Ha функций перемещений для ДвКЭ Vya будем использовать известный интерполяционный полином 3-го порядка (который применяют при построении функций перемещений для КЭ 3-го порядка формы прямоугольного параллелепипеда [6, 7]), записанный в декартовой системе координат Oxyz, рис. 3. В этом случае ДвКЭ Vya будем называть ДвКЭ 3-го порядка, то есть порядок ДвКЭ равен порядку интерполяционного полинома, построенного
52
А.Д. МАТВЕЕВ, А.Н. ГРИШАНОВ
на крупной сетке. Крупная сетка Ha по характеру расположения и количеству узлов отвечает сетке КЭ 3-го порядка в виде прямоугольного параллелепипеда [6, 7]. Узлы сетки Ha на рис. 3 отмечены точками (32 узла). На крупной сетке Ha по МКЭ строим функции перемещений ua, va, wa для ДвКЭ , используя интерполяционные полиномы 3-го порядка (записанные в декартовой системе координат Oxyz) вида
Ua, Va, Wa = ai + a2x + a3y + aAz + a5xy + a6zx + arzy + asxyz+
i 2, 2, 2, 2, 2 i 2 i 2 , 2,
+ agx + aioy + anz + ai2x y + ai3xz + auy x + a15y z + a16z x+
2 3333 3 3 3 3
+ ayjz y + aigx + aig y + a2oz + a2ix y + a22 xz + a23y x + a24y z + a25z x+
3 2 2 2 3 3 3
+ a26z y + a2rx yz + a28y xz + a2gz xy + a3ox yz + a3iy xz + a32z xy.
Функции перемещений ua , va , wa представляем в форме
32 32 32
ua = Ni ui , va = Y.nm, wa =Y1 (7)
i=i i=1 i=1
где ua, va , wa , Nij1 - перемещения и функция формы i-го узла крупной сетки Ha, i = 1,..., 32.
Пусть ось Oiyi локальной декартовой системы координат Oixiyizi КЭ Ve, см. рис. 1 (расположенного в базовом разбиении Ra ДвКЭ VП), параллельна оси Oy локальной декартовой системы координат Oxyz ДвКЭ Vna (рис. 3) и пусть между осями Oixi, Ox угол равен /Зе. Обозначим через Se вектор узловых перемещений КЭ Ve, e = 1,... ,M, отвечающий декартовой системе координат Oxyz. Вектор Se представим в виде
S
e
{uei,...,u8nve,...,v8,we1,...,w8}T
(8)
где u, vf, wf - перемещения i-го узла КЭ Ve, i = 1,..., 8, e = 1,... ,M. Векторы Si (см. представление (3)) и Se связаны соотношением [6]
Si = [Te]Se
Здесь [Te] - матрица вращений [6] размерности 24 х 24, которая в силу имеет структуру
[Te]
[Ei] [E0] [E2]
[Eo] [E] [Eo] ,
-[E2] [Eo] [Ei]
(9) (3), (8)
где подматрицы имеют размерность 8 х 8; [Eo] - нулевая и [E] - единичная матрицы, [Ei] = cos pe[E], [E2] = sin pe[E], e = 1,...,M.
Используя связь (9) между векторами S;) и Se, получим соотношения [7]
[Ke] = [Te]T [Ki][Te], Pe = [TPf Pi,
где [Ki], [Ke] - матрицы жесткости и P;), Pe - векторы узловых сил КЭ Ve, отвечающие соответственно декартовым системам координат O1x1y1z1 и Oxyz , e = 1,...,M.
Полную потенциальную энергию Wa для базового разбиения Ra запишем в форме
м
Wa =
i=i
где M - общее число КЭ Ve базового разбиения Ra .
I ST [Ke]Se - STP(
(10)
МНОГОСЕТОЧНЫЕ КРИВОЛИНЕЙНЫЕ ЭЛЕМЕНТЫ.. .
53
Обозначим через Sa вектор узловых перемещений крупной сетки Ha (отвечающий декартовой системе координат Oxyz, рис. 3). С помощью функции перемещений (7) узловые перемещения (вектора Se) КЭ Ve выражаем через узловые перемещения крупной сетки Ha, то есть вектор Se представляем через вектор Sa, в результате получим равенство
Se =[A<a\ Sa,
где [Aa] - прямоугольная матрица размерности 24 х 96, e = 1,...,M. Подставляем (11) в функционал (10) и, минимизируя его, получаем
M
M
[Ka] [Ke][Aae], Fa = P
e=1
e=1
e
(и)
где [Ka] - матрица жесткости, Fa - вектор узловых сил криволинейного ДвКЭ
V a
V n '
2.2. Вторая процедура построения ДвКЭ. Вторую процедуру рассмотрим на примере построения ДвКЭ Vnb 3-го порядка. Пусть ДвКЭ Vnb имеет размеры, форму, неоднородную (микронеоднородную) структуру и расположен в локальной декартовой системе координат Oxyz как и ДвКЭ Vna (рис. 3). В этой процедуре используем функции перемещений ua, va, wa вида (7) (см. п. 2.1), базовое разбиение Ra, мелкую ha и крупную Ha сетки ДвКЭ V£, другими словами ua , va, wa, Ha, Ra, Sa есть функции перемещений, крупная сетка, базовое разбиение и вектор узловых перемещений ДвКЭ V^, которые для удобства изложения обозначим соответственно через ub, vb, wb, Hb, Rb, Sb. На базовом разбиении Rb с помощью метода конденсации [6] строим криволинейный суперэлемент Gs, полную потенциальную энергию Ws которого запишем в виде
Ws = 2 Sf [Ks] Ss - Sf Ps, (12)
где [Ks], Ps, Ss - матрица жесткости, векторы узловых сил и перемещений суперэлемента Gs , определяемые в локальной декартовой системе координат Oxyz ДвКЭ vn.
Используя функции перемещений (7), узловые перемещения суперэлемента Gs выражаем через узловые перемещения крупной сетки Hb , то есть вектор Ss представляем через вектор Sb , в результате получим соотношение
Ss = [Abs]Sb, (13)
где [Abs ] - прямоугольная матрица.
Подставляем (13) в функционал (12) и, минимизируя его, получаем
[Kb] = [Abs]T [Ks][Abs], Fb =[Abs]T Ps,
где [Kb] - матрица жесткости, Fb - вектор узловых сил криволинейного ДвКЭ V^.
Замечание 1. Базовые разбиения (мелкие сетки) ДвКЭ V£, Vb могут иметь сколь угодно большую размерность. Это связано с тем, что для более точного учета неоднородной (микронеоднородной) структуры ДвКЭ Vn, Vnb, сложных форм полостей и отверстий необходимо использовать мелкие разбиения. Однако если размерность крупной сетки задана, то увеличение размерности мелкой сетки не приводит к увеличению числа узловых неизвестных ДвКЭ Vna , Vnb , то есть в этом случае размерности двухсеточных дискретных моделей не связаны со сложностью форм неоднородных структур панелей и оболочек.
54
А.Д. МАТВЕЕВ, А.Н. ГРИШАНОВ
Замечание 2. Если базовое разбиение ДвКЭ имеет высокую размерность, то построение суперэлемента Gs затруднительно, поскольку эта процедура связана с обращением матрицы высокого порядка. В этом случае сначала на базовом разбиении строим систему суперэлементов Gsa (а = 1,..., Ms, где Ms - общее число суперэлементов), которые покрывают область ДвКЭ VП. Затем, используя суперэлементы Gsa, строим суперэлемент Gs.
Замечание 3. Расчеты показывают, что процедура, приведенная в п. 2.2, позволяет построить более точное решение, чем решение, полученное с помощью процедуры из п. 2.1. Однако временные затраты на реализацию процедуры п. 2.2 больше, чем процедуры п. 2.1, что связано с необходимостью обращения матриц высокого порядка (в процедуре п. 2.2).
Процедуры построения криволинейных композитных ДвКЭ 1-го и 2-го порядка, имеющих такие же геометрические формы, как ДвКЭ V£ (ДвКЭ Vb), аналогичны процедурам п. 2.1 и п. 2.2.
2.3. Определение напряжений в ДвКЭ.
2.3.1. Рассмотрим процедуру определения напряжений в ДвКЭ VП (см. п. 2.1). Пусть найден вектор Sa узловых неизвестных ДвКЭ в локальной декартовой системе координат Oxyz (рис. 3). По формуле (11) определяем (в декартовой системе координат Oxyz) вектор Se узловых перемещений КЭ Ve, то есть Se = [^a]Sa, e = 1,...,M. Вектор узловых перемещений S1 КЭ Ve, отвечающий локальной декартовой системе координат Oixiyizi (рис. 1), находим по формуле (9), а именно S1 = [Te]Se, e = 1,..., M. С помощью функций перемещений (2) и соотношений (5) вектор ае функций напряжений КЭ Ve определяем по формуле ae = [^^DB^S1, e = 1,... ,M. Эквивалентные напряжения для КЭ Ve базового разбиения ДвКЭ Vna вычисляем по 4-й теории прочности в их центрах тяжести, e = 1,... ,M.
2.3.2. Рассмотрим процедуру определения напряжений в ДвКЭ Vb (см. п. 2.2). Пусть найден вектор Sb узловых неизвестных ДвКЭ Vb, отвечающий декартовой системе координат Oxyz (рис. 3). По формуле (13) определяем (в декартовой системе координат Oxyz ) вектор Ss узловых перемещений суперэлемента Gs , то есть Ss = [^s]Sb. Используя вектор Ss и соответствующие формулы метода конденсации [6], находим перемещения в узлах мелкой сетки, которые лежат внутри области суперэлемента Gs. Затем по алгоритмам, аналогичным п. 2.3.1, находим векторы функций напряжений и эквивалентные напряжения для КЭ Ve базового разбиения ДвКЭ Vb, e = 1,..., M. Так как компоненты композитных панелей и оболочек представляются КЭ Ve базовых разбиений ДвКЭ, то напряжения можно определить в любом компоненте неоднородных структур панелей и оболочек.
3. Криволинейные сложные МнКЭ с полостями и отверстиями
Основные положения процедуры построения композитного криволинейного сложного многосеточного конечного элемента (не теряя общности суждений) рассмотрим на примере построения сложного МнКЭ Vm с характерными размерами hm х h™ х hm, который имеет регулярную систему полостей сложной формы и расположен в локальной декартовой системе координат Oxoyozo (рис. 5), где hm -толщина, hm - длина МнКЭ Vm. Форма МнКЭ Vm есть прямая призма высотой hm. Область сложного МнКЭ Vm представлена криволинейными композитными ДвКЭ Vjb 3-го порядка с характерными размерами ЪП х х ЪП (рис. 3), n = = 1,..., N, N - общее число ДвКЭ V b, для рис. 5 имеем N = 27. Пусть ДвКЭ VП
МНОГОСЕТОЧНЫЕ КРИВОЛИНЕЙНЫЕ ЭЛЕМЕНТЫ.. .
55
Рис. 5. Сложный МнКЭ Vm
Рис. 6. Панель Vo
построены по процедуре, изложенной в п. 2.2. Базовое разбиение ДвКЭ Vjb состоит из криволинейных КЭ Ve 1-го порядка с характерными размерами hex х hey х hez (рис. 1). Это разбиение учитывает неоднородную (микронеоднородную) структуру и форму ДвКЭ Vn (то есть учитывает структуру и сложную форму МнКЭ Vm).
Пусть hy, hy, ae = const, e = 1,...,No; No - общее число КЭ Ve базового разбиения ДвКЭ Vb. ДвКЭ V^b расположен в локальной декартовой системе координат Oxyz (рис. 3), где ЪП = 9hex, Ъгу = 9hey, ЪП = 9hy . Пусть Ъгу, ЪП, jn = const, In - угол раствора ДвКЭ , n = 1, . . . , N. Имеем hT = 3ЪП, hmi = 3Ъгу, hT = 3ЪП . ДвКЭ Vnb армирован непрерывными волокнами с характерным сечением hex х hez, направленными вдоль оси Oy. В центре области ДвКЭ Vnb расположена полость сложной формы с характерными размерами 3hex х 3hy х 3hy, то есть полость имеет форму криволинейного параллелепипеда. На рис. 4 показано поперечное сечение ДвКЭ Vnb, сечения волокон закрашены. Отметим, что МнКЭ Vm включает множества криволинейных мелких hby и крупных Hn вложенных сеток ДвКЭ , n = 1,... ,N. Мелкие сетки hby образуют мелкую сетку hm сложного МнКЭ Vm. На мелкой сетке hm определяем крупную сетку Hm сложного МнКЭ Vm . Считаем, что угол раствора МнКЭ Vm мал.
При построении по МКЭ функций перемещений для МнКЭ Vm на крупной сетке Hm используем интерполяционный полином 3-го порядка (см. п. 2.1), записанный в локальной декартовой системе координат Oxoyozo (рис. 5). В этом случае элемент Vm будем называть сложным МнКЭ 3-го порядка. Узлы крупной сетки Hm МнКЭ Vm на рис. 5 отмечены точками (32 узла). Функции перемещений um, vm, wm МнКЭ Vm, построенные на крупной сетке Hm, запишем в виде
um
32
J^NTum,
i=1
vm
32
YjNTvT,
i=1
wm
32
^NTwy
^-•i ^T,
(14)
где NT - функция формы i -го узла сетки Hm; uT, vT, wT - значения перемещений в i-м узле (определяемые в декартовой системе координат Oxoyozo), i = 1,..., 32.
Обозначим через ST вектор узловых неизвестных (размерности 96) крупной сетки HT , то есть ST = {и"1,,.. ., uT, vT, .. ., v^^T,. .., w3[2}T . Пусть ось Oy локальной декартовой системы координат Oxyz ДвКЭ (рис. 3) параллельна оси Oyo локальной декартовой системы координат Oxoyozo МнКЭ VT (рис. 5) и пусть между осями Ox и Oxo угол равен уn. Векторы Sn, Sn узловых перемещений ДвКЭ Vjb (рис. 3), отвечающие соответственно декартовым системам координат
56
А.Д. МАТВЕЕВ, А.Н. ГРИШАНОВ
Oxyz и Oxoyozo, имеют структуры
Sn {u1, ■ ■ ■ , u32,v1, ■ ■ ■ , v32, W1,
So — UP u0 vo v0 w0
°n — lUli ■ ■ ■ , u32, v11 ■ ■ ■ ; v32,w1
iT
,w32} , ,W°°2}T ■
Между векторами 5n, бП имеем связь
Sn — К] бП,
(15)
(16)
где [Tn] - матрица вращений [6] размерности 96 х 96, которая в силу (15) имеет стРУктУРУ
' [М1] [Mo] [M2]
[T0] — [Mo] [M ] [Mo]
-[M2 ] [Mo] [M1I
Здесь подматрицы имеют размерность 32 х 32; [Mo] - нулевая и [M] - единичная матрицы, [M1] — cospn[M], [M2] — sinpn[M].
Используя связь (16) между векторами 5n и бП, получим соотношения [7]
К] — [тП]т Km, рП — [тП]т Pn,
где [Kn], [КП] - матрицы жесткости, Pn, РП - векторы узловых сил ДвКЭ Vjb, отвечающие соответственно декартовым системам координат Oxyz и Ox0y0z0 . Полную потенциальную энергию Wm МнКЭ Vm определим выражением
Wm —
Е (\(s0n)TКб - (бП)тPn) ■
(17)
Используя функции перемещений (14), выражаем узловые перемещения ДвКЭ VП через узловые перемещения крупной сетки Hm, в результате получим соотношение
бП — [Am]Sm,
где [Am] - квадратная матрица размерности 96 х 96.
Подставляя (18) в функционал (17) и минимизируя его, получим
N
N
[Km] — J2[Amf[кП][Ат], рт — Y^AmVр
0,
n
n=1
n=1
(18)
где [Km], Fm - матрица жесткости и вектор узловых сил сложного МнКЭ Vm.
4. Результаты расчетов
В качестве модельной задачи рассмотрим расчет композитной консольной прямоугольной в плане панели Vo регулярной волокнистой структуры c малым коэффициентом наполнения (рис. 6). Панель Vo состоит из пластичных материалов, имеет статическое нагружение, регулярную систему полостей сложной формы и расположена в декартовой системе координат Oxyz. Левый торец панели жестко закреплен, то есть при у — 0 имеем u,v,w — 0. Непрерывные волокна параллельны оси Оу. Базовое разбиение панели Vo, состоящее из однородных криволинейных КЭ Ve 1-го порядка с характерными размерами hex х hey х hez (рис. 1), учитывает ее неоднородную структуру, форму и порождает криволинейную мелкую сетку ha размерности 73 х 145 х 19, где e — 1, ■ ■ ■ ,M, M - общее число КЭ Ve базового разбиения панели Vo, M — 186624. Пусть hey, hez, ae — const, e — 1, ■ ■ ■ ,M. Двухсеточная дискретная модель панели Vo состоит из криволинейных ДвКЭ VП 3-го
МНОГОСЕТОЧНЫЕ КРИВОЛИНЕЙНЫЕ ЭЛЕМЕНТЫ.. .
57
порядка с характерными размерами 18hX х 18hy х 18hez, которые построены по процедуре п. 2.2 и мелкие сетки которых имеют размерность 19 х 19 х 19, n = = 1,..., No, No - общее число ДвКЭ УП, для панели Vo имеем No = 32. ДвКЭ УП расположен в локальной декартовой системе координат O\x\y\z\, причем ось Oiyi параллельна оси Oy декартовой системы координат Oxyz панели.
Базовые разбиения ДвКЭ состоят из КЭ Ve (рис. 1). Область ДвКЭ У^ состоит из 8 суперэлементов Gsa с характерными размерами 9hex х 9hy х 9hez, а = = 1,..., 8. Сначала строим суперэлементы Gsa. Затем, используя суперэлементы Gsa, строим суперэлемент Gs для всей области ДвКЭ VП (см. замечание 2). В центре каждого суперэлемента Gsa расположена полость формы криволинейного параллелепипеда с характерными размерами 3hzx х 3hy х 3hzz . На рис. 4 представлено сечение суперэлемента Gsa в плоскости, перпендикулярной оси Oiyi, при yi = 4hy . Сечение представлено сеткой базового разбиения, сечения волокон с характерными размерами hex х hez заштрихованы. Для узлов мелкой сетки введена целочисленная система координат ijk (рис. 6). В узлах p(i,j,k) мелкой сетки панели, где i, j, k - целочисленные координаты, k = 19, i = 1 + 6(а — 1), j = 73 + 6(в — 1), а = 1,..., 7, в = 1,..., 13, действуют вертикальные силы q. На рис. 6 поверхность панели, на которой задано нагружение, заштрихована (1/4 часть верхней поверхности панели).
Модуль Юнга связующего материала равен 1, волокна - 10, коэффициент Пуассона для волокна и связующего материала равен 0.3; q = 0.1. Радиус внутренней поверхности панели равен 25, внешней - 30, толщина панели h = 5, длина панели L = 40 (рис. 6). Угол раствора панели равен п/4, коэффициент наполнения панели равен 0.31.
Результаты расчетов панели Vo показывают, что максимальное перемещение (эквивалентное напряжение) двухсеточной дискретной модели панели Vo, состоящей из ДвКЭ УП1, отличается от максимального перемещения (эквивалентного напряжения) базовой дискретной модели Ro на 2.46% (на 8.8%). Размерность базовой модели Ro равна 593040, ширина ленты системы уравнений (СУ) МКЭ равна 8448. Двухсеточная модель R^ имеет 1344 узловых неизвестных (то есть имеет в 440 раз меньше, чем неизвестных в базовой модели Ro), ширина ленты СУ МКЭ равна 408 (в 20 раз меньше ширины ленты СУ МКЭ модели Ro). Реализация МКЭ для двухсеточной модели R^ требует примерно в 9100 раз меньше объема памяти ЭВМ, чем для базовой модели Ro . Эквивалентные напряжения определяются по
4-й теории прочности.
Заключение
В работе для расчета трехмерных линейно упругих композитных цилиндрических панелей и оболочек с различными коэффициентами наполнения, имеющие статическое нагружение, полости и отверстия сложной формы, предложены криволинейные ДвКЭ и сложные МнКЭ, при построении которых используются известные интерполяционные полиномы 1-го, 2-го и 3-го порядков и уравнения трехмерной задачи теории упругости, записанные в локальных декартовых системах координат данных элементов. Достоинства предлагаемых элементов состоят в следующем. ДвКЭ и сложные МнКЭ в композитных панелях и оболочках с различными коэффициентами наполнения:
• описывают трехмерное напряженное состояние;
• учитывают неоднородную и микронеоднородную структуры;
• учитывают полости и отверстия сложной формы;
• учитывают сложный характер закреплений.
ДвКЭ и сложные МнКЭ порождают двух- и многосеточные дискретные модели
58
А.Д. МАТВЕЕВ, А.Н. ГРИШАНОВ
композитных панелей и оболочек, размерности которых в 102 т 104 раз меньше размерностей базовых моделей, Отметим, что напряжения можно определить в любом компоненте неоднородных структур панелей и оболочек. Реализация МКЭ для двух- и многосеточных дискретных моделей требует в 102 т 103 раз меньше временных затрат, чем для базовых.
Работа выполнена при финансовой поддержке РФФИ (проект № 14-01-00130).
Summary
A.D. Matveev, A.N. Grishanov. Multigrid Curvilinear Elements in the Three-Dimensional Analysis of Composite Cylindrical Panels and Shells with Cavities and Openings.
The paper presents the procedure of construction of curvilinear double-grid finite elements (DgFE) and complex multigrid finite elements (MgFE) for calculating three-dimensional elastic composite cylindrical panels and shells with cavities and openings. The standard discrete models of three-dimensional composite panels and shells, taking into account their inhomogeneous (micro-inhomogeneous) structure and complex form, have very high dimension. DgFE and complex MgFE are built based on the standard discrete models of panels and shells using the interpolating polynomials and equations of the three-dimensional problem of elasticity, written in the local Cartesian coordinate systems of these elements. DgFE and complex MgFE describe the three-dimensional stress state in the composite panels and shells taking into account their structure and thus generating the discrete models of small dimensions. Calculations of panels with fibrous structure demonstrate that the maximum displacement and equivalent stress of the basic and double-grid (multigrid) discrete panel models differ by 1 T 12 %. Implementation of the finite element method for two- and multigrid discrete models of panels requires 103 T 104 times less computer memory and 102 T 103 times less time than for the basic model.
Keywords: composites, elasticity, cylindrical shells and panels, complex multigrid and double-grid curvilinear finite elements.
Литература
1. Матвеев А.Д. Некоторые подходы проектирования упругих многосеточных конечных элементов. - Красноярск: Ин-т вычисл. моделирования СО РАН. - 2000. — 30 с. — Деп. в ВИНИТИ № 2990-В00.
2. Матвеев А.Д. Многосеточное моделирование композитов нерегулярной структуры с малым коэффициентом наполнения // Прикл. механика и техн. физика. — 2004. — № 3. — С. 161—171.
3. Матвеев А.Д. Построение сложных многосеточных элементов с неоднородной и микронеоднородной структурой // Изв. Алт. гос. ун-та. — 2014. — № 1/1. — С. 80—83.
4. Матвеев А.Д., Гришанов А.Н. Двухсеточное моделирование цилиндрических оболочек и панелей переменной толщины // Вестн. Краснояр. гос. аграр. ун-та. — 2014. — № 4. — С. 90—97.
5. Матвеев А.Д., Гришанов А.Н. Одно- и двухсеточные криволинейные элементы трехмерных цилиндрических панелей и оболочек // Изв. Алт. гос. ун-та. — 2014. — № 1/1. —
С. 84—89.
6. Норри Д, де Фриз Ж. Введение в метод конечных элементов. — М.: Мир, 1981. — 304 с.
7. Зенкевич О. Метод конечных элементов в технике. — М.: Мир, 1975. — 541 с.
Поступила в редакцию
29.10.14
МНОГОСЕТОЧНЫЕ КРИВОЛИНЕЙНЫЕ ЭЛЕМЕНТЫ.. .
59
Матвеев Александр Данилович - кандидат физико-математических наук, старший научный сотрудник отдела № 5, Институт вычислительного моделирования СО РАН, г. Красноярск, Россия.
E-mail: [email protected]
Гришанов Александр Николаевич - преподаватель кафедры «Прочность летательных аппаратов», Новосибирский государственный технический университет, г. Новосибирск, Россия.
E-mail: [email protected]
