Использование множителей Лагранжа при формировании матрицы жесткости треугольного конечного элемента Текст научной статьи по специальности «Физика»

Численные методы расчета конструкций

ИСПОЛЬЗОВАНИЕ МНОЖИТЕЛЕИ ЛАГРАНЖА ПРИ ФОРМИРОВАНИИ МАТРИЦЫ ЖЕСТКОСТИ ТРЕУГОЛЬНОГО КОНЕЧНОГО ЭЛЕМЕНТА

Ю.В. КЛОЧКОВ, д-р техн. наук, профессор, А.П. НИКОЛАЕВ, д-р техн. наук, профессор, О.В. ВАХНИНА, ст.пр.

Волгоградская государственная сельскохозяйственная академия

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

Рассмотрим треугольный КЭ с катетами единичной длины в локальной системе координат 0 < 4, ^ < 1 (Рис.1) с узлами г, ],k , на который отображается произвольный треугольный элемент оболочки вращения.

Столбец узловых варьируемых параметров данного КЭ выбирается в виде [2]

к Г=1к Г к Г к Г1 (1)

4

где и и V - тангенциальные, а w- нормальная компоненты вектора перемещения. Верхний индекс «л» указывает, что рис. 1 столбец узловых неизвестных определен в локальной сис-

теме координат.

Входящие в правую часть равенства (1) подматрицы-строки имеют следующую структуру

к г 4

г ! кг г к г г к I

в 9 9

(2)

Ъ

где под в понимается компонента вектора перемещения и, V или w.

Введем в середину каждой из сторон треугольного КЭ дополнительные узлы, которые обозначим цифрами 1,2,3 и укажем в них направления внешних

нормалей (т = 1,2,3 ) (рис.2).

Рассмотрим в качестве варьируемых параметров во вновь введенных узлах 1,2,3 производные нормальной компоненты вектора перемещения по направлениям внешних нормалей ст/(т = 1,2,3 ) к сторонам треугольного КЭ.

Дополнительные узлы 1, 2, 3 треугольного КЭ (например, имеющего номер I) в то же время являются узлами для соседних элементов сетки дискретизации (рис.3).

Таким образом, например, для дополнительного узла 2 будет справедливо равенство 5 т (1) ,5 т (Щ = о (3)

582 С37 = .

0

Рис.2

X

z

0

Рис.3

X

Для отдельного треугольного КЭ на основании (3) можно записать равенство, которое можно рассматривать как дополнительное условие, необходимое для построения функционала Лагранжа

<4>

где X 2, Х3 - множители Лагранжа.

При формировании матриц жестко-стей соседних по сетке дискретизации треугольных КЭ с учетом направлений

нормалей Sn будут реализованы соотношения для всех сторон отдельного КЭ (4). Производные нормальной компоненты вектора перемещения в направлении нормалей к сторонам треугольного КЭ могут быть выражены через обычный столбец узловых варьируемых параметров (1)

^ Г и № Г & ) (5)

О \ 1x27 27x1 1x27 27x1

где & |т = | {иу |т {^ |т |т I - столбец узловых неизвестных в глобальной систе-

у )

1x27

у >

1x9

у 1x9

у 1x9

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

Правую часть равенства (4) можно представить в матричном виде

мт

О w

О 8П

О w

О 8П

О w

О 8П

3x1

:МТ

и)

1x27

{г2)1{иуу}=Мт [G]{иу}.

1x27

{г1)

(6)

Функционал, выражающий равенство потенциальной энергии деформации на возможных перемещениях для треугольного КЭ, с использованием множителей Лагранжа можно записать в следующем виде

П = } КР )dV+ {иу } Мт {!)-{ {и)т {р)<* = 0, (7)

V F

где {^р)Т = {^8^28^ ) - матрица-строка деформаций в произвольном слое оболочки, отстоящем от срединной поверхности на расстоянии ; ; {стар)=(ст11ст22ст12) -

столбец напряжений в произвольном слое оболочки; {и)Т ={uvw)- матрица-строка, содержащая компоненты вектора перемещения внутренней точки треугольного КЭ; {р)={р;р2р3)- столбец внешней поверхностной нагрузки треугольного КЭ.

Минимизируя функционал (7) по узловым неизвестным {иу )Т и множителям Лагранжа {х)т, получим следующие равенства

О П

О{иу )Т

[К]{иуу )+ [0]Т {*)-£ )= 0;

О П

[°]{иу}

= 0. (8)

Входящие в (8) матрица [К] и столбец Г} определяются стандартным для конечно-элементной процедуры образом [3]

[К] =[Рк Г {[В]т [г]т [с][г]№у [ра ]; {Г} = [РК ]т {[л]т [р]

27x27 27x27 27x6 6x3 3x3 3x6 6x27 27x27 27x1 27x27 г27х3 3x1

где [Рк ] - матрица преобразований при переходе от локальной системы коор-

27x27

динат 4, П к глобальной, например, S, 0 .

Принимая во внимание (8), можно сформировать необходимое для конечно-элементной процедуры матричное равенство

[K] [g]t

tGt 1

3x27 3x3

к }

27x1

м

27x1 {0}

3x1 ,

или [K]p {иуг } ={f }р, (9)

30x1

'р'

30x1

где [K]p, {иу } - расширенные матрица жесткости и столбец узловых неизвестных треугольного КЭ.

Процедуру определения входящих в (5) матриц-строк {tm} можно продемонстрировать на примере оболочки вращения, срединная поверхность которой в исходном состоянии описывается радиус-вектором

R0 = xi + r sin в j + r cos ek, (10)

где x- осевая координата, r- радиус вращения.

Формула, устанавливающая связь между глобальными координатами S, 6 и локальными координатами £ и q треугольного КЭ может быть записана в

виде [2] S = (1 -^-n)Si + 4Sj + nSk; 0 = (i- 4-n)0' + 40j + n0k . (11)

Уравнение гипотенузы треугольного КЭ в локальной системе координат имеет вид

q = 1-£. (12)

В результате подстановки (12) в (11) можно получить равенства

S = ^(Sj - Sk)+ Sk; 0 = ^(0J - 0k)+ 0k. (13)

Выразим из 1-го уравнения (13) локальную переменную £ и подставим ее во второе уравнение (13). В результате получим следующую зависимость глобальной координаты 6 от дуги меридиана S на границе j - k треугольного КЭ

(s - sk\ej -вк)+вк (sj - Sk)

Sj - Sk

Орт касательной к стороне j - k треугольного КЭ можно получить дифференцированием (10) по формуле полной производной, т.к. на этой линии 6 = f (s)

в = Л-Л 'k V-¿. (14)

8F 8F дв - . „г í

т

t = R =--1---= x i + r sine j + r cosek + Ir cose j - r sinek -в , (15)

,S 8S дв 8S 'S ,S ,S v ' ,S

a 0j-0k где 0 =

S SJ - Sk '

Орт нормали к стороне ] - k треугольного КЭ определим в результате векторного произведения

^=а0 x т, (16)

где а0 - орт нормали к срединной поверхности оболочки вращения

a0 =-r i + x sin 6 j + x cos 6 k.

,S ,s J ,s

3x1

Определенный по формуле (16) вектор находится в плоскости векторов

локального базиса а0, а0, касательных к срединной поверх- Да ности оболочки (рис.4).

Косинусы углов а и р (рис.4) могут быть вычислены при использовании формулы скалярного произведения

„О с2 „О с2

cosa =

■ n -; cosp = - n

а0 sn а2 sn

(17)

2

sn

Р

a;

Рис.4

Рассмотрим сетку дискретизации оболочки вращения в виде (рис.5)

S

Рис. 5

Можно убедиться, что для треугольного КЭ под номером I будут справедливы следующие соотношения

дм дм дм дм дм дм дм

Для треугольного КЭ, ориентированного подобно элементу под номером II, соотношения (18) принимаются с противоположным знаком.

Таким образом, матрицы-строки {^ ] для треугольного КЭ под номером I будут иметь вид

=--cosan---cos В. (18)

dS2 8S гдв

w = ыт {0}Т {-ф н {t3 } =]{0}Т {0}Т u i1

1x27 [1x9 1x9 1х9 [ 1x27 [1x9 1x9 1х9

{t2} =^{0}Т {0}Т ts}Tcosa + {ф,е}Т cosр)[,

1x27 ' ~

(19)

где {ф}Т ={ф 1 ф 2... ф 9}- матрица-строка функций формы, определенных в соот-

1x9

ветствии [2]. Входящие в (19) производные полиномиальных функций определяются по формулам

^ ^щ-ВД (20) ф}=Ф г I-к Г §.

В качестве примера была решена задача определения напряженно-деформированного состояния (НДС) жестко защемленного цилиндра, нагруженного внутренним давлением интенсивности q (рис.6). Были приняты следующие исходные данные:

L = 1.0 м; R = 1.0 м; Е = 2 -105МПа; V = 0.3; t = 0.02 м; д = 5МПа.

Вследствие наличия осевой симметрии рассчитывалась 1/4 часть оболочки.

1x9

Рис. 6

Расчеты были выполнены в двух вариантах: в первом варианте в качестве элементов дискретизации использовались треугольные КЭ, матрицы жесткости которых формировались стандартным образом [1,2]; во втором варианте был реализован описанный выше алгоритм, основанный на использовании множителей Лагран-жа (3)...(19).

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

Число узлов сетки дискретизации вдоль кольца было принято равным 4, а вдоль образующей варьировалось от 5 до 57.

Анализ полученных результатов показывает, что в первом варианте наблюдается неудовлетворительная сходимость вычислительного процесса, особенно в жесткой заделке (точка 1). Кроме того, в наружной и внутренней поверхностях оболочки в жесткой заделке с измельчением сетки дискретизации напряжения имеют одинаковый знак, что противоречит физическому смыслу решаемой задачи, т.к. в жесткой заделке под действием внутреннего давления возникает деформация изгиба, т.е. внутренняя поверхность растягивается, а наружная сжимается, что и наблюдается во втором варианте расчета. Анализируя численные значения напряжений, представленные в таблице № 2, можно отметить быструю сходимость конечно-элементных решений уже при достаточно редкой сетке дискретизации. Так при размере сетки 4 х17 погрешность вычислений находится в пределах 1%. Полученные значения напряжений во втором варианте хорошо согласуются с физическим смыслом решаемой задачи.

Таблица 1

Номера точек Численные значения напряжений, МПа Число узлов сетки дискретизации

4 х 5 4 х 17 4 х 33 4 х 49 4 х 57

1 ст в 333.71 149.91 108.06 93.14 88.80

стн -210.33 -24.94 17.22 31.99 36.26

ст ср 61.69 62.48 62.64 62.56 62.53

2 ст в 47.24 57.33 58.64 62.14 62.19

стн 69.08 68.81 67.87 61.34 61.53

ст ср 58.16 63.07 63.26 61.74 61.86

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

неизвестных в виде [3] {и^ }Т =|{ил }Т ^ }Т|, где }Т ^Уд^ }. (21)

1х8 I 1х4 1х4 \

Здесь под д понимается меридиональная или нормальная компонента вектора перемещения, а локальная координата х изменяется в пределах -1 < х < 1.

Результаты конечно-элементного решения рассматриваемой оболочки с использованием в качестве элемента дискретизации одномерного КЭ с размером матрицы жесткости 8 х 8 представлены в табл. № 3, структура которой аналогична предыдущим табл. 1 и 2.

Номера точек Численные значения напряжений, МПа Число узлов сетки дискретизации

4 х 5 4 х 17 4 х 33 4 х 49 4 х 57

1 в 408.16 475.49 479.62 480.40 480.57

СТн -286.17 -356.06 -360.26 -361.05 -361.21

ср 61.00 59.72 59.68 59.68 59.68

2 в 69.62 67.15 67.05 67.03 67.02

СТн 49.65 52.20 52.31 52.33 52.33

ср 59.63 59.68 59.68 59.68 59.68

Таблица 3

Номера точек Численные значения напряжений, МПа Число узлов сетки дискретизации

1х5 1х 17 1х 33 1х 49 1х 57

1 в 368.38 475.52 479.62 480.39 480.56

СТн -241.19 -356.01 -360.22 -361.01 -361.17

ср 63.24 59.76 59.70 59.69 59.69

2 в 73.04 68.01 67.90 67.89 67.88

СТн 46.18 51.38 51.49 51.50 51.51

ср 59.60 59.69 59.69 59.69 59.69

Как видно из табл. № 3, численные значения контролируемых параметров НДС практически совпадают с данными таблицы № 2, начиная с сетки узлов 1 х 17 , что позволяет сделать вывод о достоверности полученных результатов.

Основываясь на анализе табличного материала, представленного в таблицах 1,2,3 можно сделать окончательный вывод о высокой эффективности разработанного алгоритма формирования матрицы жесткости треугольного КЭ с использованием множителей Лагранжа.

Л и т е р а т у р а

1. Клочков Ю.В. Сравнительный анализ эффективности использования конечных элементов треугольной и четырехугольной форм в расчетах оболочек вращения// Ю.В. Клочков, А.П. Николаев, Н.А. Гуреева// Изв. вузов. Строительство. - 2004. - № 3. -С.103-109.

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

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

LAGRANGIAN COEFFICIENTS USING IN CASING CALCULATION WITH TRIANGULAR FINITE ELEMENT

Y.V.Klochkov, A.P.Nikolaev, O.V.Vakhnina

In this article, the calculation of hard pinched cylinder casing mode of deformation with one-dimensional and triangular finite element using was presented. Comparative analysis of gotten results was done/ Calculations improvement while using triangular finite element as discretization element method with Lagrangian coefficient using was offered.