CC BY
31
Вычислительные технологии
Том 10, № 3, 2005
СРАВНЕНИЕ РАЗЛИЧНЫХ СПОСОБОВ АППРОКСИМАЦИИ ПЕРЕМЕЩЕНИЙ НА ТРЕУГОЛЬНОМ ЭЛЕМЕНТЕ В РАСЧЕТАХ ОБОЛОЧЕК
Ю. В. Клочков, А. П. Николаев, Н.А. Гуреева Волгоградская государственная сельскохозяйственная академия, Россия
e-mail: [email protected]
Vector interpolation of displacements was used to obtain the interpolation dependencies for a triangular finite element with nine degrees of freedom in the node. These interpolations express components of the vector of displacement of an internal point of the discrete element along with their derivatives through the column of varying parameters, which include values of all components of the node displacement vectors and their derivatives. Calculations of shells with significant gradients of curvatures of the middle surface or shells allowing stiff displacements have revealed principle advantages of the proposed method of the vector interpolation of displacements in the triangular finite element in contrast with generally accepted interpolation procedure.
Введение
В настоящее время в конечно-элементном анализе оболочек широкое применение находят дискретные элементы, узловыми неизвестными которых выбираются компоненты вектора перемещения и их первые производные по криволинейным координатам [1-3]. Это обстоятельство вызвано, по всей видимости, тем, что конечные элементы (КЭ), число узловых варьируемых параметров которых ограничено компонентами вектора перемещения и только их первыми производными, имеют существенно меньший (по сравнению с высокоточными КЭ) размер матрицы жесткости, что в свою очередь позволяет снизить требования к объему оперативной памяти используемого компьютера. Кроме того, ряд исследователей [3, 4] указывают на возможные проблемы совместности, которые могут возникнуть с применением в качестве элементов дискретизации высокоточных КЭ с набором узловых неизвестных, включающим производные высших порядков. Однако при использовании теории оболочек [5, 6] для определения напряженно-деформированного состояния требуются вторые производные нормальной компоненты вектора перемещения срединной поверхности. Поэтому для достижения приемлемых точности и сходимости вычислительного процесса требуется достаточно мелкая сетка элементов со столбцом узловых параметров, включающим только первые производные компонент вектора перемещения. Но в ряде случаев (когда градиенты кривизны срединной поверхности достигают существенных значений или оболочка имеет закрепления, допускающие ее смещения как абсолютно твердого
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.
тела) даже весьма густая сетка элементов с первыми производными перемещений (КЭ малого порядка) не позволяет получить удовлетворительный результат. Увеличение числа элементов дискретизации до бесконечности также невозможно в связи с накоплением ошибки округления, которая пропорциональна количеству КЭ.
Кардинально повысить эффективность применения такого рода КЭ к расчету оболочек можно за счет использования способа векторной интерполяции перемещений [7], принципиально отличного от общепринятой интерполяционной процедуры в методе конечных элементов. В настоящей статье на примере треугольного конечного элемента, узловыми неизвестными которого выбраны компоненты вектора перемещения и только их первые производные, показана высокая эффективность векторной интерполяции перемещений при расчете оболочек сложной конфигурации или допускающих в процессе своей эксплуатации смещения как жесткого тела.
1. Варианты способов аппроксимации полей перемещений
В качестве элемента дискретизации выбирается треугольный фрагмент срединной поверхности с узлами г,], к, который для удобства численного интегрирования отображается на прямоугольный треугольник с катетами единичной длины. Столбец узловых варьируемых параметров данного треугольного КЭ в локальной (0 < х, у < 1) системе координат выбирается в виде
кг= ( к-г к-г кг
1x27 ^ 1x9 1x9 1x9
где
ЫТ = {?У Яг,уЯ3,у9,у} •
1x9
Здесь и ниже под з понимаются тангенциальные г'1,!'2 или нормальная г-компоненты вектора перемещения. Согласно общепринятой интерполяционной процедуре [3] каждая компонента вектора перемещения внутренней точки треугольного КЭ интерполируется через узловые значения соответствующей компоненты
г1 = МТ К}, г2 = Мт К}, г = МГ {гу}, (2)
1x9 9^ 1x9 9^ 1x9 9x1
где {^}Т = {<^1, <^2, •••, ^9} — функции формы, полученные на основе полного полинома 1x9
третьей степени [8].
Для вычисления вторых производных нормальной компоненты вектора перемещения необходимо дважды выполнить дифференцирование последнего равенства (2) по глобальным криволинейным координатам а и в:
г,аа {^,аа} {гу} , г,ДО {^,ДО} {гу} , г,ав {^,ав} {гу} • (3)
Согласно формулам (3) вторая производная в каком-либо узле элемента будет равна алгебраической сумме первых производных всех узлов элемента.
В отличие от общепринятой интерполяционной процедуры, векторный способ аппроксимации предполагает варьирование не скалярной величины — отдельной компоненты
вектора перемещения, а непосредственно вектора перемещения точки срединном поверхности
V = МТ {Щ}, (4)
1x9 9x1
где
{®;}т = {«'«‘гХйХгкЙ^г';,}.
Столбец {V;} выражается через столбец векторных узловых неизвестных {й^} в глобальной системе координат соотношением
Jr'
Jy) L^J ll.
9x1 9x9 9x1
Ю = [L]KI • (5)
Здесь [L] — матрица вида
[L] =
1
1
1
да дх дв дх
да дх дв дх
да дх дв дх
да ду дв ду
да ду дв ду
да ду дв ду
Вектор перемещения и его производные в точке глобальной системы координат могут быть представлены в локальном базисе
v = v1a1 + v2a2 + va1,
—* 1 —* 2 —* —* v,a = taa1 + taa2 + taa1,
12 v,e = t^a1 + t@a2 + t@a,
12 v,aa taaa1 + taaa2 + taaa;
12 v,ee = teea1 + teea2 + tввa,
v,ae = taea1 + ^ (6)
где многочлены ta, ta, ta, te, te, te содержат компоненты вектора перемещения и их первые производные, а taa, t2aa, taa, t^g, t^g, tpp, t^, t2ap, taf3 включают и вторые производные компонент вектора перемещения.
Столбец {V,} с учетом (6) может быть представлен матричным произведением
К}
9x1
А
9x27
{ту }, 27x1
(7)
где {ту}Т = {Л2*^V2'VV1*V2*V* <*г^Цча*0"$$ ...4} ;
у
1x27
да
А
матрица ви-
9x27
от
{0'}
{а* }
ОТ
{а}
\ак }
ОТ
{0 }
К }
Входящая в структуру матрицы А строку, содержащую узловые векторы локального базиса
{о™} = {<0^0™} ,
подматрица {о™} представляет собой матрицу-
1x3
где индекс т обозначает узел треугольного конечного элемента %, ], к.
Столбец {ту} может быть выражен через стандартный столбец узловых варьируемых параметров в глобальной системе координат
{ту} = [Я] {Цуг},
27x1 27x27
(8)
где {иУ}Т = {КгГ КТ КГ \; МГ = чк^$} ■
I Т
у
1x27
^У 1x9
■Уу
1x9
,Т( ГЛ Т
V,
у
1x9
ГТ у
1x9
Вектор перемещения внутренней точки срединной поверхности элемента (4) с учетом соотношений (5) и (8) запишется в виде
V — {р}т [I] А [Я] {Цуг} = МТ А [С] [Я] {Цуг}
9x9 9x27 27x27
9x27
27x27 27x27
Здесь матрица [С] сформирована так, чтобы выполнялось равенство [I]
9x9
А
9x27
А
9x27
(9)
[С],
27x27
что достигается без принципиальных затруднений.
Векторы Ао1 и Ао2 локального базиса произвольной точки срединной поверхности элемента являются касательными к срединной поверхности, а вектор Ао направлен по нормали к ней. Векторы локальных базисов узловых точек А0т, 0^ и о™ выражаются через векторы локального базиса внутренней точки КЭ соотношениями
&2 — Сц&х + сі2^2 +
а — + ^і2^2 + ^І^а,
а1 — Ь3паі + 612^2 + Ь3иа и т. д.
(10)
С учетом (10) матрица А может быть представлена матричной суммой
А — [аі [Аі] + а2 [А2] + а [А3]] •
(11)
Принимая во внимание (6) и (11), соотношение (9) запишем как
ї]1аі + V2а2 + ьа — (аі [Аі] + сі2 [А2] + а [А3]) [С] [Н] {иГ} • (12)
Из равенства (12) можно записать искомые интерполяционные выражения, например
В результате дифференцирования (4) по криволинейным координатам можно получить первые и вторые производные вектора перемещения, например
гочленов 1^, Ь'О,, ■■■, *аа, , *ав, в состав которых входят искомые первые и вторые произ-
водные компонент вектора перемещения, например
Таким образом, можно отметить, что при векторном способе интерполяции перемещений каждая компонента вектора перемещения зависит от полного набора узловых варьируемых параметров, в состав которого входят узловые значения всех его компонент. Вычислительный алгоритм также претерпевает существенные и принципиальные изменения. Вместо одного, по сути, оператора, реализующего общепринятую интерполяционную процедуру, в подпрограмму формирования матрицы жесткости КЭ включается целый блок векторной интерполяции перемещений. Причем данный блок присутствует как на этапе формирования матрицы жесткости КЭ, так и на этапе вычисления напряжений (для вычисления вторых производных нормальной компоненты вектора перемещения), т. е. дважды.
2. Численные расчеты
Высокая эффективность применения векторной интерполяции перемещений при использовании КЭ низкого порядка подтверждается решением ряда конкретных задач.
Пример 1. В качестве примера была решена задача по определению напряженно-деформированного состояния компенсатора, находящегося под воздействием внутреннего давления интенсивности д (рис. 1). Радиус вращения задавался в виде
(13)
(14)
Из равенств (14) с учетом (5)—(11) можно получить интерполяционные формулы для мно-
(15)
г — А + В ооб(х/с).
г
Рис. 1.
Вследствие наличия осевой симметрии для дискретизации оболочки использовалась одна полоска элементов, ориентированная вдоль меридиана. Первоначально были приняты следующие исходные данные: (I = 0.2 МПа; А =1.3 м; В = 0.4 м; С = 0.48 м; £ = 0.01 м; Е = 2 • 105 МПа; V = 0.3. Координата х изменялась в пределах 0 < х < 1.44п м.
Расчеты были выполнены в двух вариантах: в первом варианте при формировании матрицы жесткости треугольного конечного элемента размером 27 х 27 использовалась традиционная независимая интерполяционная процедура (2), (3); во втором варианте был реализован предложенный способ векторной интерполяции перемещений (13), (15). Результаты расчета представлены в табл. 1, где приведены величины кольцевых ак и меридиональных ам напряжений (в мегапаскалях) во внутренних волокнах оболочки (для точек 1, 2 рассчитываемой конструкции) в зависимости от густоты сетки дискретизации. Анализ результатов расчета показывает, что и в первом, и во втором варианте расчета наблюдается быстрая сходимость вычислительного процесса. Значения расчетных величин в точках 1, 2 компенсатора практически совпадают в обоих вариантах расчета (что и следовало ожидать при решении осесимметричной задачи).
Выбранная расчетная схема, при которой край оболочки (точки 1,2 на рис. 1) остается свободным, позволяет сравнить полученное численное значение меридионального напряжения с точным, которое по физическим соображениям должно быть равно нулю. Из табл. 1 видно, что в обоих вариантах расчета контролируемое напряжение ам близко к нулю.
Если частоту волн компенсатора увеличить, например, в 2.67 раза (с = 0.18 м, 0 < х < 0.54п м), то градиенты кривизны срединной поверхности компенсатора в меридиональном направлении существенно возрастут. Результаты расчета оболочки при вышеупомянутом изменении исходных данных представлены в табл. 2. Анализ конечно-элементных решений показывает, что результаты расчета по вариантам принципиально различаются (в отличие
Таблица1
Точка Напряжение Вариант I Вариант II
2 х 13 2 х 25 2 х 49 2 х 13 2 х 25 2 х 49
1 ^к 18.20 17.51 17.52 17.35 17.56 17.52
^м 1.008 -0.042 -0.080 0.472 0.095 -0.067
2 ^к 18.20 17.51 17.52 17.33 17.55 17.52
^м 1.008 -0.042 -0.080 0.466 0.091 -0.067
от табл. 1).
Как видно из табл. 2, сходимость вычислительного процесса в первом варианте практически отсутствует. Контролируемое напряжение ам на краю оболочки остается весьма далеким от нуля, несмотря на достаточно густую сетку элементов дискретизации (216 КЭ) рассчитываемой конструкции. И, напротив, во втором варианте наблюдается быстрая сходимость вычислительного процесса при относительно редкой сетке конечных элементов (по сравнению с первым вариантом расчета). Меридиональное напряжение ам близко к нулевому значению.
В правой крайней колонке таблицы приведены результаты расчета компенсатора, полученные при использовании в качестве элемента дискретизации высокоточного четырехугольного конечного элемента с размером матрицы жесткости 72 х 72 [7]. Алгоритм формирования матрицы жесткости четырехугольного конечного элемента также включал процедуру векторной интерполяции перемещений.
Сравнительный анализ конечно-элементных решений второго варианта расчета и результатов, представленных в крайней правой колонке, показывает, что численные значения контролируемых параметров напряженно-деформированного состояния компенсатора практически совпадают, что также подтверждает достоверность результатов численного анализа оболочки во втором варианте расчета.
Таким образом, можно сделать вывод о том, что векторная интерполяция перемещений в треугольном конечном элементе позволяет кардинально повысить эффективность использования данного типа элемента при расчете оболочек сложной геометрии со значительными градиентами кривизны срединной поверхности. Применение треугольных дискретных элементов со стандартной интерполяционной процедурой к решению данного типа задач не приводит к получению приемлемых по точности результатов.
Пример 2. В качестве примера была рассчитана цилиндрическая оболочка, загруженная вдоль образующей равномерно распределенной нагрузкой интенсивности д. На диаметрально противоположной образующей цилиндра относительно приложения нагрузки д оболочка имеет пружинные опоры с переменной жесткостью, позволяющие смещаться ей вертикально вниз как абсолютно твердому телу (рис. 2).
Рис. 2.
Т аблица2
Точка Напряжение Вариант I Вариант II Высокоточный четырехугольный КЭ с размером матрицы жесткости 72 х 72 [7]
2 х 19 2 х 31 2 х 37 2 х 55 2 х 73 2 х 109 2х 19 2 х 31 2 х 37 2 х 55 2 х 55
1 &к 7.24 -102.4 -97.5 -49.86 -21.26 -0.24 10.77 11.47 11.57 11.60 11.41
& м 73.03 -345.6 -343.3 -191.8 -103.0 -36.55 -2.19 -0.81 -0.45 -0.12 0.049
2 &к 7.14 -102.9 -98.1 -50.49 -21.78 -0.58 10.75 11.45 11.55 11.57 11.41
&м 72.81 -346.0 -343.6 -191.7 -101.6 -36.09 -2.20 -0.82 -0.46 -0.135 0.049
ТаблицаЗ
Смещения оболочки как жесткого тела, м Вариант I Вариант II
&вн, Н/м2 ст|н, Н/м2 &нар, Н/м2 ст|н, Н/м2
0 194.8 ■ 104 194.4 ■ 104 -195.2 ■ 104 -195.2 ■ 104
0.01 233.8 ■ 104 154.2 ■ 104 -195.2 ■ 104 -195.2 ■ 104
0.05 385.4 ■ 104 -1.83 ■ 104 -195.2 ■ 104 -195.2 ■ 104
0.10 565.5 ■ 104 -187.2 ■ 104 -195.2 ■ 104 -195.2 ■ 104
54 Ю. В. Клочков, А. П. Николаев, Н.А. Гуреева
Были приняты следующие исходные данные: радиус цилиндра г = 0.1 м; толщина оболочки Ь = 0.001 м; модуль упругости материала Е = 2 ■ 1011 Н/м2; коэффициент Пуассона
V = 0; интенсивность равномерно распределенной линейной нагрузки д = 10 Н/м.
Расчеты были выполнены также в двух вариантах: в первом при формировании матрицы жесткости треугольного конечного элемента 27 х 27 применялась традиционная независимая процедура (2), (3); во втором варианте расчета была реализована аппроксимация полей векторов перемещений (4)—(15). Результаты расчетов оболочки представлены в табл. 3, где приведены численные значения кольцевых напряжений при сетке узлов 2 х 19 в зависимости от величины смещения цилиндра как абсолютно жесткого тела. Как видно, результаты расчетов по вариантам существенно отличаются друг от друга. Если во втором варианте кольцевые напряжения в наружных волокнах оболочки в точках 1 и 2 остаются неизменными, то в первом варианте погрешность вычисления кольцевых напряжений во внутренних волокнах оболочки в наблюдаемых точках с увеличением смещения цилиндра как абсолютно твердого тела стремительно нарастает до совершенно неприемлемых значений. Даже при отсутствии жесткого смещения (у = 0) напряжения в точках 1 и 2 различаются на 0.4-104 Н/м2, хотя при выбранной расчетной схеме они должны полностью совпадать (что наблюдается во втором варианте расчета).
Таким образом, можно сделать вывод, что предложенная интерполяция полей векторов перемещений треугольного дискретного элемента с размером матрицы жесткости 27 х 27 позволяет в неявном виде корректным образом в полной мере учесть смещения конечного элемента как жесткого тела.
Список литературы
[1] ГоловАнов А.И., Корнишин М.С. Введение в метод конечного элемента статики тонких оболочек. Казань, 1990. 269 с.
[2] Хечумов Р.А., Кепплер Х., Прокофьев В.Н. Применение метода конечных элементов к расчету конструкций. М.: Изд-во АСВ, 1994. 351 с.
[3] Зенкевич О. Метод конечных элементов в технике. М.: Мир, 1975. 542 с.
[4] Секулович М. Метод конечных элементов. М.: Стройиздат, 1993. 645 с.
[5] Новожилов В.В. Теория тонких оболочек. Л.: Судпромгиз, 1962. 432 с.
[6] Черных К.Ф. Линейная теория оболочек. Л.: Изд-во ЛГУ, 1962. Т. 1. 374 с. (1964. Т. 2. 395 с.).
[7] Клочков Ю.В., Николаев А.П., Киселев А.П. Применение четырехугольного элемента с матрицей жесткости 72 х 72 для расчета оболочечных конструкций // Изв. вузов. Сер. Строительство. 1998. № 4-5. С. 36-41.
[8] Клочков Ю.В., Николаев А.П., Киселев А.П. О функциях формы в алгоритмах фор-
мирования матриц жесткости треугольных конечных элементов // Изв. вузов. Сер. Строи-
тельство. 1999. № 10. С. 23-27.
Поступила в редакцию 25 февраля 2004 г., в переработанном виде —16 августа 2004 г.
