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

И. В. Михайлов, А. В. Синельщиков

МЕТОДИКА ВЫВОДА МАТРИЦ ЖЕСТКОСТИ ЛИНЕЙНО-УПРУГИХ ОБЪЕМНЫ1Х КОНЕЧНЫ1Х ЭЛЕМЕНТОВ

В практической деятельности инженеров используется перечень программных продуктов, реализующих расчет напряженно-деформированного состояния несущих систем и структур различного технического назначения: в автомобильной промышленности, судостроении, краностроении (ЛИРА, AnSys, Nastran, CAN и др.). Однако их использование стесняет инициативу инженеров, когда требуется развитие алгоритмических построений.

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

- сейсмических воздействиях;

- транспортных нагрузках рабочего состояния для автомобильных кранов;

- при ветре нерабочего состояния и др.

Известны способы получения матриц жесткости четырехузловых тетраэдрических конечных элементов, предложенные Галлагером [1]. Однако разработке строгого математического аппарата и верификации решения по получению матриц жесткости тетраэдрических конечных элементов с 4-мя или 10-ю узлами посвящено ограниченное число работ.

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

Сущность алгоритма заключается в следующем.

Введем местную систему координат O^nZ, координатные оси которой совпадают с тремя ребрами тетраэдра (рис. 1). Назовем эту систему местной системой координат первого рода (МСК1). Введем также прямоугольную декартову систему координат O2xyz, имеющую следующие особенности:

O2Х ° , 02 ° 0^ (1)

(O2xy)°(о1Хл). ( )

Эту систему назовем местной системой координат второго рода (МСК2) (рис. 2).

Рис. 1. Косоугольная система координат МСК1

Рис. 2. Относительное расположение координатных систем МСК1 и МСК2: /. 7, к — единичные векторы МСК2; /1, у'1, к1 — единичные векторы МСК1

Аппроксимирующие функции для компонентов перемещения (в МСК1) по объему тетраэдра:

и(1) =аі +а2 -Х + аз 'Л + а4 'С у(1) =а5 +а6 'Х + а7 л + а8 'С w(1) = а9 + аш 'Х + аіі 'Л + аі2 'С

или в матричной форме

{у (і) }=иН4

(2)

(3)

где

[А] =

і X л С 0 0 0 0 0 0 0 0

0 0 0 0 і X л С 0 0 0 0

0 0 0 0 0 0 0 0 і X л С

(4)

(1)

(5)

{а}=(аі а2 аз а 4 аз аб а7 а8 а

9 аі0

а

аі2Г.

(6)

За основные неизвестные принимаются узловые значения перемещений:

и

V

4і

{„ }=(и<

(!) л/1)

w

(1)

1

,(1) „(1)

*2

2

w

(1) (1)

...У

(1)

w

(1) (1)

(1)

4

w

(1)

2

Г

(7)

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

{я(1) }=(я

■1$

к Я

Ящ Я1^ я

Я3С Я4$

2$

Я2ц Я2с1

я

3$...

Я

4^

(8)

Дальнейшая задача состоит в установлении связи между перемещениями (7) и реакциями (8):

{я'" }=к(1) IVі},

(9)

где к о ] - искомая матрица жесткости конечного элемента. Учитывая свойства (1) системы координат МСК1, имеем:

$1 = о Л1 = 0 Сі = 0

Л2 = 0 С 2 = 0

$3 = 0 С 3 = 0

$4 = 0 Л4 = 0 Тогда из (2) и (10) следует:

(10)

и а) = а"

у(1) = а

w = а9

= а1 +$2 а

= а5 + $2 а

+ а = а1

= аі +П3 а

= а5 +^3 а

= а9 +^3 а1

= а1 + С 4 а

= а5 + С4 а

4 = а9 + С4 а1

(11)

или в матричной форме

{?(1) }=М-Ы.

(12)

Матрицу В легко получить, продифференцировав каждое из выражений (11) по аг- (і = 1... 12):

1

3

4

4

B =Э?®' (і3)

в=Ю7. (13)

Тогда коэффициенты из (2) и (6) будут:

{а}= [в]_" -{д(" }. (14)

Подставив (14) в (3), получим аппроксимирующие функции:

U(1) }= [л] - [в]-1 -{?(1)}. (15)

Введем два вектора, определяющих положение произвольной межуз-ловой точки тетраэдра в МСК1 и МСК2 соответственно:

{/»}=($ я С)'; (16)

{г(2) }=(х У 2)' . (17)

Векторы (16) и (17) связаны между собой следующей зависимостью:

{г(2) }=['і_2 ]-{(1) }, ("8)

где [Т"_2 ] — матрица преобразования координат

Tl-2 ] =

1 sin(0) cos(P)

0 cos(0) cos(g)

0 0 д/і - cos2(b) - cos2(g)

(19)

cos (0) - направляющий косинус оси n в МСК2; cos (P) и cos (у) -направляющие косинусы оси Z в МСК2 (см. рис. 2).

Рассмотрим компоненты ы<2\ v(2), w(2) перемещения произвольной межузловой точки тетраэдра в декартовой прямоугольной системе координат МСК2:

{U(2) }=(u(2) v(2) w(2) J; (20)

{и(2) }= Т;_2 ]-{и(1) }= [7;_2 ] -[л] -[в]"1 -fe(1)}. (21)

Подставим в матрицу [А] формулы (21) выражения для координат межузловой точки в МСК1 (т. е. от n, Z) из формулы

{r (1) }=[Т1-2 ]-1 -{Г(2) } (22)

и, таким образом, выразим компоненты перемещения в МСК2 межузловой точки через координаты межузловой точки в МСК2 (т. е. через x, y, z).

Компоненты деформации удобнее искать в МСК2. Запишем уравнения Коши:

є = -

Є уу =

ди(2)

дх ду^

дУ

дw(2)

Є 22 дг

ди(2) ду(2)

У =---------------1--------

їху ду дх

ду(2) дw(2)

У уг = У ху =

уг дг ду

дw(2) ди(2)

- + -

дх

дг

где бхх, Еуу, егг, Уху, Ууг, угх - компоненты тензора деформаций среды в межузловой точке в МСК2. Запишем равенства (23) в матричной форме:

где

{е(2) }=(е

{є(2) }=[Адиф ]7 •{и(2)},

_д_ дх

0

0 0 — 0 — —

уу Є гг У ху У уг У гх

0 0 д ду 0 д " дг

д 0 д д 0

ду дх дг

0 д 0 д д

дг ду дх

(24)

(25)

(26)

После выполнения всех операторов дифференцирования формула (24) в общем случае дает зависимость вектора {е(2)} (25) от координат межузловой точки в МСК2 (т. е. от х, у, г). Далее необходимо найти зависимость вектора {е(2)} от координат межузловой точки в МСК1 (т. е.

от £,, п, С). Для этого в компоненты вектора {е(2)} нужно подставить

выражения для х, у и г из формулы (18).

Запишем обобщенный закон Гука для компонентов деформации в МСК2:

°хх = (2 • & + 1) ехх + 1 • Єуу + 1 • Єг.

°уу = 1 • Єхх + (2 • & + 1) • Єуу + X • Є г

°гг = 1 • Єхх + 1 • Єуу + (2 • & + I) • Єг;

*ху = &• Уху

* уг = &• У у2

* гх = &• У гх

Є

где О - модуль сдвига, I - параметр Лямэ, о**, от огг, т^, V, тгх -компоненты тензора напряжений среды в межузловой точке в МСК2. В матричной форме

где

{с<2) }=[£,]-|е<2) },

УУ

ху

<28)

<29)

№

+ о 1 1 0 0 0

1 2-О + 1 1 0 0 0

1 1 2-О + 1 0 0 0

0 0 0 О 0 0

0 0 0 0 О 0

0 0 0 0 0 О

Удельная потенциальная энергия деформации:

п0 = 2-{*"’ Г-{а121 }= |е<!)Г •[£,]•{="’}.

Потенциальная энергия деформации тетраэдра с учетом <31):

п= п 0= Ш п о**[Гь2 ] - ¿Х ¿л - ¿С.

V V

Из теоремы Клапейрона о работе внешних сил следует:

п= 1-{«">Г •{?,1)}= 2-{?<цГ -И] •{?<1)},

откуда

к(1) =

I л

д2 п

<30)

<31)

<32)

<33)

<34)

Очевидно, <34) является решением поставленной задачи.

Предложенный авторами алгоритм апробирован на численном примере для 4-узлового тетраэдра: в программном комплексе Лп8у8 был построен тетраэдр как единый <4-узловой) конечный элемент, закрепленный в трех вершинах по всем степеням свободы, а в четвертой

т

вершине — по всем, кроме перемещения вдоль одного из примыкающих ребер. В этом направлении вершине было задано единичное усилие (узловая сила). После выполнения расчета программой было определено узловое перемещение, действующее в этом направлении, и сравнено с соответствующим компонентом матрицы, обратной к матрице жесткости [К(1)]. Расхождение результатов составило 0,005 %, что объясняется, по всей видимости, использованием в программе Лп8у8 точно такого же конечного элемента.

Матрица жесткости рассмативаемого здесь четырехузлового тетраэдра, таким образом, имеет вид:

К« ]=

" К,1 К 1,2 К1,3 К1,4 К1,5 К1,6 К1,7 К1,8 К1,9 К1,10 К1,11 К1Д2 '

К 2J К 2,2 К2,3 К2,4 К2,5 К2,6 К2,7 К2,8 К2,9 К2,10 К2,11 К2,12

К3,1 К 3,2 К3,3 К3,4 К3,5 К3,6 К3,7 К3,8 К3,9 К3,10 К3,11 К3Д2

К 4Д К 4,2 К4,3 К4,4 К4,5 К4,6 К4,7 К4,8 К4,9 К4,10 К4,11 К4,12

К 5Д К 5,2 К5,3 К5,4 К5,5 К5,6 К5,7 К5,8 К5,9 К5,10 К5,11 К5,12

К 6,1 К6,2 К6,3 К6,4 К6,5 К6,6 К6,7 К6,8 К6,9 К6,10 К6,11 К6,12

К 7Д К7,2 К7,3 К7,4 К7,5 К7,6 К7,7 К7,8 К7,9 К7,10 К7,11 К 7,12

K8J К8,2 К8,3 К8,4 К 8,5 К8,6 К8,7 К8,8 К8,9 К8,10 К8,11 К8,12

К 9,1 К 9,2 К9,3 К9,4 К9,5 К9,6 К9,7 К9,8 К9,9 К9,10 К9,11 К9,12

К10Д К10,2 К10,3 К10,4 К10,5 К10,6 К10,7 К10,8 К10,9 К10,10 К10,11 К10,12

Кп,1 К11,2 К11,3 К11,4 К11,5 К11,6 К11,7 К11,8 К11,9 К11,10 К11,11 К11,12

К12,1 К12,2 К12,3 К12,4 К12,5 К12,6 К12,7 К12,8 К12,9 К12,10 К12,11 К12,12

(35)

Ввиду ограниченности объема настоящей статьи приведем только один элемент матрицы жесткости (35):

К (‘) = и

008'

(в)-Пз

■ cos(5) • Z4 • G + cos(0) • h3

• 008'

(8К4 •l

6- X2

cos

(8) • h3 • Z 4 • G + cos(d) • X4 • G • sin(0) + h • Z 4 -G- cos2 (g)

6 • X2 • cos(q)

3 • cos(q)

6 • X2 • cos(0> cos(8)

cos

(0) • h3 • Z4 • G • cos2 (g) h3 • Z4 • G • sin(q) • cos(g) • cos(b)

6 • X2 • cos(8) cos(0)^3 •Z4 • G • cos2 (b) X4

3 • X2 • cos(8) • cos2 (g) • G • sin(0)

6 • X2 • cos(8)

3 • cos(0) • cos(8)

+

(36)

Z4 • G • cos(g) • cos(b) h3 • G • sin(0) • cos(g)

3^cos(8)

3 cos(8) cos^h •G^cos1

3 cos(8) cos(g) • G • X

(b). X2 •cos(8)^Z4^

6^ cos^)^ •^(0)^3

6•cos(0)•hз ^cos(d)

3 • 008(5) 6 • С4 • 008(5)

Из рассмотренного следует, что методика (2)-(34) позволяет получить матрицу жесткости произвольного конечного элемента с целью конечноэлементной аппроксимации машиностроительных конструкций различного назначения.

+

+

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

Конечноэлементная модель зубчатого зацепления разработана также для расчета остаточного ресурса портального крана типа «Ганц 16/22,5».

СПИСОК ЛИТЕРАТУРЫ

1. GallagerR. H., Padlog I., BijlardP. P. Stress Analysis in Heated Complex Shapes // J. Aerospace Sci. - 1962. - № 29. - P. 700-707.

2. Панасенко Н. Н., Левин А. И., Юзиков В. П. Расчет на сейсмические нагрузки машиностроительных конструкций из тонкостенных стержней // Изв. Северо-Кавказ. центра высш. шк. Техн. науки. - 1988. - № 3. - С. 75-82.

3. Панасенко Н. Н., Левин А. И., Юзиков В. П. Статический деформационный расчет пространственных тонкостенных стержневых систем произвольного вида // Изв. Северо-Кавказ. центра высш. шк. Техн. науки. - 1988. - № 4. -С. 134-138.

4. AnSys Theory Reference. 001242. Eleventh Edition. - SAS IP, Inc. 1999.