Алгоритм учета смещения как жесткого целого треугольного элемента оболочки при использовании множителей Лагранжа Текст научной статьи по специальности «Физика»

АЛГОРИТМ УЧЕТА СМЕЩЕНИЯ КАК ЖЕСТКОГО ЦЕЛОГО ТРЕУГОЛЬНОГО ЭЛЕМЕНТА ОБОЛОЧКИ ПРИ ИСПОЛЬЗОВАНИИ МНОЖИТЕЛЕЙ ЛАГРАНЖА

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

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

При осуществлении конечно-элементного анализа оболочек неизбежно возникает весьма сложная проблема по учету смещений используемого элемента дискретизации как жесткого целого. На наличие данной проблемы указывают как отечественные [1, 2], так и зарубежные [3, 4] исследователи. Имеющиеся к настоящему времени предложения по решению указанной проблемы [5, 6, 7] сводятся к попыткам учесть жесткие смещения элемента дискретизации в явном виде за счет включения в интерполяционные полиномы дополнительных членов, расширения матрицы жесткости путем конгруэнтного преобразования, строгого соблюдения совместности между соседними элементами и так далее. При этом приводились примеры расчета оболочки несложной геометрии, как правило, цилиндрической панели. В то же время в [8] предложено решение проблемы учета смещений конечного элемента как жесткого целого в неявной форме за счет использования векторного способа аппроксимации полей перемещений.

Целью настоящей работы является рассмотрение вопроса применения векторного способа аппроксимации перемещений [8] в алгоритмах расчета оболочек вращения, разработанных на основе треугольных элементов дискретизации, матрицы жесткости которых формировались с использованием корректирующих множителей Лагранжа, существенно повышающий точность конечно-элементных решений оболочек [9].

Рассмотрим элемент дискретизации в форме треугольного фрагмента, произвольно ориентированного на срединной поверхности оболочки вращения (рис.1) с узлами i, /, k . Столбец узловых варьируемых параметров данного элемента в глобальной системе координат (длина дуги меридиана) и в (угол, отсчитываемый от образующей против часовой стрелки) выбирается в виде

К г

1x27

(1)

где

{Г) ( I / к I / к i /к

Чу } = \ЧЧ Ч Ч5 Ч5 Ч 5 Ч в Чв Чв

ВГ Г КГ }'КГ }>

1x9 1x9 1x9 J

}. Здесь под ч понимаются тангенциаль-

ные V , V или нормальная V компоненты вектора перемещения.

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

границе соседних элементов дискретизации. Для устранения указанной рис 1 погрешности предлагается применять

1x9

корректирующие множители Лагранжа Л в дополнительных узлах 1, 2, 3, расположенных в серединах сторон треугольного конечного элемента (рис. 1):

(

Л

(т)

(т)

(т)

Л

дБ

V пт

дБ ,

пт

= о,

(2)

где верхний индекс т обозначает номер вышеупомянутого дополнительного срединного узла. Для отдельного треугольного конечного элемента (КЭ) равенство (2) может быть преобразовано к виду [9]

Л

(1)

5у

(1)

+ Л( 2 >--+ Л( 3 -= 0.

ЗSn дБп дБп

п1 п2 п3

( 2 )

>( 3 )

¿V

(3)

(3)

Входящие в (3) производные тУдБп могут быть выражены через производные нормальной компоненты вектора перемещения по глобальным криволинейным координатам

дv(m) дь>(т) (т, дv(m) -=-cos а + —

дБ,,

дБ

дв

cos Р

(т)

(4)

где а и Р - углы между ортом нормали пт и векторами локального базиса а°

- о

и а2, касательными к срединной поверхности в узле т.

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

дv дБ дv дв

(к-Г |+КГ 2]К}'

И |+<>Т £)(*} •

(5)

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

1x9

полные полиномы третьей степени [10]; £ п - локальные координаты, изме-

няющиеся от 0 до 1; } = {

, i ] к i ] к = ^VV V V,^ V,^ V,^ V ,„ V ,„ V ,„ }-

к i ] к \ г

VVV,ц } - столбец узловых значе-

1x9

ний нормальной компоненты вектора перемещения в локальной £ п системе координат. Равенство (3) с учетом (4) и (5) может быть записано в матричном виде [9]:

{Л}Г [} = о, где К} = ]{У}' {.у }' ^}

1x9

1x3 3x27

27x1

1x27

{V?

(6)

- столбец узловых неизвестных треугольного КЭ в локальной системе координат;

где

{Я}Т ,{/2}Т ,{/3}Т

1x9 1x9 1x9

{0}Т {0}Т {/1}Т

1x9 1x9 1x9

[^]= {0}Т {0}Т {/2 }Т

3x27 1x9 1x9 1x9

{0}Т {0}Т {/}

_ 1x9 1x9 1x9

-матрицы-строки, содержащие функции-формы, определяемые согласно (4), (5). 42

+

1x9

1x9

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

К Ыл К Г

(7)

27x1

27x1

где ] - матрица преобразования.

27x27

(8)

С учетом (7) равенство (6) может быть записано в виде

И' ['][',К} = М' [т]{иуГ}- 0.

1x3 3x27 27x27 2М 1x3 3x27 27к1

Функционал, выражающий равенство энергии деформации для треугольного КЭ с учетом (8) запишется следующим образом

* - К}' [КК} - {иуГ}' {р} + {я}' [']{иуГ} - 0. (9)

Минимизируя (9) по узловым неизвестным {иу } и значениям корректирующих множителей Лагранжа в дополнительных узлах на серединах на серединах сторон треугольного КЭ {X}', можно получить систему уравнений

дф/д{иГ}' -[К]{иуГ}-{р} + [']' {х}-0;

д*дЩт -[']{иУ}-0.

(10)

Систему (10) можно п ровке

редставить в удобной конечно-элементной формули-

(11)

где

[ к, }

30x30

[К] [Т]т

27x27 27x3 ['] [0]

Кг к }-{рр г

0x30 30x1 30x1

; (и;Г --К}Т {х}г}; [Рг}' --{{Р}' {0}Т}.

1x30 I 1x27 1x3 } 1x30 11x27 1x3 }

Для того чтобы в описанном выше алгоритме реализовать векторный способ аппроксимации перемещений [8], равенство (3) необходимо записать в векторной форме

X1 — + Я2 — + Я3 — - 0, где Я(т> - Х(т>а0(т> + Х2(т>а20(т> + Х(т>а0(т>; (12)

дБ,,

дБ,,

дБ,,

п3

дБ

дБ

( т >

- cos а +

дв

cos Р

(т >

(12а)

Входящие в (12а) производные вектора перемещения по глобальным криволинейным координатам Б и в определяются равенствами

дУ

дБ

т

дУ дв

,1(т> ~0(т > ,2(т> ~0(т> .(т > ~0(т >

-11 а^ + ^ а2 + ц а ;

- 1 т >а 0( т > +12( т >а 0( т > +1(т >а 0( т >

(13)

Л( т > ,2( т > .( т > Л( т > ,2( т > .( т >

где ^ , tl , tl , 12 , 12 , t2

многочлены, содержащие компоненты век-

тора перемещения и их производные в узле т, расположенном в середине стороны треугольного КЭ [8], например,

*(т> - у(т> +и°(т>,Мт>. *(т> _ у(т> +и°(т>,,2(т> (14)

<1 - У,б V , <2 - у'в 22 У ' V14'/

3x27

3x3

2

тт

ду ду

(15а)

Здесь b°, b°2 - компоненты тензора кривизны. Для того чтобы размерность

расширенной матрицы жесткости Кр не превышала порядка 3° х 3° при выполнении скалярного умножения в (12) можно ограничиться нормальной компонентой множителя Лагранжа, т.е. воспользоваться равенством

rhl<1> гК<2> г?<3>

1< 1 }a0<1 > • — + х<2}a°<2> • — + х<3}a0<3> • — = °, (15)

dS„ dSn dSn

«1 n2 «3

или с учетом (12а) и (13)

х<1; (41; cos «<1;+41; cos р< 1;)+А< 2; (42; cos «<2; + 42; cos Р< 2;)+

+Х(3; [4Ъ>cos а<3; +123>cos р<3; ) =

Входящие в (15а) многочлены t[m), t(2m> (m = 1,2,3) определяются на основе формул векторной аппроксимации перемещений [8]

г = ыт №} = ыт №г} = ыт ИИ [ N] {К} =

1х9 9х1 1х9 9х9 9х1 1х9 9х9 9х27 27х27 9х1

= мт №] [n ]{<},

1х9 9х27 27х27 27х27 9х1

где {V/ }Т =\vlv]vkvl^v^v^vlnv,]nv!:n}; }Т = {v'V]VkV!SVjVksVjeV,jVвв}

1х9 1х9

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

[l] [A] = [A] [g] . Здесь [A] - квазидиагональная матрица, содержащая локаль-

9х9 9х27 9х27 27х27

ные векторы базисов узловых точек треугольного КЭ , а°!, а°, 1,..., а°k , которые могут быть выражены через векторы локального базиса внутренней

а1 , а2 , а

{аot} = [Y]{а0}, где {аot}Т = {a1°ta2°ta°t}; {а0}Т ={а°а°°а°}; (17)

3х1 3х3 3х1 1х3 1х3

верхний индекс t обозначает узлы треугольного КЭ i, j, k . С учетом(17) можно сформировать равенство [A] = а1° [A1 ] + [A2 ] + а° [A3 ]. (17а)

Выражение (16) с учетом (17а) запишется следующим образом

I? = {<р}Т | а° [A1 ]+ [A2 ]+ а ° [A3 ]| [g] {Zy }, где {Zy }=[ N ]{u^ } . (18)

1х9 V 9х27 9х27 9х27 J 27х27 27х1 27х27

Дифференцируя (18) по глобальным криволинейным координатам S ив, можно получить производные вектора перемещения

,1 - ° 2 ? ° , ? ° V ,S = t1 а1 +11 а2 + ^а =

к } § + KF § J(«-1° [A1 М [A2] + а0 [A3])[g]{Zy};

f} ,1 ? О ,2 - 0 , ? ° V,в = t2(X1 + t2 а2 + t2(X =

{^}Т I + W SJ^[A1 ] + «2°[A2] + а°[A3])[G]{Zy}.

Из соотношений (19) можно получить зависимости для многочленов

(т) (т)

'1 >'2

'1 =1кГ 5+кГ Аз 11° К2-}'

£ I - \г>п ( - I ^з

ч ^ дS 1 ц дS Г 3 } д д\ (2°) =}{^}Т ^+{<р'ч}} £][ Аз }

Равенство (15а) с учетом (2°) может быть записано в матричном виде

1МТ

1x27

{1}Т [в 1«} = 0, (21) где [В]=&}'. (21а)

1x3 3x27 27х1 3x27 1x27

(Ьз}Т

\1x27 _

Входящие в (21а) матрицы-строки Ьт определяются согласно (15)...(20), например,

\Т д£ ( дц

М'

1=0.5 ц=0.0

}{ д^ { }Т дц] (1)

!{р'*} дs + {р,Ц} дs)с°'а +

^ дS 1 ц дS)

(22)

[ А31[° 1{2у }.

(КГ I +мТ 5 ] с-'1)

Сопоставляя равенства (6), (7) и (21), (21а), следует отметить весьма существенные отличия данных соотношений, которые обусловлены принципиально различными способами аппроксимации полей перемещений в треугольном КЭ. Так, например, матрица [F] содержит только третью часть ненулевых элементов, в то время как матрица [В] является полностью заполненной. Этот факт объясняется тем, что в соответствии с традиционной интерполяционной процедурой каждая компонента вектора перемещения и ее производные зависят от узловых значений только этой же компоненты и ее производных, как это видно из (5). В соответствии с векторной аппроксимацией полей перемещений [8] отдельная компонента вектора перемещения и ее производные зависят от полного

набора узловых варьируемых параметров треугольного КЭ {и^}, в структуру

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

Дальнейшая процедура получения расширенной матрицы жесткости [Кр] и столбца узловых нагрузок {Рр} при векторном способе аппроксимации аналогична (9)___(11) с заменой матрицы \Т] на [В] и соответствующей заменой интерполяционных формул, входящих в матрицы [К] и {Р} [8].

Пример расчета. Была рассчитана открытая с торцов оболочка в форме эллипсоида вращения, загруженная внутренним давлением интенсивности q, имеющая на левом краю пружинные опоры, позволяющие смещаться ей под действием заданной нагрузки в меридиональном направлении как жесткому телу (рис. 2). Радиус вращения вычислялся по формуле

Ь Г~2 2 г = — у/а - х

а

где Ь = 0,9 м; а = 1,3м; осевая координата х изменялась в пределах 0 < х < 1.2 м.

Были приняты следующие исходные данные: q = 5МПа; Е = 2 • 105МПа; V = 0,3; ' = 0,02м. Расчеты были выполнены в двух вариантах: в первом варианте был реализован алгоритм, описанный соотношениями (2)_(11); во втором варианте при формировании матрицы жесткости треугольного КЭ с множителями Ла-гранжа использовалась векторная аппроксимация перемещений (12)_(22).

Вследствие наличия осевой симметрии оболочка была представлена одной лентой треугольных КЭ, ориентированных в меридиональном направлении. _Таблица 1

Величина смещения оболочки как жесткого тела, м Вариант 1

Точка 1 Точка 2

в а —, МПа н 7 ам в а , МПа н ак в а —, МПа н ам в а —^, МПа н ак

0.00 0.02 169.2 0.76 171.6

0.71 167.2 0.09 164.9

0.09 6.2 58.5 7.1 66.1

-20.6 46.2 -24.4 37.7

0.90 62.2 -937.8 64.3 -883.1

-212.6 1042.7 -245.1 -1106.1

9.00 618.5 -10840.8 632.7 -10318.0

-2121.2 -11866.3 -2438.7 -12476.0

Величина смещения оболочки как жесткого тела, м Вариант 2

Точка 1 Точка 2

в а —, МПа н ам в а —^, МПа н ак в а —, МПа н ам в а —^, МПа н ак

0.00 -0.17 169.2 0.81 171.6

0.87 167.1 0.50 164.9

0.09 -0.17 169.2 0.81 171.6

0.87 167.1 0.50 164.9

0.90 -0.17 169.2 0.81 171.6

0.87 167.1 0.50 164.9

9.00 -0.17 169.2 0.82 171.6

0.87 167.2 0.49 165.0

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

Верхние индексы «в» и «н» в обозначениях напряжений указывают на внутреннюю и наружную поверхности оболочки соответственно.

Анализ результатов показывает, что при отсутствии жесткого смещения численные значения напряжений практически совпадают. Однако с возникновением жесткого смещения результаты первого варианта нельзя признать удовлетворительными.

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

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

Л и т е р а т у р а

1. Голованов А.И., Тюленева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций. - М.: Физматлит, 2006. - 391 с.

2. Скопинский В.Н. Напряжения в пересекающихся оболочках. - М.: Физматлит,

2008. - 399 с.

3. Зенкевич О.М. Метод конечных элементов в технике. - М.: Мир, 1975. - 542 с.

4. Оден Дж. Конечные элементы в нелинейной механике сплошных сред : Перев. с англ. - М.: 1976. - 464 с.

5. Кантин Г. Смещение криволинейных элементов как жесткого целого // Ракетная техника и космонавтика. - 1970. - №7. - С. 84-88.

6. Кузнецов Ю.М., Голованов А.И. Элементы с явным выражением жестких смещений в расчете тонких цилиндрических оболочек // Прочность и устойчивость оболочек. Труды семинара. - Вып. XIX. Часть II. - Казань - 1986. - №2. - С.83-93.

7. Косицын С.Б. Метод построения базисных функций для искривленных конечных элементов с учетом жесткого смещения // Исследования по строительным конструкциям и их элементам. М.: ЦНИИСК. - 1982. - С. 17-27.

8. Николаев А.П., Клочков Ю.В., Киселев А.П., Гуреева Н.А. Расчет оболочек на основе МКЭ в двумерной постановке. - Волгоград: ИПК ФГОУ ВПО ВГСХА «Нива».

2009. - 196 с.

9. Клочков Ю.В., Николаев А.П., Вахнина О.В. Использование множителей Лагран-жа при формировании матрицы жесткости треугольного конечного элемента// Строительная механика инженерных конструкций и сооружений. - 2009. - № 2. - С. 38-43.

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

THE ALGORITHM OF THE ACCOUNT OF DISPLACEMENT AS A RIGID WHOLE TRIANGULAR SHELL ELEMENT WITH THE HELP OF LAGRANGE MULTIPLIERS

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

A question of application of a vector way of approximation of movings in the algorithms of calculation of environments of rotation developed on the basis of triangular elements of digitization which matrixes of rigidity were formed with use of the correcting multipliers of Lagrange essentially raising accuracy of certainly-element decisions of environments