CC BY
68
Проблемы теории упругости
МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ В РЕШЕНИИ ТРЕХМЕРНЫХ ЗАДАЧ ТЕОРИИ УПРУГОСТИ
А. П. КИСЕЛЁВ, канд. техн. наук, доцент
Волгоградская государственная сельскохозяйственная академия
Для решения трехмерных задач по расчету на прочность инженерных конструкций требуются высокоточные конечные элементы. В зависимости от формы конструкции, представляемой дискретной моделью это могут быть тетраэдры, треугольные призмы или восьми узловые шестигранники. В настоящей работе на основе соотношений механики сплошной среды разработаны алгоритмы необходимые для формирования матриц жесткости высокоточных конечных элементов в виде треугольной призмы и объемного восьми узлового элемента. За узловые неизвестные конечных элементов выбирались компоненты вектора перемещения и их первые производные [1,2].
Предложен вариант формирования функций формы на основе полных алгебраических полиномов в комбинации с полиномами Эрмита для конечного элемента в виде призмы с треугольным основанием. Положение произвольной точки сплошной трехмерной среды в недеформированном состоянии в декартовой системе координат может быть определено радиус-вектором (рис. 1)
К°=хк(вт%, {к,т =1,2,3), (1)
где хк, 4 - переменные и орты декартовой системы координат, соответственно; 0" - переменные криволинейной системы координат, отличной от декартовой.
Дифференцированием (1) по криволинейным координатам вт определяются ковариантные векторы локального базиса в исходном состоянии
(2)
1 дв'
Ковариантные компоненты метрического тензора определяются скалярными произведениями базисных векторов
О -О -О /оч
ац ~ ' а] ■ О)
Производные векторов локального базиса (2) определяются из выражения
п да0: „ „.
дв*
Символы Кристоффеля второго рода Гд определяются через символы
Кристоффеля первого рода
Г%=а°«Г°=Г1 (5)
¡к - " * ¡¿к ~1 к/
+
V
дхк дх] дх1
Г°к: - символы Кристоффеля первого рода.
а;
Положение точки сплошной среды после деформации будет определяться
радиус-вектором & = + У. (6)
Ковариантные векторы локального базиса в деформированном состоянии определяются путем дифференцирования (6)
(7)
Компоненты метрического тензора в деформированном состоянии определяются скалярным произведением векторов
а„ =3,3, = (8)
Вектор перемещений произвольной точки сплошной среды может быть
представлен в виде
Г = а?и + а°2у + а°3м>, (9) где и, V, ц> - проекции вектора перемещений.
Производные вектора перемещений V в криволинейной системе координат:
V
ди"
.к
да
Рис. 1. Перемещение произвольной точки сплошной среды в результате деформации
двК "' дВк Ла1+/к2а2+/к%,(к= 1,2,3), (10)
ггтр {} /'.2 f? — мнпгпчпенм сп-
J к и к и к ---------------» " -
держащие компоненты вектора перемещений и их первые производные. Ковариантные компоненты тензора деформаций определяются на основе соотношений механики сплошной среды [3]
-в;)/2. (П)
Контравариантные компоненты
тензора напряжении для упругого тела определяются через ковариантные компоненты тензора деформаций выражением
сгу =А'*ИетИ, (и ,т,п = 1,2,3)
= Л аща™" + м (ашаор + аша*]т);
/ути _ 5 „01/ Отп
(12) (13)
где
- контравариантные компоненты метрического тензора, полученные в исходном состоянии. На основе (13) можно записать матричную зависимость между вектором контравариантных компонент тензора напряжений и вектором ковариантных компонент тензора деформаций
И=№}, (14)
где {<т}г = {о-и,а22,сг33,сги,а23,а31}- вектор напряжений;
15 £22' % >2% >2^23 >2^31,} - вектор деформаций; [с] - матрица упругости размером 6x6.
При представлении трехмерного континуума дискретной моделью часто используется элементы в виде треугольной призмы (рис. 2). Для выполнения численного интегрирования по объему конечного элемента произвольная треугольная призма в глобальной системе координат х!, х2, х3 отображается на треугольную призму в локальной системе координат £ т], С, (рис. 3) с координатами узлов /(0;0;-Цу(1;0;-1), Л(0;1;-1), /(0;0;1), ю(1;0;1), и(0;1;1).
Зависимость между координатами х', х2, х3 и локальными £ т], £ конечного элемента определяется соотношением
где
координаты узловых точек в системе х', х2, х3; а = 1,2,3.
х-.
Рис. 2. Конечный элемент в виде треугольной призмы в системе координат
Из (15) определяются производные х1, х2, х3 координат в локальной системе координат и производные локальных координат в координатной системе хх ,х2,хг.
За узловые неизвестные конечного элемента выбирались компоненты вектора перемещений и, V и ы в направлениях осей координат х1, х2 и х3
их первые производные и/, и 2, и3, V/, Чл у.1, Ы з в угловых точках. Таким образом, вектор перемещений узловых точек конечного элемента в глобальной системе координат будет иметь вид
{г/ -
1X72 (16)
и'2, и'2, ил, и '2, и™, и"2, и з, Ы}ъ, Из, М3, И 3 , И " , V*,.. •, V,...}.
Для аппроксимации поля перемещений внутренней точки конечного элемента через его узловые значения в треугольной призме используются два типа функций формы [4].
Для аппроксимации полей перемещений внутренних точек в координатах £ т1 конечного элемента использовался полный алгебраический полином третьей степени, который в этом случае будет содержать десять членов
£(£/7) = *:, +к2% + кгг] + к^2 +к^г] + к6т]2 +к^2т] + кд%т]2 +к]0т]\ (17)
где к], ...,кю - неизвестные коэффициенты.
Для получения дополнительного условия, необходимого для определения всех коэффициентов аппроксимирующих функций в вектор узловых неизвестных, кроме узловых перемещений и их первых производных добавлялась смешанная производная перемещения узла
q'¿r]. Вектор узловых неизвестных в
этом случае будет иметь следующую структуру [4]
' * «-/ л.' «У „к „1 „] к
А (18)
1x10
В результате количество неизвестных оказывалось равным числу уравнений полученной системы.
Перемещение внутренней точки элемента можно выразить через узловые неизвестные матричным соотношением
1*10 1*10
Смешанную производную перемещения первого узла д'^ с использовани-
(ОН'.О)
Рис. 3. Конечный элемент в виде треугольной призмы в системе координат г}, С,
ем метода конечных разностей в локальной системе координат можно выразить через первые производные узловых перемещений
С использованием (20) вектор-столбец узловых перемещений с десятью компонентами {<?* можно определить через вектор-столбец
1ч10
КГ = У.^Л.^^'^'^п'Уп}
1x9
имеющий только девять составляющих компонент
I<21>
1x10 10x9 9x1
где [Ь] - матрица преобразования.
Принимая во внимание (21) выражение (19) примет вид
Ч = {<р)Т{ь\ (22)
1*9 9x1
гле
. Пеоемещение внутоенней точки (22) в развернутом виде
Ш0 10x9
будет определяться выражением
д = Ох пУ+С2 {§, т/У + (§, г]У + + <?4 чЦ + £, пУ4 + с6 & пЦ + (23)
+о, {§, пУ, + о, пУ + <?9 (с, пУп,
где 01, ..., Од - аппроксимирующие функции для треугольной области конечного элемента. С учетом (23) перемещение точек, расположенных внутри призматического конечного элемента, будет определяться через перемещения узловых точек в локальной системе координат
Ч = {§, пУн Ш + С2 (СУ + 03 {4, //К {сУ + +с, (£, 7)л2 {сУ + 02 (<?, п)ь2 (сУ + с3 тт)и2 (сУ +
+о4 (<?, {су + о5 ^ (су + о6 (#, {су + (25)
+ с4 (4,[су +05& 1т)и2 {су + с6 [§, п)ьг (су +
+а7 (£ п)ъ2 {сУ„ + с8 (<г, ^ + с9 п)к2 (сУ + + О, £, + & + С3 & ^ Шп +
+о, (4, + + с3 {4, - ,
где /I, И2 (С), К (с), кА (С) - полиномы Эрмита третьей степени:
кх{с)Мс'-X+ н2{с)=-Нсъ-ъс-2}
(26)
Мс)=~{с3-с2-с+\)к{сЬ~{с3 +с2-с А
В качестве объемного восьми узлового конечного элемента выбирался произвольный шестигранный элемент с узлами
1(1),20),3(к),4(1),5(т),6(п),7(р),8(И) (рис. 4).
Положение узловых точек конечного элемента определялось тремя координатами - х', х2, х3. Для выполнения численного интегрирования по объему ко-14
немного элемента произвольный восьми узловой элемент отображался на куб в системе координат £77,«¿'(рис. 5). Локальные координаты кубического элемента изменялись в пределах -1 < < 1.
Зависимость между глобальными координатами х',х2,х3 и локальными £ г], С, конечного элемента определяется выражением
X = X,
8
■ + х"
8
8 ' 8
— + -
+
(27)
8
8
где а ~ \,2,3;х",х",х", х", х„, х", хар, х" - координаты узловых точек конечного элемента в глобальной системе координат.
Дифференцированием (27) определялись производные х'\х2,х' - координат в локальной системе координат и производные локальных координат в глобальной системе. Вектор перемещений узловых точек конечного элемента в глобальной системе координат будет иметь вид:
ш Ь
/
Рис. 4. Произвольный восьмиугольный конечный элемент в системе координат х*. X3
Рис. 5. Восьмиугольный конечный элемент в системе координат
4.Г1.С
(28)
96X1
Перемещения внутренних точек конечного элемента определяются через перемещения узловых точек с помощью аппроксимирующих функций
Ч = У)Т{яу}, (29)
1*32 32x1
где {^¿/}7 - вектор-строка тридцати двух аппроксимирующих функций, осно-
1x32
ванных на полиномах Эрмита третьей степени.
Матрицы конечных элементов формировались на основе условия равенства работ внутренних и внешних сил. Интегрирование по объему конечных элементов выполнялось численно.
Пример 1. В качестве примера расчета исследовалось напряженно- деформированное состояние объемного тела при действии на него сосредоточенной
силы (рис. 6). Были приняты следующие исходные данные: величина сосредоточенной силы Р = 39.2-104 Н; сторона куба а = 3.05 м. Использование симметрии позволяет записать краевые условия для перемещений в виде: и = V = м> = 0.0 на грани АВСО, и = 0.0 на АЕНБ и V = 0.0 на АЕББ, все другие границы свободные.
Расчет выполнялся в двух вариантах: с использованием конечного элемента в виде треугольной призмы и с использованием восьми узлового шестигранного конечного элемента. При дискретизации объема отдельными треугольными
призмами вначале пространство х5 разбивалось на отдельные объем-
ные прямоугольники, которые затем разбивались на две призмы.
н
N 1 V* 1 "S. |
Рис. 6. Задача о действии сосредоточенной силы на объемное тело
Рис. 7. Защемленная панель Таблица 1
№ Г1.П. Размер конечно-элементной сетки Величина вертикального перемещения куба под силой \ух10"6м
Треугольная призма Восьмиугольный элемент
1 1x1x1 3.3 3.5
2 2x2x2 6.8 6.5
3 3x3x3 10.5 10.2
4 4x4x4 14.2 13.9
5 5x5x5 17.9 17.5
6 5x5x6 18.0 17.8
7 5x5x10 18.4 18.1
8 6x6x6 19.2 18.9
Результаты расчетов приведены в таблице 1. В таблице даются численные значения величины вертикального перемещения куба под силой, полученные при различном количестве конечных элементов дискретной сетки. В 3-й колонке таблицы приводятся результаты, полученные с использованием элемента в виде треугольной призмы, в 4-й с использованием восьмиугольных конечных элементов. В обоих вариантах проверялась сходимость вычислительного процесса с различным числом точек интегрирования и при увеличении числа дискретных элементов.
Пример 2. Рассчитывалась защемленная с краю прямоугольная панель, загруженная линейной распределенной нагрузкой (рис.7). Были приняты следующие исходные данные: величина распределенной нагрузки # = 39.2-104 Н/м; длина панели Ь = 1,1 м и ширина Ъ = 1,0 м; толщина А = 0,395 м. Величину про-
гиба панели можно определить аналитически.
Таблица 2
№ Дискретная КЭ в виде Восьмиузловой
П.П. сетка треугольной призмы конечный элемент
NNxNMxNL
1 lxlxl 0.000345 0.000357
2 1x10x1 0.000359 0.000359
3 1x15x1 0.000362 0.000361
4 1x20x1 0.000363 0.000362
5 1x25x1 0.000363 0.000362
Аналитическое решение: 0.000359 м.
В таблице 2 приводятся численные величины максимального прогиба панели под силой, полученные с использованием разработанных элементов. Как видно, результаты расчетов, полученные различными конечными элементами, практически совпадают.
В обоих примерах расчета наблюдается сходимость вычислительного процесса с увеличением числа дискретных элементов, а также совпадение значений полученными элементами и результата вычисленного аналитически (пример 2), что говорит о достаточной точности разработанных алгоритмов.
С использованием предложенных элементов можно значительно повысить точность расчета инженерных конструкций с ярко выраженными зонами концентраторов напряжений. Разработанные конечные элементы также могут эффективно использоваться в программных комплексах научно- исследовательскими и проектно-конструкторскими организациями, деятельность которых связана с проектированием и эксплуатацией сложных инженерных конструкций для уточненного расчета на прочность конструктивных элементов нефтегазового и химического оборудования, строительных и других конструкций. Это позволит проектировать экономически наиболее выгодные конструкции с обеспечением их надежной эксплуатации.
Литература
1. Киселёв А. П., Николаев А. П. Объемный конечный элемент в виде треугольной призмы с первыми производными узловых перемещений// Известия вузов, сер. «Строительство». - 2006. - № 1.
2. Киселёв А. П., Николаев А. П., Юшкин В. Н. Восьмиугольный конечный элемент для расчета толстостенных оболочек вращения// Сб. трудов Межд. на-учно-техн. конференции «Актуальные проблемы механики оболочек». - Казань, 2000. - С. 27-30.
3. Седов Л. И. Механика сплошной среды. - Т. 1-2. - М.: Наука, 1976.
4. Николаев А. П., Клочков Ю. В., Киселёв А. П. О функциях формы в алгоритмах формирования матрицы жесткости в треугольном конечном элементе// Известия вузов, сер. «Строительство». - 1999. - № 10. - С. 23-27.
FINITE ELEMENT METHOD FOR SOLUTION OF THREE-DIMENSIONAL PROBLEMS OF THEORY OF ELASTICITY
A. P. Kiselyov
In the present work, on the basis of relations of mechanics of a continuous medium the necessary algorithms for forming stiffness matrixes of triangular finite elements and volumetric eight node elements are developed.
The variant of forming of functions of the form is offered on the basis of complete algebraic polynomials in a combination with Hermite polynomials for a finite element as a prism with the triangular base.
