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

Вычислительные технологии

Том 16, № 1, 2011

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

Ю.В. Клочков, А.П. Николаев, О. В. Вахнина Волгоградская государственная сельскохозяйственная академия, Россия e-mail: [email protected], [email protected]

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

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

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

В [1, 2] без приведения численных результатов указана возможность обеспечения совместности значений производных нормального перемещения на границах треугольных КЭ. Для решения этой задачи применяются метод корректирующих множителей Лагранжа и метод штрафа [1]. В первом случае к функционалу Лагранжа добавлены слагаемые, представляющие интегралы от произведения функции Лагранжа на разность алгебраических выражений производных прогибов соседних элементов по их границам. Метод штрафа предполагает применение в функционале Лагранжа штрафного слагаемого с единым множителем, регулирующим совместность деформирования на границах соседних элементов [1].

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

Конечный элемент выбирается в виде криволинейного треугольного фрагмента срединной поверхности оболочки вращения с узлами i,j, k. Столбец узловых варьируемых параметров рассматриваемого треугольного элемента дискретизации в глобальной системе координат содержит компоненты вектора перемещения и их первые производные

по криволинейным координатам £ (длина дуги меридиана) и 9 (угол, отсчитываемый от образующей против хода часовой стрелки) [4]

{Щ} = ^ к}тк}тк}т

1x27 1x9 1x9 1x9

(1)

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

К У = № 4 як я,8 я,8 я,к8 я,в я,в я,в}

19

(2)

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

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

дш® дш^

дБ3

к

+

дБ3':

- = 0,

(3)

где верхние индексы I, II указывают на номера смежных элементов дискретизации.

Однако равенство (3) в силу несовместности по производным компонент вектора перемещения не выполняется, поэтому для его корректного соблюдения рассмотрим интегральное равенство

/

Х

{дБ'п

дш дш

+

к

дБЗ,

-к'

Ж-к = 0.

(4)

Здесь р к длина дуги ]—к дискретного элемента; ХР к функция множителя Лагранжа в произвольной точке дуги ] — к; №-к дифференциал дуги ] — к.

к

Рис. 1. Смежные треугольные элементы дискретизации в глобальной системе координат

Рис. 2. Треугольный конечный элемент в локальной системе координат

Для отдельного треугольного конечного элемента равенство (4) может быть трансформировано к виду

У д&п3 У д&пк У 1 }

¡г-] Ц-к ¡к-г

где Хг-* Х-к и Хк-г — функции множителей Лагранжа в произвольных точках соответствующих сторон дискретного элемента.

Рассмотрим треугольный элемент дискретизации в локальной системе координат 0 < £ + п < 1 (рис. 2). Связь между глобальными Б, 9 и локальными п координатами осуществляется с помощью соотношений

Б = (1 - £ - п)Бг + £Б* + пБк, 9 = (1 - £ - п)9г + £9* + п9к. (6)

Множители Лагранжа на границах треугольного КЭ с учетом (6) могут быть выражены через их узловые значения

Х-* (п = 0) = (1 - £)Хг + , Х-к (£ = 1 - п) = (1 - п)Х + пХк,

Хк-г(£ = 0) = (1 - п)Х + пХк. (7)

Соотношения (7) могут быть представлены в матричном виде

{Ат-п}Т =[7(£,п)] (Хг}, (8)

3x1 3x3 3x1

где

1 - £ £ о

{Хт-П}Т = {Х-х*-кхк-г}, Х}т = (ХгХХк}, [7(£, п)] = ( 0 1 - п п

1 - п 0 п

1x3 1x3 3x3

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

дк дк дк

дБП-* у дБП-к У дБк-г

где f1(£), /2(п), /3(п) _ функции, зависящие от координат на соответствующих сторонах локального треугольника. Входящий в (9) столбец узловых неизвестных в локальной системе координат имеет следующую структуру:

и }т = { {< }т {< }т К П, (10)

1x27 ^ 1x9 1x9 1x9

где {$}Т = {дг о> дк д^ д,к^ д^ д^ дк }.

1x9

Соотношение (5) с учетом (8) и (9) может быть представлено в матричном виде

Х }т [7 (£,п)]т {/п }{Щ} = 0, (11)

1x3 3x3 3x27 27x1

где

Ш

3x27

/1 (£) ¡2 (п) /з (п)

Столбцы } и } представляют собой набор узловых констант для каждого треугольного элемента дискретизации. Матричное произведение [7(£, п)]т(/п} зависит от локальных координат (£ ми п) точек на сторонах треугольного КЭ, выбранных в качестве точек интегрирования.

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

кам:

1- у (1 - £)Ш) ЯГ3 + ] (1 - п)/з(п) ¿1

Ц-] ¡к-г

2 - I ) ¿Г- +1(1 — п)/2(п) №-к

¡г-] ¡]-к

3 / ^ Г \ Л]-к I / (г„\ лк-г

к—г

П/2(п) ¿Р-к + п/з(п) Л-

(12)

¡]-к ¡к-г (12)

мену переменных £ и п через тараметр —1 < Ь < 1 по формулам

£=-¿4—, т? =4—, с1£=-сИ, ¿71 = -(Й, 2 2 2 2 2 2

(13)

а затем воспользоваться квадратурой Гаусса

/ (Ь) ¿1 / (1г) • Щ,

(14)

-1

г=1

где — координаты и веса квадратуры Гаусса, п — количество точек интегрирова-

ния.

Дифференциал дуги, например, стороны I — ] дискретного элемента определяется выражением [4]

( ¿Б

¿0

где

сИ 2 сИ 2 ;

Равенство (11) с учетом (12) и (13) можно представить в виде

(Ау}Т [Р] (Щ} = 0,

1x3 3x27 27X1

где

[Р ]

3x27

Г №} 1

1x27

< ш >

1x27

№}

^ 1x27 )

(15)

(16)

1

п

2

2

Входящие в (16) подматрицы-строки {Тп} (п = 1, 2, 3) определяются согласно (13) (14),

1x27

Равенство (16) может быть записано в следующем виде:

{\}Т [Ря] и} = {^}т [Я] и} = 0, (17)

1x3 3x2727x27 27x1 1x3 3x27 27x1

где [Ря] — матрица преобразования, формируемая на основе соотношений, получаемых

27x27

дифференцированием (6) по локальным координатам £ и ц:

дБ/д£ = Б3 - Б\ дБ/дп = Бк - Бг,

дв/д£ = в3 - вг, дв/дп = вк - вг. (18)

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

Если в соотношении (5) вместо функций корректирующих множителей Лагранжа Аг-3, А3-к, Ак-г использовать их значения А1, А2, А3 в дополнительных узлах 1, 2, 3, расположенных в серединах сторон треугольного элемента дискретизации (рис, 3), то полу чим выражение

А1 [ ^Г- + А2 [ + А3 [ = 0. (19)

¡г-] р-к 1к-г

Вычисление интегралов, входящих в (19), можно выполнить по формулам, аналогичным (12)—(14), с учетом того, что из подынтегральных выражений (12) необходимо исключить сомножители вида (1 - £), (1 - п)-, £> П- Для реализации данного варианта используется сетка дискретизации с узлами, расположенными в серединах сторон треугольных элементов дискретизации (рис, 4),

В каждом из основных узлов сетки дискретизации, отмеченных на рис, 4 кружками, имеется до девяти (в зависимости от заданных граничных условий) узловых неизвестных (2), В дополнительных узлах, расположенных в серединах сторон треугольных КЭ

(отмеченных квадратами), в качестве узловых варьируемых параметров принимаются

А1 А2 А3

Рис. 3. Треугольный элемент дискретизации с множителями Лагранжа в серединах сторон

Рис. 4. Сетка дискретизации

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

а1!^+л2!!1+л8!^0- <2о>

дБп1 дБп2 дБп3

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

1 2 3 1 2 3

столбец узловых неизвестных (10)

дьт

дБп

г 1Т д£ , , Т дп \

Ш, (21)

где ат и вт ~ углы между вектором пт и касательными векторами локального базиса а1 и а2 узла т (ш = 1, 2, 3), расположенного в середине стороны треугольного элемента на срединной поверхности оболочки,

( /п}

деляютея путем подстановки в интерполяционные полиномы третьего порядка [5] 123

Таким образом, процедура интегрирования по сторонам треугольного элемента дискретизации (12)—(15) исключается.

Выражение (20) с учетом (17) и (21) может быть записано в виде

(АС}Т [Ф] [Ря] (Щ} = (Ас}Т [У] (Щ}, (22)

ЧЧ = (Ас] ^^

1x3 3x2727x27 27x1 1x3 3x27 27x1

где

Г (Ь1}

(Ас}Т = (А1,А2, А3}, [Ф] = <^ (12}

3x27 1 (13}

Равенство работ внешних и внутренних сил на возможном перемещении для треугольного КЭ можно представить в виде выражения

Ф = и}Т [К] и} —и}Т (Я} +(Ас}Т [У] и} = 0. (23)

27x1 27x27 27x1 1x27 27x1 1x3 3x27 27x1

Варьируя функционал (23) то узловым неизвестным (Щ}Т и значениям множителей (Ас}Т

д(и} 27x27 27>а 27x1 27x3 3x1

<9Ф _ = (24)

Т 3x27 27x1

д{Хс}Т - М Ю = О-

(24)

вид

[К] (Щм} = Я}, (25)

30x30 30^ 30x1

т

[K ] [Y ] т

27x27 27x3

[Y] [0]

, {Щы}т = [Щ}т[Xc}T\, R}т = {R}T{0}T } расти-

1x30 1x27 1x3 1x30 L 1x27 1x3

[K ]

30x30

3x 27 3x 3

репная матрица жесткости треугольного КЭ, вектор искомых узловых неизвестных и расширенный вектор узловых усилий соответственно.

Процедуру определения входящих в правую часть (21) матриц-строк стоящих перед столбцом узловых варьируемых параметров [Щ} можно продемонстрировать па примере оболочки вращения, срединная поверхность которой в исходном состоянии описывается радиус-вектором

R = xi + r sin 9j + r cos 9k, (26)

где x — осевая координата, r = r(x) — радиус вращения. Дифференцированием (26) по координатам S, 9 определяются векторы локально го базиса а1; а2 с нормаль ю а = а1 х а2/ |а 1 х а2|.

Связь между глобальными координатами S, 9 и локальными координатами £ и п

(6)

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

П =1 - £• (27)

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

S = £ S - Sk) +Sk, 9 = £ (9j - 9k) + 9k• (28)

Выразим из первого уравнения (28) локальную переменную £ и подставим ее во второе (28)

ты 9 от дуги меридиана S на границе j — k треугольного КЭ:

е = (s -sk) (вз - ñ+ek -sk) ,29) Si - Sk ' 1 '

Орт касательной к гипотенузе треугольного КЭ можно получить дифференцированием (26) по формуле полной производной, так как на этой линии 9 = f (S)

т т д9

t = R,S = тг77 + = x,sf 1 + г,<? sill 9j +

dS д9 dS

+r,S cos 9k + (r cos 9j — r sin 9k) 9,S , (30)

9j — 9k

j — k

в результате векторного произведения

n = SП = a xt• (31)

Найденный по формуле (31) вектор S2n находится в плоскости векторов локального 12

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

а1 ^п о а2 ^п /00\

cosa = i i Ig2|> cos ¡3 = (32)

• fin1 21 • fin1

Для треугольного КЭ для узла 2 будет справедливо соотношение

дт дт дт

Таким образом, матрица-строка узла 2 для треугольного КЭ имеет вид

(33)

21 Т

О2}

1x27

{0}Т{0Г ( }т сов а + }Т сов в}

1x9 1x9 \ 1x9

Т

(34)

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

д£

Т

дп

дв'

дБ'

(35)

В качестве примера расчета была решена задача по определению напряженно-деформированного состояния цилиндра, загруженного двумя противоположно направленными сосредоточенными силами (рис. 5). Приняты следующие исходные данные: модуль упругости Е = 7.38 • 106 Н/см2, коэффициент Пуассона V = 0.3, радиус цилиндра Я = 12.58 см, толщина оболочки t = 0.24 см, сил а Р = 453.6 Н.

Расчеты были выполнены в трех вариантах. В первом случае реализован алгоритм с множителями Лагранжа в вершинах треугольного КЭ (соотношения (4)—(18)), во втором использовались множители Лагранжа в серединах сторон треугольника в сочетании с процедурой интегрирования по сторонам треугольного элемента дискретизации (соотношение (19)), в третьем варианте процедура интегрирования по границам треугольного элемента была исключена (соотношения (20)-(35)). В качестве контрольного реализован стандартный для треугольного КЭ алгоритм без использования корректирующих множителей Лагранжа [5]. В табл. 1 для всех вариантов приведены значения прогибов (в см) в точке приложения сосредоточенной силы (точка А на рис. 5) в зависимости от густоты сетки дискретизации, а также результат известного решения [6]. Анализ табличных данных показывает, что варианты I и II дают существенно заниженные (на несколько порядков) значения прогибов, поэтому их следует признать неудовлетворительными.

В варианте III наблюдается устойчивая сходимость вычислительного процесса, а численные значения прогиба практически совпадают с известным решением [6] (погреш-

Рис. 5. Расчетная схема цилиндра, нагруженного двумя сосредоточенными силами

Таблица 1. Значение прогиба в точке А приложения сосредоточенной силы P

Число узлов сетки Вариант расчета

дискретизации Контрольный I II III

5 х 5 -0.3058 -0.0130 -0.0118 -0.2715

8x8 -0.3061 -0.0092 -0.0045 -0.2770

10 х 10 -0.3033 -0.0074 -0.0028 -0.2773

11 х 11 -0.3021 -0.0068 -0.0023 -0.2773

Решение [6] -0.2883

ность при сетке узлов 11 х 11 составляет 3,8%), В контрольном варианте также наблюдается удовлетворительная сходимость вычислительного процесса, а значения прогиба

11 х 11

Таким образом, анализируя данные табл. 1, можно отметить, что варианты III и контрольный дают достаточно близкие величины прогибов в точке приложения силы P. Однако следует отметить, что расчеты па прочность требуют анализа не столько величин прогибов, сколько значений напряжений па внутренней и наружной поверхностях оболочки.

Значения напряжений в точке A приложения силы P могут быть вычислены с использованием узловых варьируемых параметров треугольного элемента как (I), так и элемента (II) (рис, 6), Результаты вычисления меридиональных ам и кольцевых ак на-

A

табл. 2, Дня сравнительного анализа численных значений напряжений были использованы вариант III с множителями Лаграпжа (соотношения (20)-(35)) и контрольный

P

численные с привлечением узловых неизвестных треугольного КЭ (I) и элемента (II), должны быть одинаковыми, что и наблюдается в варианте III расчета с множителями Лаграпжа, В контрольном варианте численные значения меридиональных напряжений, вычисленные с использованием узловых варьируемых параметров смежных треугольных КЭ, различаются в несколько раз, что указывает па неприемлемость погрешности расчета по стандартной конечно-элементной процедуре в треугольном КЭ без множителей Лаграпжа, Весьма интересным является анализ напряжений, вычисленных в точке B (рис, 5) оболочки, отстоящей от точки A на угол п/2 вдоль окружности. Численные значения напряжений в точке B приведены в табл. 3, в которой представлены значения

Р

Рис. 6. Сетка дискретизации цилиндра, на- Рис. 7. Деформирование цилиндра под дей-груженненч) сосредоточенными силами ствием сосредоточенных сил

Таблица 2. Значения напряжений в точке А приложения сосредоточенной силы P

Численные Вариант расчета

значения Контрольный III

напряжений, Число узлов сетки дискретизации

в точке А, 5x5 5x5 8x8 8x8 5x5 5x5 8x8 8x8

Н/см2 (I) (П) (I) (II) (I) (II) (I) (II)

в "м 5958 1313 11074 2565 2620 2626 2578 2583

е "м -7248 -2602 -13205 -4697 -1986 -1993 -2452 -2458

в U к 11427 12554 13139 15297 7660 7665 7172 7175

н UK -8671 -9798 -14353 -16513 -6049 -6054 -7226 -7230

Таблица 3. Значения напряжений в точке B

Численные Вариант расчета

значения Контрольный III

напряжений Число узлов сетки дискретизации

в точке В, 5x5 6x6 8x8 10 х 10 11 х 11 5x5 6x6 8x8 10 х 10 11 х 11

Н/см2

в "м 592 1098 1492 1605 1629 -1111 -1042 -952 -907 -894

е "м -916 -1338 -1680 -1773 -1789 1100 1089 1036 998 847

в U к -2290 -1725 -1228 -1045 -994 -4390 -4347 -4264 -4223 -4210

н UK 2179 1721 1276 1105 1060 3743 3827 3823 3801 3792

меридиональных и кольцевых напряжений во внутренней и наружной поверхностях оболочки в зависимости от густоты сетки дискретизации. Как следует из таблицы, меридиональные напряжения в контрольном варианте имеют противоположные знаки по сравнению с меридиональными напряжениями, полученными при использовании варианта III с множителями Лагранжа, Если проанализировать процесс деформирования оболочки под действием двух сосредоточенных сил (рис, 7), то можно отметить, что в точке B внутренняя поверхность будет сжата, а наружная растянута, что и наблюдается в варианте III с множителями Лагранжа, В контрольном же варианте адекватную картину напряженно-деформированного состояния получить не удалось.

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

Список литературы

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

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

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

[4] Погорелов A.B. Дифференциальная геометрия. М.: Наука, 1969. 252 с.

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

[6] Ashwe ll D.G., Sabir А.В. A new cylindrical shell finite elements based on simple independent strain function // Internat. J. Mech. Sci. 1972. Vol. 14, No. 3. P. 171-183.

Поступила в редакцию 31 августа 2009 г., с доработки — 21 мая 2010 г.