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

Вычислительные технологии

Том 11, № 4, 2006

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

Ю. В. Клочков, А. П. Николаев, А. Ш. Джабраилов Волгоградская государственная сельскохозяйственная академия, Россия

e-mail: [email protected]

A comparative analysis of the efficiency of use of isoparametrical finite elements and elements which geometrical characteristics are calculated using exact formulas is presented for the case of shells allowing rigid displacements under the influence of loading. Triangular and quadrangular finite elements are considered. Presented calculations for the case of an ellipse of revolution and a cylinder leaning on a spring support, have allowed to draw a conclusion on necessity of use of the finite elements which geometrical parameters are calculated using analytical formulas combined with vector interpolation of translations.

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

1. Алгоритм формирования матрицы жесткости конечного элемента

При формировании матрицы жесткости конечного элемента и столбца внешней нагрузки использовалось равенство работ внешних и внутренних сил на возможном перемещении [1]:

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

Введение

где }Т = (єііє22єі2І — матрица-строка деформаций в произвольном слое оболочки,

отстоящем от срединной поверхности на расстоянии £;

а

{аав} = { а

її

22

— матрица-столбец напряжений в произвольном слое оболочки;

а

12

{и} = (и у '} — матрица-строка, содержащая компоненты вектора перемещения

точки срединной поверхности оболочки;

( Р \

(Р} = < Р2 > — матрица-столбец внешней нагрузки.

I рз

Интеграл по объему конечного элемента представляет собой работу внутренних сил, а интеграл по площади поверхности конечного элемента (так как внешняя нагрузка является поверхностной) — работу внешних сил на возможном перемещении.

Применяя традиционную интерполяционную процедуру [1], столбец (и} можно представить в виде матричного произведения

{и} = [Л]{ЦЛ},

(2)

где {ил}т = {{и^}^^}^^}} — столбец узловых варьируемых параметров конечного элемента в локальной системе координат, а матрица [А], содержащая аппроксимирующие функции <^, имеет следующий вид:

[А]

мг {0}т {0}т

{0}т Мг {0}т

{0}т {0}т мт

Используя геометрические и физические соотношения теории оболочек [2], можно сформировать матричные зависимости

{єав }

[Г]{Єав}, {аав} = [С]{4,в},

Єав = [С][А]{и?},

(3)

где {бав}Т верхности;

{єіі є22 Є12 Хп Х22 Х12} — столбец деформации точки срединной по-

[Г]

3x6

(4)

100 с о о 0100 с о 0 0 2 0 0 2(

Входящие в (3) матрица упругости [С] и матрица дифференциальных операторов [Д имеют следующий вид:

Е Ей

[С ]

3x3

1 — V2 1 — V2

Ей Е

1 — V2 1 — V

2

Е

2(1 + V)

(5)

0

0

0

0

Здесь Б и в — глобальные криволинейные координаты точки срединной поверхности оболочки вращения (Б — длина дуги меридиана; в — угол, отсчитываемый от образующей против хода часовой стрелки); г — радиус вращения; х — глобальная осевая координата. С учетом (2)-(6) функционал (1) примет вид

{и;}т I [в]т [Г]т [С ив ] лу {и;}т = {и;}т | [А]т (р } лр. (7)

Выполняя операцию минимизации функционала (1), можно записать следующее матричное уравнение:

[к ]{и;} = {Я}, (8)

где [К] = f[B]т[Г]т[С][Г][В] ЛУ — матрица жесткости конечного элемента в локальной

V

системе координат; {Я} = /[А]т(Р} — столбец внешней нагрузки конечного элемента.

Матрицу жесткости и столбец внешней нагрузки конечного элемента в глобальной системе координат получают в результате умножения [К] и {Я} на матрицу преобразований:

[К ]г = [Т ]т [К ] [Т ], (Я}г = [Т ]т {Я}, (9)

где элементы матрицы [Т] определяются на основании соотношений

ддш = ддш дБ + ддт дв

~о£ = + ~ЖЖ’

ддш = ддт дБ ддт дв

дп дБ дп + дв дг/'

Здесь под д понимается компонента вектора перемещения и, V или '; £ и п — локальные координаты конечного элемента, а верхний индекс т указывает на номер узла используемого дискретного элемента.

2. Аппроксимация перемещений и геометрических величин

В настоящей статье рассматриваются тре- и четырехугольные конечные элементы, узловыми варьируемыми параметрами которых выбираются компоненты вектора перемещения и их первые производные. Столбы неизвестных в узле КЭ в локальной и глобальной системах координат имеют вид

{иул}т = (и и? и,^ V v>v ' };

1x9

(иу}Т = (и и,5 и>0 V V)0 ' }, (11)

1x9

где запятая означает операцию дифференцирования по соответствующей координате.

Перемещение внутренней точки КЭ определяется через свои узловые значения с помощью интерполяционной зависимости следующего вида:

д = Мт № (12)

где матрица-строка {^}т содержит функции формы.

Структура матрицы-строки {^}т и столбца узловых значений компонент вектора перемещения в локальной системе координат дл зависит от типа используемого дискретного элемента, например, для треугольных КЭ с девятью степенями свободы в узле можно записать

Мт = {^2 -"Ы,

1x9

(дЛ}т = (дУ дк д?« й д,« д?п дП дкп}, (13)

1x9

где ^1 - - - ^9 получены с использованием полных двумерных полиномов третьей степени [3]; г,], к — узлы треугольного КЭ; локальные координаты £ и п изменяются в пределах 0 < £, П < 1.

При использовании в качестве элемента дискретизации четырехугольного КЭ столбцы {^} и дл принимают следующий вид:

Мт = М ^2 '-'^12},

1x12

12/

1x12

(14)

где <^1... ^12 представляют собой произведения полиномов Эрмита третьей степени, а локальные координаты £ и п изменяются в диапазоне —1 < £, п < 1.

При изопараметрической аппроксимации геометрические величины внутренней точки КЭ, входящие в матрицу [Д соотношения (6), определяются через свои узловые значения по формулам, аналогичным (12). Так, например, для радиуса вращения можно записать интерполяционную зависимость

г = мт к}, (15)

где структура столбца |г^} подобна структурам {д^} из соотношений (13) и (14).

Столбец {г^} может быть выражен через столбец узловых значений радиуса вращения и его производных в глобальной системе координат

К} = [РД]{гГ},

(16)

где {Гу/ = {гггЙг г^г^г^г^} для треугольного конечного элемента;

{гГ}т = {гггЙгкг1 г^г^У^гЙвгквг|в} для четырехугольного КЭ.

Матрица преобразования [РД] для тре- и четырехугольного конечных элементов имеет размеры 9 х 9 и 12 х 12 соответственно.

С целью сокращения объема статьи ниже приводятся соотношения, справедливые только для треугольного конечного элемента.

Если в качестве глобальных координат а и в использовать длину дуги меридиана Б и окружную координату в, то матрица РД для треугольного конечного элемента будет иметь следующий вид:

[РД]

9x9

1

1

1

-и^ Б •и^ в

Б,с •и^

•и^ Б •и^

Б,п в,п

Б,п в,п

Б,п в

(17)

Глобальные координаты Б и в внутренней точки треугольного КЭ выражаются через координаты узлов соотношениями

Б =(1 — £ — п)Бг + £БЙ + пБк, в = (1 — £ — п)в" + £вЙ + п'к, (18)

где £, п изменяются в пределах 0 < £, п < 1.

Входящие в (17) производные Б?, Б,^, в,£ и в,^ определяются дифференцированием (18):

Б,? = БЙ — Б", Б,п = Бк — Б", в,? = вЙ — в", в,п = вк — вк. (19)

Зависимость (15) с учетом (16) примет вид

г = Мт [РД]{г£}. (20)

Первые, вторые и третьи производные радиуса вращения в точке интегрирования определяются соответствующим дифференцированием (20) по глобальной координате Б:

г,- = ({^5ГО + Ып}ТП,)

г,-- = ({^,55}Т(Ы2 + {^,пп}Т(П,-)2 + 2{^,5п}Т£,-П,-) [рД]{гГ}, г,--- = ({^,555 }ТС3- + {^,55п }ТС2-П,- + {^,пп5}Т П,2£,- + Ыппп}Т П3-+

+2С-П,-({^,5п5}Т^,- + {^,5пп}ТП,-)) [Р^]{гГ}- (21)

С учетом (20) и (21) формируется матрица жесткости изопараметрического треугольного конечного элемента размером 27 х 27. Данный тип элементов использовался для дискретизации оболочек, допускающих жесткие смещения под действием заданной нагрузки.

При этом в названном элементе геометрические параметры узловых значений векторов

г

гу определяются по аналитическим зависимостям рассматриваемой поверхности, а уже в каждой точке интегрирования с локальными координатами £, п геометрические параметры определяются аппроксимирующими соотношениями (21). Операция использования соотношений (21) при формировании матрицы жесткости конечного элемента приводит к значительным погрешностям.

Координаты точки интегрирования £ и п в локальной системе координат являются заданными величинами. По заданным координатам точнее определять с использованием аппроксимирующих соотношений (18) глобальные координаты Б и в, после чего по аналитическим зависимостям рассматриваемой поверхности находить геометрические параметры в точке интегрирования конечного элемента, отказавшись от названия “изопара-метрический элемент”.

Результаты конечно-элементных решений, полученных при использовании в качестве элементов дискретизации изопараметрических конечных элементов, сравнивались с результатами численного анализа, выполненного с применением конечных элементов, в которых геометрические величины в точке интегрирования по площади конечного элемента в локальной системе координат вычислялись по точным аналитическим формулам после определения глобальных координат Б, в по соотношениям (18).

Пример 1. Рассчитан фрагмент эллипсоида вращения, нагруженный внутренним давлением интенсивности д (рис. 1). Исходные данные имели следующие значения: Е = 2 • 105 МПа, ^ = 0.3, д = 5 МПа, і = 0.02 м, параметры эллипса а =1.3 м, Ь = 0.9 м, координата х изменялась в пределах 0 < х < 1.2 м. На правом краю оболочка имеет пружинные опоры, позволяющие ей смещаться в меридиональном направлении как абсолютно твердому телу под действием заданной нагрузки. Вследствие наличия осевой симметрии эллипсоид представлен одной полоской треугольных КЭ, ориентированной в меридиональном направлении.

Расчеты выполнены в трех вариантах. В первом варианте для дискретизации оболочки использовались изопараметрические конечные элементы; во втором и третьем вариантах геометрические характеристики в точке интегрирования вычислялись по аналитическим формулам, определяющим форму срединной поверхности рассчитываемой оболочки. Второй и третий варианты различались способом аппроксимации перемещений внутренней области треугольного КЭ. Во втором варианте использовалась традиционная интерполяционная процедура (12). В третьем варианте реализован векторный способ аппроксимации перемещений в треугольном КЭ [4].

Рис. 1.

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

Анализ результатов, представленных в табл. 1, показывает, что во всех трех вариантах наблюдается сходимость вычислительного процесса к точному решению, полученному по формулам сопротивления материалов, однако скорость сходимости вычислительного процесса в первом варианте ниже, чем во втором и третьем.

Если жесткость пружинных опор на правом краю оболочки постепенно уменьшать, то оболочка получит возможность смещаться в меридиональном направлении как абсолютно жесткое тело. Причем усилие пружины прямо пропорционально величине смещения оболочки как жесткого целого.

Результаты расчета оболочки при различных значениях величин смещения оболочки как жесткого целого представлены в табл. 2, где приведены меридиональные и кольцевые напряжения в концевом сечении оболочки в точках 1, 2 (рис. 1) рассчитываемой конструкции (при сетке узлов 2 х 33).

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

По результатам анализа приведенного табличного материала можно сделать однозначный вывод, что вычисление геометрических параметров срединной поверхности оболочки в каждой точке интегрирования по аналитическим функциям в зависимости от глобальных координат, найденных по соотношениям (18), гораздо проще и эффективнее изопара-метрической параметризации (15).

Таблица 1. Количество элементов дискретизации для трех вариантов расчета

х, м а, МПа Вариант I Вариант II Вариант III Точное решение [5]

8 16 32 48 8 16 32 48 8 16 32 48

0.00 аМ 96.08 95.92 95.88 95.87 95.82 95.85 95.86 95.86 95.94 95.89 95.88 95,уу87 95,уу86

ак 179.17 179.08 179.09 179.09 179.07 179.06 179.08 179.09 179.03 179.05 179.08 179.09 179.06

1.20 аМ 10.99 3.68 1.11 0.53 2.37 0.71 0.16 0.06 3.64 0.90 0.21 0.09 0.00

ак 174.0 169.51 168.05 167.73 170.90 168.82 167.84 167.63 169.70 168.70 167.86 167.65 167.82

Таблица 2. Величина смешений оболочки как жесткого целого, для трех вариантов расчета, м

х, м 0, рад а, МПа Вариант I Вариант II и р а В ант III

0.00 0.09 0.90 9.00 0.00 0.09 0.90 9.00 0.00 0.09 0.90 9.00

аМ 1.11 512.64 5108.0 50331.1 0.16 -3.02 -31.68 -318.0 0.21 0.21 0.21 0.21

0.00 ак 168.05 290.19 1387.6 12185.4 167.84 117.53 -335.34 -4858.6 167.86 167.86 167.86 4167.85

1.20 п аМ 1.06 553.01 5512.4 54308.2 0.10 -2.19 -22.814 -228.8 0.12 0.12 0.12 0.12

60 ак 168.04 301.66 1502.2 13315.0 167.83 117.78 -332.7 -4832.8 167.86 167.85 167.85 167.84

К вопросу о неадекватности изопараметрической параметризации ... 61

Пример 2. Решена задача по определению напряженно-деформированного состояния цилиндра, загруженного в середине сосредоточенной силой, опирающегося на пружинную опору (рис. 2). Исходные данные: Е = 7.38• 104 МПа, ^ = 0.3125, Р = 45.36 Н, £ = 0.24 см, Я = 12.58 см, Ь = 26.26 см. Вследствие наличия плоскостей симметрии рассчитывалась четвертая часть цилиндра. Расчеты, как и в примере 1, выполнялись в трех вариантах, причем в качестве элемента дискретизации использовался четырехугольный конечный элемент.

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

Рис. 2.

Таблица 3. Результаты расчетов при жесткости пружины, равной бесконечности

Сетка элементов дискретизации Вариант I Вариант II Вариант III

5х 5 0.1612x2.54 0.1612x2.54 -0.0157x2.54

7x7 -0.02076x2.54 -0.02076x2.54 -0.0198x2.54

9x9 -0.02207x2.54 -0.02207x2.54 -0.0213x2.54

11x11 -0.02244x2.54 -0.02244x2.54 -0.0220x2.54

антах расчета наблюдается устойчивая сходимость вычислительного процесса. При этом значения прогиба в первых двух вариантах полностью совпадают, а в третьем величина нормального перемещения отличается в меньшую сторону всего на 2%.

Если жесткость пружины постепенно уменьшать, то оболочка получит возможность смещаться вертикально вниз как абсолютно твердое тело. Моделирование пружинной опоры можно осуществить следующим образом. Если жесткость пружины принять равной бесконечности, то данный факт будет соответствовать запрещению вертикального перемещения оболочки, что легко выполняется путем зануления соответствующего граничного условия. Так, например, для узла 2 (рис. 2) в матрице-строке граничных условий в колонке 7 табл. 4 стоит “нуль”, указывающий на то, что вертикальное перемещение данного узла запрещено. Цифра 1 означает, что перемещение или его производные разрешены.

Чтобы оболочка получила возможность смещаться как жесткое целое в вертикальном направлении, необходимо разрешить данное перемещение и одновременно увеличить соответствующий элемент матрицы жесткости КЭ на величину, равную жесткости пружины. Так, например, если необходимо задать оболочке жесткое смещение, равное единице, то величина жесткости пружины должна быть равна величине сосредоточенной силы, под действием которой цилиндр получит вертикальное единичное перемещение. Номер ячейки матрицы жесткости КЭ, в которую необходимо сделать соответствующую добавку, определяется физическим смыслом коэффициента матрицы жесткости КЭ. В рассматриваемом примере это ячейка под седьмым номером строки и столбца. Таким образом, в указанную ячейку матрицы жесткости КЭ необходимо добавить слагаемое, характеризующее жесткость пружины:

Я77 = Я77 + с,

где Я7,7 — коэффициент матрицы жесткости КЭ; с — жесткость пружины.

Результаты расчета цилиндра при различных значениях жесткости пружины представлены в табл. 5, где приведены значения кольцевых напряжений в наружных волокнах оболочки в точке 3 при сетке узлов 11 x 11 (рис. 2). Анализ данных таблицы показывает, что результаты расчета по вариантам существенно различаются между собой. При уменьшении жестокости пружины погрешности расчета в первом и во втором вариантах резко возрастают до совершенно неприемлемых значений. В третьем варианте наблюда-

Таблица 4. Матрица-строка граничных условий

и и,8 и,0 V V,.? w W)s w,в

1 2 3 4 5 6 7 8 9

0 1 0 0 0 1 0 0 0

Таблица 5. Параметры напряженно-деформированного состояния для трех вариантов расчета, Н/см2

Жесткость пружины Смещение оболочки как жесткого целого, см Вариант I Вариант II Вариант III

то 0 -608.1 -608.1 -623.8

43.36 1.0x2.54 -576.5 -576.5 -623.8

4.366 10.0x2.54 -339.6 -339.6 -623.8

0.4336 100x2.54 47.7 47.4 -623.4

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

Анализируя результаты расчетов, представленные в примерах 1 и 2, можно сделать следующие выводы.

1. Изопараметрическая форма параметризации срединной поверхности не позволяет в какой-либо мере учесть смещения КЭ как жесткого целого.

2. Для корректного учета жестких смещений в конечно-элементном анализе оболочек необходимо использовать векторную аппроксимацию перемещений [4].

Список литературы

[1] Постнов В.А., Хархурим И.Я. Метод конечных элементов в расчетах судовых конструкций. Л.: Судостроение, 1974. 344 с.

[2] Новожилов В.В. Теория тонких оболочек. Л.: Судпромиздат, 1962. 432 с.

[3] Клочков Ю.В., Николаев А.П., Киселев А.П. О функциях формы в алгоритме формирования матриц жесткости треугольных конечных элементов // Изв. вузов. Сер. Строительство. 1999. № 10. С. 23-27.

[4] Клочков Ю.В., Николаев А.П., Гуреева Н.А. Сравнение различных способов аппроксимации перемещений на треугольном элементе в расчетах оболочек // Вычисл. технологии. 2005. Т. 10, № 3. С. 47-55.

[5] ГоловАнов А.И. Новый конечный элемент для расчета произвольных тонких оболочек // Строительная механика и расчет сооружений. 1986. № 4. С. 21-23.

[6] Скоплинский В.Н. Об особенностях напряженно-деформированного состояния в области пересечения цилиндрических оболочек // Строительная механика и расчет сооружений. 1986. № 2. С. 19-22.

[7] Беляев Н.М. Сопротивление материалов. М.: Наука, 1976. 608 с.

Поступила в редакцию 22 февраля 2006 г., в переработанном виде —17 апреля 2006 г.