ТРЕУГОЛЬНЫЙ КОНЕЧНЫЙ ЭЛЕМЕНТ С ШЕСТЬЮ СТЕПЕНЯМИ СВОБОДЫ В УЗЛЕ Текст научной статьи по специальности «Физика»

Транспортное машиностроение. 2022. № 05(5). С. 22-29. ISSN 2782-5957 (print) Transport Engineering. 2022. no. 05(5). P. 22-29. ISSN 2782-5957 (print)

Научная статья

Статья в открытом доступе

УДК 519.63

doi: 10.30987/2782-5957-2022-5-22-29

ТРЕУГОЛЬНЫЙ КОНЕЧНЫЙ ЭЛЕМЕНТ С ШЕСТЬЮ СТЕПЕНЯМИ СВОБОДЫ В УЗЛЕ

Александр Викторович Яковлев10, Олег Дмитриевич Казаков2 Алексей Петрович Болдырев3

1 2 Брянский государственный инженерно-технологический университет, Россия, г. Брянск 3 Брянский государственный технический университет, Россия, г. Брянск

1 jav.05@mail.ru

2 kod8383@mail.ru

3 apb.tubryansk@gmail.com

Аннотация

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

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

нечного элемента с 12 степенями свободы в узле. Для тестирования полученной матрицы жесткости, была написана программа, основанная на MKЭ, с помощью которой выполнен расчет консольной балки. Погрешность расчета перемещений составила всего 0,22 %.

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

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

Ссылка для цитирования:

Яковлев А.В. Треугольный конечный элемент с шестью степенями свободы в узле /А. В. Яковлев, О. Д. Казаков, А. П. Болдырев // Транспортное машиностроение. - 2022. - № 5. - С. 22-29. doi: 10.30987/2782-5957-2022-522-29.

Original article Open Access Article

TRIANGULAR FINITE ELEMENT WITH SIX-DEGREE-OF-FREEDOM NODE

Aleksander Viktorovich Yakovlev1H, Oleg Dmitrievich Ka-zakov2 Aleksey Petrovich Boldyrev3

1 2 Bryansk State University of Engineering and Technology, Russia, Bryansk 3 Bryansk State Technical University, Russia, Bryansk

1 jav.05@mail.ru

2 kod8383@mail.ru

3 apb.tubryansk@gmail.com

The work is devoted to obtaining the stiffness matrix of a high-precision, flat finite element with 6 degrees of freedom in the node, for solving plane elasticity problems by the finite element method.

The scientific literature describes higher order elements. However, the theoretical results of these studies are quite far from their practical application. This paper gives a very detailed derivation of the stiff-

ness matrix of a high-precision finite element. Similarly, a stiffness matrix of a tetrahedral finite element with 12-degree-of-freedom node can be obtained. To test the obtained stiffness matrix, a program based on the FEM was written, with the help of which the cantilever beam is calculated. The error in calculating displacements is only 0.22%.

22

© Яковлев А.В., Болдырев А.П., Казаков О.Д., 2022

Conclusion: the stiffness matrix presented in the Keywords: problem, theory, elasticity, dis-

paper can be used with great success in numerical placement, deformation, stress, calculation, cantilever

methods for calculating the stress-strain state. beam.

Reference for citing:

Yakovlev A.V., Kazakov O.D., Boldyrev AP. Triangular finite element with six-degree-freedom node. Transport Engineering. 2022;5:22 - 29. doi: 10.30987/2782-5957-2022-5-22-29.

Введение

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

На рис. 1 показан высокоточный конечный элемент с 6 степенями свободы в узле, в котором, в качестве степеней сво-

боды используются перемещения и производные от перемещений. Такой конечный элемент будет обладать повышенной точность, поскольку связывает условиями непрерывности не только поля перемещений, но и поля деформаций.

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

Рис. 1. Треугольный конечный элемент, с шестью степенями свободы в узле Fig. 1. A triangular finite element with six degrees of freedom at the node

Поля перемещений

Обозначим узлы треугольного конечного элемента (рис 2) буквами /, } и к. Каждому узлу треугольника соответствуют координаты У^), (XУ])

(Хк, Ук).

и

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

u(L) = a1L3 + + + + + k^i + L[ + &qLj L fc + j Lfc + ai0LiLjLk

v(L) = ßiLf + ß2L] + ß3L3k + ß4L2iLj + ß5L]Lk + ß6L2kLt + ß7L)Lt + ß8LjL2k + ß9L2tLk + ßioLtLjLk

(1)

или более коротко:

{u(L)} = Ко] ■ {Lw}.

(2)

Рис. 2. Система координат Si,Sj,Sk Fig. 2. Coordinate system Si,Sj,Sk

Для того, чтобы в качестве степеней свободы использовать производные от перемещений вдоль сторон Бг(г = ¿, }, к)

д{и(Ь)}

треугольника, выполним дифференцирование зависимостей (1).

дБг дБг

Видим, что дифференцирование векторной функции {и(Ь)} свелось к дифференцированию векторной функции {^о}. Найдем соответствующие производные. Для чего, каждой стороне треугольника 1г(г = ¿, }, к) поставим в соответствие

= -^[aio№io] = [aio]diLlo}

(3)

х= fSi+Xj ц +

У

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

д{Ью} _ д{Ью} дЦ д{Ью} дБг дЬ^ дБг дLj

дБг '

координату Бг с началом в младшем узле (рис. 2).

Рассмотрим сторону /Л: треугольника. Связь между координатой и декартовой системой х, у задается очевидными соотношениями:

(4)

нате Бг (г = ¿, }, к) запишем следующим образом:

dLj | djLio) dLk dSr dLk dSr

(5)

Функции Ьр(р = ], к) входящие в (5), являются сложными функциями от координат х, у. Поэтому:

дЬр _ дЬр дх дЬр ду

dSr дх dSr

ду дБг

Имея в виду известные зависимости (4):

(6)

т.

дБг 2-Л-1г

Для треугольника известны соотношения

фк - ск^1 = С;ЪI - с1Ъ] = скЪ{ - с]Ьк = щ + а, + ак = 2Л, учитывая которые, для различных сочетаний р и г из набора /,}, к:

I-1 при г = Ь, р = к; г = }, р = г, г = к, р = ] _1 при г = Ь, р = }) г = }, р = к; г = к, р = /. при г = р

Подставляя эти зависимости в (9), получаем производные по конкретным переменным , , и Бк:

dSr { 1г 0

>

дЦю}

dS дЦю} dS

как

дЦю)

Эти зависимости можно представить произведение некоторой матрицы

i k\ i ij( =4

d{L

dLj d{Lio} dLk d{Lio} dLi

10} d{L10

}

+

+

d Lk d{Lio dL

d{Lio} dL

AJ

[Тг](г = I, ), к), строение которой очевидно, на вектор {Ь6}:

d{Lio} d S

1

= Т [ Tr]{ L6}, где rEi' J' к.

(7)

Где:

{ЬбУ = {Ь2, I'/■, Ь^Ь], ЪкЬ^ Подставив эти зависимости в (3) получаем:

д[и(Ь)}

dSr

= Т[ Я1о][ТгШб}, где г е i, j , к.

(8)

Векторы узловых перемещений

Введем вектор узловых перемещений [5]. С этой целью, выполним некоторые

д[и(Ь)}

преобразования (8). Умножим (8) на длину стороны треугольника Iг. Получим:

1Г = [ а10][Тг][Ьб]. (9)

Введем определение производной от функции [и(Ь)} вдоль некоторой стороны треугольника I г.

{и(I)} г = • и

dSr

dSr

В этом случае, зависимость (9) примет вид:

м/)Ь = [аю][ТгШб1

(10)

Теперь можно определить векторы узловых перемещений. Введем два локальных вектора узловых перемещений с компонентами, упорядоченными вдоль направлений х и у.

Uj

Uj„k Vj„k

UjJ

Vj„i

Uk Vk

ukj

Vk„i

Uk,j} Vkj}

(11)

Новое понятие производной обеспечивает независимость компонент матрицы [а10] от размеров стороны треугольника и делает зависимость (10) универсальной для любого конечного элемента при дифференцировании вдоль любой из его сторон.

{ и}Т = {и ии и1як

Индекс 5 здесь указывает на то, что производные взяты вдоль соответствующих сторон треугольника по формуле (10). Вместе эти два вектора образуют глобаль-

{ и}Т = {{ и}Т { У}Ц

Введем также три локальных вектора узловых перемещений для узла р треугольного конечного элемента

т

" иу,х Ур УриГ УРи1}. (11)

нечного элемента, компоненты которого

ный вектор перемещений для треугольного конечного элемента, с компонентами, упорядоченными по направлениям х и у.

U

упорядочены по узлам

{ ир)т = { ир"г

где р Е Ь, }, к; г Е Ь, }, к; t Е Ь, }, к. Вместе эти три вектора образуют глобальный вектор узловых перемещений для ко-

{ ир}Т = {{и}т { идт { ик}Т}.

Между векторами {и}5 и {ир} установим связь с помощью матрицы [Е-^]:

3 { и}5 = [Е-]{ 0р}5. 25

Строение матрицы [Ец] очевидно.

Вектор узловых перемещений с производными вдоль сторон треугольного конечного элемента мы будем использовать для определения компонент матрицы [а1с]. Использовать же этот вектор для вывода матрицы жесткости и расчете

{и(Ь)}л - д{и{1)]

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

{и(Ь)},2 =

д{и(Ь)}

(15)

д^ ,2 ду

В соответствии с (15) определим три локальных вектора узловых перемещений для узла

р треугольного конечного элемента:

{ иР}1 - { ир

и

и

Р,1 ир,2 у Ур,1 Ур,2}. (16)

где: р Е ¿, }, к. Индекс х в (16) указывает на то, что производные взяты вдоль координатных осей х и у. Из (9):

г гт -1 [д{и^)}дх , д{и(Ь)}ду\

{и(Ь)},,г= дБг 1г = 1г( дх дБг+ ду д5г)

У„

_ /д{и(Ь)} _ д{и(Ь)} \ Н дх Сг ду °ГУ

(17)

ду г>

Теперь можно установить связь между векторами { ир} из (11) и { "р}х из (16) с помощью матрицы [Мр] :

" Ы = [Мр\х{ир)х

где:

\Wihx [Об] [Об] [06] Щ5х [Об] [Об] [Об] [М^

Здесь [Об] нулевая матрица размером 6*6.

[МрЪХ =

1 О О 1 О О 1 О О

Шзх = О С] -Ъ] - х О к -Ък - х О Сг -Ъг

О Ск -Ък О —ЪiJ О ] - Ъ]

Определим коэффициенты а1, (где I Е 1, 1, ..., 1О) интерполяционных полиномов (1) Первые 9 коэффициентов определим из условий непрерывности перемещений и производных вдоль сторон треугольника в его узлах.

Для определения коэффициентов а10 и р10 необходим ещё один узел. Этот узел можно задать в центре тяжести треуголь-

ного элемента ствует много

{и(Ь)} - [П(Ь)] \{^с}Т

ДОю}7

Где {и(ь)} = ота, Ш)] =

Введем матрицу функций формы:

[М(Ь)] - [П(Ь)]

Тогда окончательно:

что нежелательно. Суще-способов доопределения «лишних» неизвестных, которые довольно подробно излагаются в литературе. В результате, можно получить матрицу преобразования [2], связывающую компоненты вектора перемещения внутри конечного элемента с узловыми перемещениями следующей зависимостью:

[г] [О]1 [О] [г]] {О1с}Т1 {^сЛ

№5. { й}5-

{ и} {У}5

[г]

[О]

[О]

[г]

[ит = [адя ор}. (18)

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

ные от перемещений вдоль координатных осей х и у, для этого, подставим зависимость (17) в (18):

'X

понентами которого являются производ-

[иШ} = №ажм]зх{Ор}х. (19)

Деформации и напряжения

Введем вектор деформаций:

[г}7 = [ех тху еу} = [^11 £12 £22}. и запишем уравнения Коши в виде:

[е(Ь)} = Шит, (20)

Подставим (19) в (20):

[а(ь)} = тмажм]5х{Ор}х.

В результате действия матрицы [5] на матрицу [N(1)] будет порождена матрица градиентов [В(Ь)]:

[г(Ь)} = [В(1)][М]3х{Ор}х. (21)

Уравнение (21) дает искомую зависимость между деформациями и узловыми перемещениями.

Введем вектор напряжений:

[а}Т = [^х *ху ау} = [а11 а 12 ^22}. и запишем связь между напряжениями и деформациями в виде:

[а(Ь)} = [Э][Е(Ь)}, (22)

чим, если в (22) подставить уравнения (21):

где [Э] - матрица упругости размером 3 X 3 компоненты, которой определяются законом Гука Искомую зависимость полу-

[а(Ь)} = [0][В(1)][М]5х{0р} . (23)

Система разрешающих уравнений МКЭ

Введем обозначения: А - работа внешних сил

П - потенциальная энергия деформированного тела Л - энергия системы внешних и внутренних сил

Для вывода разрешающих уравнений МКЭ воспользуемся формулой Лагранжа [6]:

8Л = 0.

Л = П -А. (24)

Потенциальная энергия определяется формулой Клайперона [6]:

П=±Цг}Т[а}с1У. (25)

Найдем потенциальную энергию Пу для конечного элемента с номером V, подставив в (25) уравнения (21) и (23):

Пу = 1 !у{Ор}Т [Му]Тх [Ву(Ь)]Т[Бу] [Ву (Ь)] [Му]3х{Ор}хйУ. (26)

Работа внешних сил для узла р конечного элемента с номером V:

Т

ААр = {Ор} {рУ}, где: { Оу}Т = [Ор Ору}, {рру}Т = [Ррх РрУу}.

С помощью матрицы [Р2], введем в зависимость (26) вектор { Ор} из (16)

Т

Ар = Ш^рру},

где:

и = [

1 0 0 0 0 0

ООО1ОО

Тогда, для конечного элемента с номером :

Т

Ар - Т^-Лир}^]^ Выполнив суммирование по всем коне ч2н7ым элементам, получим:

Л -1тЕУ-11у{0р}Тт[мр]Тх[вр(ь)]Т[ор][вр(ь)][мр]5х{0р}йУ-

JVl up}x[Mv

т

n=iY?v=i{uv} [Р2][РЛ

Минимизируя полученный функционал приходим к разрешающей системе уравнений

МКЭ:

lUL[Mv]L[Bv(L)]T[Dv][Bv(L)][Mv]sx{up}dv = ZUZU^PVl

(27)

Матрица жесткости конечного элемента

Интеграл в левой части (27) представляет собой матрицу жесткости конечного элемента [ к]. Считая, что толщина

(Ь)=ЬгЬг +И]Ь]+Ьк1к.

где: Иг, К], Ик - толщина конечного элемента в узлах / , }, к и учитывая, что компоненты матрицы [Му]зх являются посто-

конечного элемента изменяется линейно, запишем:

янными величинами, получаем матрицу жесткости конечного элемента:

[kv] = [Mv]TsxUBv(L)]T[Dv][Bv(L)]dV[Mv]sx.

(28)

Тестовый расчет

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

консольной балки (Рис. 3). Расчет имеет чисто математический характер, поэтому величины и размерности физических величин не приводятся.

\ \ \ \ \ \ --< Р

7

Рис. 3. Триангуляция консольной балки для выполнения тестового расчета

Fig. 3. Triangulation of a cantilever beam for performing a test calculation

Расчет перемещений для точки А балки (Рис. 3), по точным формулам, составил величину 2.744, а при расчете по

методу конечных элементов получено перемещение 2.738. Таким образом, погрешность расчета составила величину 0.22%.

Заключение

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

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

СПИСОК ИСТОЧНИКОВ

1. Корнеев В.Г. Схемы метода конечных элементов высоких порядков точности. - Л.: Изд-во Ленингр. ун-та, 1977. 208 с.

2. Галлагер Р. Метод конечных элементов. Основы. - М.: Мир, 1984. - 428 с.

3. Сегерлинд Л. Применение метода конечных элементов. - М.: Мир, 1979. - 392 с.

4. Zienkiewicz O.C., Taylor R.L. The Finite Element Method. Vol. 1: The Basis. - Butterworth Heinemann, 2000. - 707p.

5. Баландин, М.Ю., Шурина Э.П. Векторный метод конечных элементов: Учеб. пособие. — Новосибирск: Изд-во НГТУ, 2001. — 69 с.

6. Ern A., Guermond J.L. Theory and practice of finite elements. Applied Mathematical Sciences. 2004;159.

1. Korneev VG. Schemes of the finite element method of high order accuracy. Leningrad: Publishing House of the Leningrad University; 1977.

2. Gallagher R. Finite element method. Fundamentals. Moscow: Mir; 1984.

3. Segerlind L. Application of the finite element method. Moscow: Mir; 1979.

4. Zienkiewicz OC, Taylor RL. The finite element method: The basis. Butterworth Heinemann; 2000.

5. Balandin, MYu, Shurina EP. Vector finite element method. Novosibirsk: Publishing house of NSTU; 2001.

Информация об авторах:

Яковлев Александр Викторович, кандидат технических наук, доцент, преподаватель Брянского государственного инженерно -технологического университета, тел. +7(919)197-62-11. Казаков Олег Дмитриевич, кандидат экономических наук, доцент, проректор по цифровизации

Yakovlev Aleksander Viktorovich, Candidate of Technical Sciences, Associate Professor, Lecturer of Bryansk State University of Engineering and Technology, phone: +7(919)197-62-11. Kazakov Oleg Dmitrievich, Candidate of Economic Sciences, Associate Professor, Vice-Rector for Digital-

7. Zienkiewicz O.C., Taylor R.L. The Finite Element Method. Vol. 2: Solid Mechanics. - Butterworth Heinemann, 2000. - 459 p.

8. Shames I.H., Cozzareli F.A. Elastic and Inelastic Stress Analysis, revised edition. - Washington: Taylor & Francis, DC, 1997. - 187 p.

9. Di P.D.A., Ern A. Mathematical aspects of discontinuous Galerkin methods. Mathématiques et Applications. 2012; 69.

10. Стринг Г., Фикс Дж. Теория метода конечных элементов. - М.: Мир, 1977. - 350 с.

6. Ern A, Guermond JL. Theory and practice of finite elements. Applied Mathematical Sciences. 2004; 159.

7. Zienkiewicz OC, Taylor RL. The finite element method: Solid mechanics. Butterworth Heinemann; 2000.

8. Shames IH, Cozzareli FA. Elastic and inelastic stress analysis, revised edition. Washington: Taylor & Francis, DC; 1997.

9. Di PDA, Ern A. Mathematical aspects of discontinuous Galerkin methods. Mathématiques et Applications. 2012; 69.

10. String G, Fix J. Theory of the finite element method. Moscow: Mir; 1977.

Брянского государственного инженерно -технологического университета, тел. +7(920)604-65-36.

Болдырев Алексей Петрович, доктор технических наук, профессор Брянского государственного технического университета, тел. +7(910)331-20-00.

ization at Bryansk State University of Engineering and Technology, phone: +7(920)604-65-36. Boldyrev Aleksey Petrovich, Doctor of Technical Sciences, Professor of Bryansk State Technical University.

Вклад авторов: все авторы сделали эквивалентный вклад в подготовку публикации. Contribution of the authors: the authors contributed equally to this article.

Авторы заявляют об отсутствии конфликта интересов. The authors declare no conflicts of interests.

Статья опубликована в режиме Open Access. Article published in Open Access mode.

Статья поступила в редакцию 22.12.2021; одобрена после рецензирования 21.03.2022; принята к публикации 21.04.2022. Рецензент - Алгабачиев А.Ю., доктор технических наук, зав. отделом Института машиноведения им. А.А. Благонравова РАН, член редколлегии журнала «Транспортное машиностроение».

The article was submitted to the editorial office on 22.12.2021; approved after review on 21.03.2022; accepted for publication on 21.04.2022. The reviewer is Algabachiev A.Yu., Doctor of Technical Sciences, Head of the Department at Mechanical Engineering Research Institute of the Russian Academy of Sciences, member of the Editorial Board of the journal Transport Engineering.