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

УДК 517: 539.4

Трехмерные композитные многосеточные конечные элементы оболочечного типа

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

'Институт вычислительного моделирования СО РАН (Красноярск, Россия), 2Новосибирский государственный технический университет (Новосибирск, Россия)

Three-Dimensional Composite Multigrid Finite Shell-Type Elements

A.D. Matveev1, A.N. Grishanov2

''Institute of Computational Modeling of the Siberian Branch of the RAS (Krasnoyarsk, Russia)

2Novosibirsk State Technical University (Novosibirsk, Russia)

Предложены процедуры построения криволинейных трехмерных композитных многосеточных конечных элементов (МнКЭ) оболочечного типа для расчета напряженного состояния упругих цилиндрических оболочек, имеющих неоднородную (микронеоднородную) структуру и статическое нагружение. МнКЭ проектируются в локальных декартовых системах координат на основе мелких (базовых) разбиений (моделей) оболочек. При построении МнКЭ (без увеличения их размерности) можно использовать сколь угодно мелкие базовые разбиения оболочек, что позволяет в рамках микроподхода учитывать их неоднородную и микронеоднородную структуру, сложную форму, сложный характер нагруже-ний и закреплений. Напряженно-деформированное состояние в МнКЭ описывается соотношениями трехмерной теории упругости (без введения дополнительных упрощающих гипотез). Перемещения аппроксимируются степенными и лагранжевыми полиномами различных порядков, которые учитывают смещения МнКЭ как жесткого целого. Лагранжевые полиномы эффективно используются при проектировании МнКЭ оболочечного типа. Предлагаемые МнКЭ образуют дискретные модели малой размерности (в 103 ^ 106 раз меньше размерностей базовых моделей) и порождают приближенные решения, которые быстро сходятся к точным, что дает возможность строить при небольших временных затратах решения с малой погрешностью. Для верификации МнКЭ используется известный численный метод. Разработаны и численно исследованы трехсеточные конечные элементы (ТрКЭ) оболочечного типа. Приведен пример расчета многослойной оболочки с применением разработанных ТрКЭ и базовой модели, которая имеет около 1,4 миллиарда узловых неизвестных метода конечных элементов.

Ключевые слова: упругость, оболочки, композиты,

криволинейные композитные многосеточные конечные

элементы.

БОТ 10.14258Лгуа8и(2017)4-22

Procedures for developing curvilinear three-dimensional composite multigrid finite elements (MFE) of a shell-like type for calculating the stress state of elastic cylindrical shells having an inhomogeneous (microinhomogeneous) structure and static loading have been proposed. MFE are developed in local Cartesian coordinate systems on the basis of small (basic) shell partitions (models). When constructing MFE (without increasing their dimensionality), arbitrarily small basic shell partitions can be used. Thus, it is possible to take into account their inhomogeneous and microinhomogeneous structure, irregular shape, complex nature of loading and fastening within the micro-approach. The stress-strain state in MFE is described by the formulas of the three-dimensional theory of elasticity (without introducing any additional simplifying hypotheses). The displacements are approximated by power and Lagrange polynomials of various orders, which take into account the displacements of the MFE as a rigid whole. Lagrangian polynomials are effectively used while developing shell-type elements. The proposed MFE yield the small dimensional discrete models (103 ^ 106 times less than the dimensions of the reference models) and generate some approximate solutions that quickly converge to exact ones, which enable the construction of solutions with a high accuracy for a short time. A known numerical method is used to verify the MFE. Three-grid finite elements of a shelllike type have been developed and numerically studied. An example of a multilayer shell calculation using the developed three-grid finite elements and a reference model that has about 1.4 billion nodal unknowns of finite element method has been given.

Key words: elasticity, shells, composites, curvilinear composite multigrid finite elements.

Введение. В основе технических теорий композитных цилиндрических оболочек лежат гипотезы, накладывающие определенные ограничения на поля перемещений, деформаций и напряжений [1, 2]. Это порождает неустранимую погрешность в решениях и ограничивает области применений этих теорий. Метод конечных элементов (МКЭ) [3, с. 169; 4, с. 255] широко используется при исследовании напряженно-деформированного состояния упругих оболочек [5-8]. Построение криволинейных конечных элементов (КЭ) при расчете оболочек встречает различные трудности [5, с. 19], в частности, связанные с учетом смещений КЭ как жесткого целого. Расчет композитных цилиндрических оболочек с применением КЭ в постановке трехмерной задачи теории упругости [9, с. 40] с учетом их структуры приводит к построению базовых разбиений очень высокой размерности, для которых применение программ расчета ANSYS, МАЗТЯАК и др. [5, с. 7] затруднительно. В работах [10, 11] проведен расчет композитных цилиндрических панелей и оболочек с помощью МнКЭ, построенных с применением для аппроксимации перемещений только степенных полиномов 1-го, 2-го и 3-го порядков.

В данной работе изложены процедуры построения трехмерных композитных (однородных) многосеточных конечных элементов (МнКЭ) оболочечного типа с применением полиномов Лагранжа для расчета напряженного деформированного состояния линейно упругих композитных цилиндрических оболочек. Разработаны и численно исследованы трехсеточные оболочечные КЭ. Для верификации разработанных КЭ используется известный численный метод [4, с. 175]. Расчеты показывают, что реализация МКЭ для многосеточных дискретных моделей требует небольших временных затрат и в 103 ^ 106 раз меньше объема памяти ЭВМ, чем для базовых моделей.

1. Композитные криволинейные двухсеточные КЭ. Отметим, что в расчетах конструкций с применением МнКЭ используются последовательности базовых {Е°п} и многосеточных [Кп} дискретных моделей оболочки. Базовая модель К" оболочки состоит из однородных односеточных КЭ V, 1-го порядка с характерными размерами кехп хкеуп хк^ (рис. 1), где е — порядковый номер КЭ V, , п - порядковый номер базовой модели Е°п, а,п - угол раствора КЭ V, , ей— ось оболочки. Процедура построения КЭ V" в локальной декартовой системе координат Ох подробно изложена в работе [10]. Процедуру построения криволинейных композитных двухсеточных КЭ (ДвКЭ) с применением лагранжевых полиномов покажем на примере ДвКЭ Т'па 3-го порядка размерами 9кеш х9кеуп х9кет (рис. 2), расположенного в локальной декартовой системе координат О2х2у2г2, а — порядковый номер КЭ , п — порядковый номер многосеточной модели К . Пусть между компонентами неоднородной структуры ДвКЭ связи идеальны.

Рис. 1. Однородный КЭ VI

Рис. 2. ДвКЭ Тпа

Базовое разбиение КЭ состоит из КЭ V, (рис. 1), которые учитывают в ДвКЭ неоднородную и микронеоднородную структуры, сложный характер нагружения и закрепления. Функционал полной потенциальной энергии П" для базового разбиения КЭ ТУ" запишем в виде

е=1 2

(1)

где [К, ] — матрица жесткости; Р,, бе — векторы узловых сил и перемещений КЭ УД которые отвечают системе координат О2 Х2 уг г2; М — общее число КЭ V"; Т — транспонирование.

Рассмотрим построение полиномов Лагранжа в локальной криволинейной системе координат О2£пя на узловой сетке КЭ (рис.2). Пусть узел Р(гу],к) имеет координаты £, п^, Ск, на рисунке 2 I = j = 3 , к = 4. Заметим, что у = п , при малых углах раствора ДвКЭ можно считать, что х & £ , г & С . Принимаем

х = £, у = П, г = С .

(2)

Базисную функцию А.к для узла Р(1, к) в системе координат 02 х2 у2 г 2 с помощью полиномов Лагранжа Ь (х), Ь (у), Ьк (г) [4, с. 187] запишем в виде

(х, у, г) = Ь (х)Ь. (у)Ьк (г),

(3)

где = П

х — х_

п 1 .ПТ!

ш = п

у-у.

Ш= П

х1, у., гк — координаты узла Р(1,к) в системе координат 02 х2 у2 г2.

Для точки с координатой £ , лежащей на цилиндрической поверхности радиуса Я , имеем £ = аЯ , а — угол для координаты £ (рис. 2). Учитывая (2) и соотношения вида £ = аЯ , £ = аЯ в (3), получаем А.к (а,п,0 = Ь (а)Ь. (г))Ьк (О,где Ь (а), Ь (г,), Ьк (О — полиномы Лагранжа, имеющие следующий вид:

4(«)= П

а-а

40= п

ш= п „

у-ч,

Пп

- (4)

В расчетах удобно использовать полиномы Лагранжа (4). Функции перемещений иа, Уа, wa КЭ Ш", построенных с помощью полиномов Лагранжа (4), представим в форме

Замечание 1. В силу (6) размерность вектора 8а (т.е. размерность КЭ ) не зависит от общего числа М КЭ V", представляющего КЭ . Следовательно, можно использовать сколь угодно мелкие базовые разбиения, что позволяет учитывать неоднородную и микронеоднородную структуру в КЭ .

Пусть найден вектор 8а. По формулам (6) находим векторы 8е перемещений КЭ V" (е = 1,...,М ) в системе координат 02х2у2г2. Затем, используя векторы 8е, определяем векторы перемещений 8^ КЭ V" в локальной системе координат 01х1 у1г1 (см. рис. 1). Используя 8е, вычисляем напряжения в КЭ V", е = 1,...,М .

2. Композитные трехсеточные КЭ оболочечно-го типа. Процедуру построения композитных криволинейных трехсеточных КЭ (ТрКЭ) с применением полиномов Лагранжа рассмотрим на примере ТрКЭ О" оболочечного типа, расположенного в локальной декартовой системе координат 03 х3 у3 г3 и имеющего характерные размеры 81кех11 х 81кеуп х 9кехп (рис. 3), p — порядковый номер КЭ О"р , п — порядковый номер многосеточной модели Яп. Область КЭ О" состоит из ДвКЭ , где а = 1,...,N , в данном случае N = 81. Узлы крупных сеток КЭ образуют мелкую сетку КЭ О"р , на которой определяем крупную сетку Нр КЭ О" . На рисунке 3 узлы сетки Нр (32 узла) отмечены точками. Полную потенциальную энергию П"р КЭ О" запишем в виде

Д. Л

(9)

иа = Й Чф , 0 а = Й 4.0 , ^а = Й Ав^ , (5)

в=1 в=1 в=1

где qu, ц0, qw, N — перемещения и функция формы в -го узла узловой сетки КЭ , "0 = """ , в рассматриваемом случае "0 = 64 (рис. 2).

Порядок КЭ определяется порядком полинома, построенного на его узловой сетке. Используя (5), вектор узловых перемещений 8е КЭ V" выражаем через вектор узловых перемещений 8а КЭ в системе координат 02 х2 у2 г2

(6)

где [А"] - прямоугольная матрица, е = 1,...,М . Подставляя (6) в (1) и следуя принципу минимума полной потенциальной энергии для КЭ , т.е. 5Я"(8о)/б8а = 0, получаем соотношение [^Га]8а=Ро, которое отвечает равновесному состоянию КЭ Ш",

где [Ка ] - матрица жесткости и Pa , 8а - векторы узловых сил и перемещений КЭ , отвечающие системе координат 03х3у3г3, а = 1,...,N .

Рис. 3. Оболочечный КЭ О"

р

Функции перемещений ир

кэ ор„

по-

строенные на крупной сетке Нр с помощью полиномов Лагранжа, представим в форме

[Ка ] = Й[А" ]Г [Ке ][А ] , Fa ^ [Ае ]' Pe . (7)

е=1 е=1

где [Ка ] - матрица жесткости, Fa - вектор узловых сил КЭ Ш".

ир=Й Ав ь, ур=ЙЙ Ав ь, wp=Й Ав ь, (10)

в= в=1 в=1

где , , qW, - перемещения и функция формы в -го узла узловой сетки КЭ О"р , "0 = """, в рассматриваемом случае "0 = 32 (рис. 3).

р

Используя функции (10), узловые перемещения вектора 6а КЭ выражаем через узловые перемещения вектора Ьр крупной сетки КЭ Ср . В результате получим равенство

(11)

где [Ар] - прямоугольная матрица, а = 1,...,N .

Используя (11) в (9) и минимизируя функционал П"р (5р) по перемещениям бр, для КЭ С"р получаем соотношение которое отвечает его равновесному состоянию,

где [К, ] = £ [Ар ]г [К ] [Ар ], Бр = £ [Ар ]г Ра , (12)

а=1 а=1

[Кр ] - матрица жесткости, Бр - вектор узловых сил ТрКЭ впр.

Замечание 2. Согласно равенству (11) размерность вектора (т.е. размерность ТрКЭ Спр ) не зависит от общего числа ДвКЭ ^р и в силу (6) от общего числа базовых КЭ Ур, составляющих ТрКЭ. Следовательно, разбиение ТрКЭ С на ДвКЭ шр и односеточные КЭ Ур может быть сколь угодно мелким, что позволяет в рамках микроподхода учитывать сложную неоднородную структуру ТрКЭ.

Процедура определения напряжений в ТрКЭ аналогична процедуре определения напряжений в ДвКЭ. Используя ТрКЭ, по процедуре, аналогичной п. 2, строим 4-сеточные КЭ и т -сеточные КЭ [12], т > 4. Однородные МнКЭ проектируются по процедурам, аналогичным процедурам п.п. 1, 2.

3. Результаты численных экспериментов. Рассмотрим в декартовой системе координат Охуг модельную задачу деформирования 9-слойной цилиндрической оболочки У0 при локальном нагружении (рис. 4). Левый торец оболочки жестко закреплен, т.е. при у = 0 имеем и = V = ш = 0.

Внутренний радиус оболочки равен 200, толщина h = 9, длина L = 650. Слои (толщиной hc = h /9) являются однородными изотропными телами. Модули Юнга девяти слоев (начиная с нижнего) равны: 1, 3, 5, 10, 15, 20, 30, 40, 50. Коэффициент Пуассона равен 0,3. На внешней поверхности S оболочки длиной L /2 и с углом раствора а = п /2 (рис. 4), симметричной относительно плоскости Oyz, действует равномерно распределенная радиальная нагрузка q = 0,01. В расчетах используем половину оболочки. Базовая дискретная модель R (половины оболочки) состоит из однородных КЭ Ур 1-го порядка (рис. 1) с характерными размерами h'xnxh'pxh'zn, определяемыми по формулам

he = h". / р , he = he. / р , h" = h\ / р, р = 1,...,8 , (13)

xn x1 > ур y1 > zn z1 ' ' ' ' v '

где h'xl = а^, hey1 = L/324, h'tl = h/9, ae = п/324, R1e1 — размер для КЭ Ур при р = 1.

Сетка базовой дискретной модели R имеет размерность m1 хтр хтр,

где

т1 = 324р +1, т2 = 324р + 1,

т = 9р +1, р = 1,...,8 .

(14)

Рис. 4. Расчетная схема оболочки У

1 « 2 тр — размерность по круговой координате; тп —

по оси Оу, тр — по оси Ог .

Многосеточная дискретная модель Яп состоит из ТрКЭ Спр оболочечного типа с характерными размерами 81Н' х81Не х9Н' , рис. 3 (см. п. 2).

А хп уп гп 7 А 4 '

Трехсеточные элементы ор состоят из ДвКЭ с характерными размерами 9Н"ш х 9кеп х 9кегп (см. рис. 2, п. 1). Результаты расчетов даны в таблице 1, где шп, ап — максимальные прогиб и эквивалентное напряжение многосеточной дискретной модели Яп (напряжения определяются по четвертой теории прочности). Значения 6ап (%), 6шп (%) (п > 2) находим по формулам

6 (%) = 100% х | а — а 1 |/а ,

а ,п 4 ' I п п—11 п >

6 (%) = 100% х|ш — ш Л/ш . (15)

ш ' I п п—1 I п у '

Максимальные прогибы и эквивалентные напряжения моделей Rt¡

Таблица 1

R п R1 R2 R3 R4 R5 R6 r7 Rs

w п 148,62 196,91 208,64 213,48 215,90 217,27 218,14 218,71

6 (%) w р v J - 24,52 5,62 2,27 1,12 0,63 0,40 0,26

Яр 2,715 3,616 4,021 4,291 4,483 4,628 4,741 4,830

6 (%) арv J - 24,92 10,07 6,29 4,28 3,13 2,38 1,84

Характер изменения величин 6шп (%), 6ап (%) (п = 2,...,8) демонстрирует быструю сходимость напряжений ап и перемещений шп к точному решению. Поскольку величины 6ш 8 = 0,0026 , 6а 8 = 0,0184 малы (п = 8), то считаем, что перемещение ш8 = 218 ,71 и напряжение а8 = 4,83 определены с малой погрешностью. Максимальные эквивалентные напряжения ап возникают в девятом слое оболочки при х = 0, в окрестности ее крепления, максимальные перемещения шп — на свободном торце оболочки У0 .

Размерность базовой дискретной модели Я (для половины оболочки) равна 1471532832 (более 1,4 млрд. неизвестных), ширина ленты системы уравнений (СУ) МКЭ - 567873. Многосеточная модель Я8 имеет 249696 узловых неизвестных, ширина ленты СУ МКЭ равна 7863. Реализация МКЭ для многосеточной модели Я8 требует в 425 тыс. раз меньше объема памяти ЭВМ, чем для модели Я0. Построение решений для моделей {Я }8=1 на персональном однопроцессорном компьютере требует 4 часа 25 минут временных затрат.

4. Численный метод верификации КЭ. Для исследования сходимости приближенных решений, построенных с применением новых КЭ, широко используют известный численный метод [4, с. 175], краткая суть которого заключается в следующем. С помощью предлагаемых КЭ решается аналогичная (тестовая) задача с известным аналитическим (точным) решением и . Пусть ошибка 6к = || и0 — ик 10 при к ^ 0 , где h — характерный размер КЭ, иИ — решение тестовой задачи, построенное с помощью предлагаемых КЭ. Тогда считают, что приближенные решения, построенные с помощью предлагаемых КЭ и для дру-

гих конкретных задач, аналогичных тестовой задаче, сходятся в пределе (при к ^ 0) к точному решению.

Рассмотрим в качестве тестовой модельную осе-симметричную задачу о деформировании девятислой-ной цилиндрической оболочки У1, которая имеет такие же размеры, условия закрепления и модули упругости, как и оболочка У0 (см. п. 3). Пусть оболочка У1 находится под действием осесимметричной радиальной равномерной нагрузки q = 0,1 при у > Ь /2. Как известно, последовательность решений (при к ^ 0) осесимметричной задачи для данной оболочки, построенных по МКЭ с применением стандартных КЭ [3, с. 87], которые являются однородными кольцами с квадратным поперечным сечением, сходится к точному решению. Для дискретной модели Qn оболочки, п = 1,...,8 , состоящей из стандартных КЭ, определяем максимальные прогиб шап и эквивалентное напряжение а0п. Для полученных ш0п, а0п находим погрешности 80°п (%), 80ап (%) по формулам, аналогичным (15). Аналогичные параметры ш'п, а'п, 6'ш п (%), 6га п (%) для осесимметричной задачи определяем с помощью КЭ Опр (см. п. 3), используя для половины оболочки у законы измельчения (13), (14) базовых моделей. Значения погрешностей ¿>ст°я(%), ¿£„(%) показывают быструю сходи-

мость напряжений а0п, агп и перемещений ш0п, ш'п к точным решениям (см. табл. 2, 3). В силу малости величин 6° 8 = 0,000000044 , 6° 8 = 0,000087 ,

ш, 8 у а, 8 у

6гт ,8 = 0,0000017, 6а,8 = 0,0000439 прогибы и напряжения при п = 8 считаем точным решением. Для стандартных КЭ получаем ш^ = 25,775620, а80 = 6,087667, для ТрКЭ — ш8 = 25,777645 , а\ = 6,087704.

Таблица 2

Максимальные прогибы и эквивалентные напряжения моделей Qn

п Размерность модели Qn 6° (%) ш ,п 4 ' 0 ап 6а, п (%)

2 1297х19 25,775620842 0,0023632 6,0766437388 0,2352582

4 2593х37 25,775488560 0,0001529 6,0839633054 0,0403152

6 3889х55 25,775452136 0,0000553 6,0864290250 0,0162303

8 5185х73 25,775440899 0,0000044 6,0876676879 0,0087283

Таблица 3

Максимальные прогибы и эквивалентные напряжения многосеточных моделей, состоящих из ТрКЭ

п шп 6ш ,п (%) Г ап 6а, (%)

2 25,974091346 1,5732658 6,1216614786 1,0042443

4 25,796347946 0,1717619 6,0860925957 0,0414311

6 25,781094987 0,0052265 6,0862833577 0,0185614

8 25,777645991 0,0001745 6,0877048793 0,0043964

Относительные погрешности

6w (%) = 100% х| w80 - wr8 ,

(%) = 100% х\ - ага \/а°а

соответственно равны 0,00061 и 0,00132%.

На основе полученных результатов можно сделать вывод, что предлагаемые ТрКЭ О" порождают решения а", wrр, которые в пределе (при " ^ ж) стремятся к точному решению данной осесимметричной задачи. Тогда, следуя методу верификации, можно

считать, что приближенные решения, построенные для задачи п. 3 с применением ТрКЭ О", при " ^ ж сходятся к точному решению.

Заключение. Предложены многосеточные конечные элементы оболочечного типа для расчета напряженного состояния упругих композитных цилиндрических оболочек. Разработаны и численно исследованы трехсеточные КЭ оболочечного типа. Пример расчета многослойной оболочки с применением предложенных КЭ показывает высокую эффективность их применения.

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

1. Болотин В.В., Новиков Ю.Н. Механика многослойных конструкций. — М., 1980.

2. Голушко С.К., Немировский Ю.В. Прямые и обратные задачи механики упругих композитных пластин и оболочек вращения. — М., 2008.

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

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

5. Голованов А.И., Тюленева О.И., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. — М., 2006.

6. Клочков Ю.В., Николаев А.П., Шубович А.А. Анализ напряженно-деформированного состояния оболочек вращения в геометрически нелинейной постановке при различных вариантах интерполяции перемещений. — Волгоград, 2013.

7. Киселев А.П. Расчет тонких оболочек на прочность в трехмерной постановке без упрощающих гипотез // Известия вузов. Строительство. — 2008. — № 1.

8. Киселев А.П., Гуреева Н.А., Киселева Р.З. Расчет многослойных оболочек вращения и пластин с использованием объемных конечных элементов // Известия вузов. Строительство. — 2010. — № 1.

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

10. Матвеев А.Д., Гришанов А.Н. Одно- и двухсеточ-ное криволинейные элементы трехмерных цилиндрических панелей и оболочек // Известия Алтайского гос. унта. — 2014. — № 1/1. DOI:10.14258/izvasu(2014)1.1-19.

11. Матвеев А.Д., Гришанов А.Н. Многосеточные криволинейные элементы в трехмерном анализе цилиндрических композитных панелей с полостями и отверстиями // Ученые записки Казанского гос. ун-та. — 2014. — Т. 156, кн. 4.

12. Матвеев А.Д. Построение сложных многосеточных элементов с микронеоднородной структурой // Численные методы решения задач теории упругости и пластичности : тезисы докладов XXIII Всероссийской конференции (Барнаул, 26-28 июня 2013 г). — Новосибирск, 2013.