Конечные элементы повышенной точности для решения трехмерных задач теории упругости Текст научной статьи по специальности «Физика»

МАШИНОСТРОЕНИЕ

УДК 539.3

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

© 2003 г. П.П. Гайджуров

1. Введение. На современном этапе развития науки и техники резко возросла роль информационных технологий (САЭ/САМ/САЕ-программ) при проектировании новых машин, зданий и сооружений. Наиболее важное значение, определяющее прочность и несущую способность проектируемой конструкции, имеет комплекс программ конечно-элементного анализа (БЕА-пакет). Бесспорными лидерами среди коммерческих БЕА-пакетов являются системы А№У8, М8С^А8ТЯАМ и Рго/ЕМвВМЕЕК [1]. Не вдаваясь в подробности сервисных и вычислительных возможностей данных вычислительных комплексов, отметим, что при исследовании трехмерных объектов в них, как правило, используются четырехгранные (тетраэдаль-ные) конечные элементы (КЭ). Причем программа автоматической генерации сетки работает таким образом, чтобы все элементы были равновеликими. В качестве масштаба (меры) шага сетки принимается характерный минимальный размер тела, например толщина листа или диаметр наименьшего отверстия. Отсюда высокие порядок результирующей системы уравнений и вычислительные затраты. Подчеркнем, что задание мелкого шага при дискретизации пространственной конструкции - процедура оправданная, так как позволяет учесть смещения «как жесткого целого» при деформировании конечно-элементного ансамбля и, тем самым, обеспечить требуемую точность вычислений.

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

Разработке теории метода конечных элементов (МКЭ), позволяющей получать высокую точность на разреженных сетках, посвящены работы [2-4]. С помощью этой методики, получившей название «мо-ментная» схема МКЭ, или сокращенно МСКЭ, строятся матрицы жесткости объемных КЭ в форме шести- и пятигранников. Суть данного подхода сводится к отбрасыванию или минимизации некоторых членов разложения деформаций в степенной ряд, реагирую-

щих на «жесткие» смещения. Уравнения связи между деформациями и перемещениями устанавливаются в результате следующих формальных процедур (рассмотрим на примере полилинейного восьмиузлового КЭ).

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

ит = ат + я® * 1 + а{т} *2 + х3 + а^ х 1 х2 + к +

+ а(8) х 1 х 2 х 3 >

где а) - коэффициенты полинома; х 1, х 2, х 3 -

локальные координаты КЭ (с целью упрощения местные и глобальные оси совмещаем). Здесь и далее индексы, обозначенные латинскими буквами, если не оговорено специально, изменяются от 1 до 3.

2. Выражения для и т подставляются в формулы

для деформаций

1

£ ij = 2

д ui д u + ■

д x ,■ д x

3. Деформации разлагаются в ряд Маклорена.

4. Функция перемещений дополняется до поликвадратичного полинома (пробный закон аппроксимации)

ит = и т + А и m,

где А и т = Ь(9 х2 + Ь™ х22 + Ь™ х3 + Ь™ х2 х 2 + к + + Ь(2°) х 1 х 2х\; Ь) - дополнительные коэффициенты полинома.

5. Повторяется пункт 2 с и т .

6. Сравниваются выражения для е,, полученные

на основе принятого и и пробного ии законов

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

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

На базе МСКЭ разработана процедура построения матрицы жесткости [О] объемного криволинейного шестигранного КЭ. Результирующее выражение в тензорно-матричной форме, приведенное в работе [4], имеет вид

[ О ] = ПО"']] „х„ ; (1)

(3их3и) (3x3)

[ О"'] = [ А ]т [^ ]т [ ][ А ];

(3х3)

[Е'1к1 ] = }}}О'к1 { у (1])}{ у (k¡)}тл[gdxl йх2 дхъ,

-1

где [А ] - матрица, устанавливающая связь между коэффициентами аппроксимирующего полинома и узловыми перемещениями; [^ ] - матрица, связывающая ковариантные компоненты тензора деформаций и коэффициенты аппроксимирующего полинома; [С1 ■>к1 ] - контравариантные компоненты тензора упругих постоянных; {у (1;)} - вектор, составленный

из координатных функций с учетом принципа МСКЭ (индексы (1 ) указывают на зависимость размерности вектора от индексов /, ]); g - определитель метрического тензора (dv = dx1dx2dx3); п - число узлов

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

К сожалению, формулы для вычисления матриц [А ], [^ у ] и вектора {у (1;)} в цитируемых источниках не приводятся. Поэтому программная реализация выражения (1) проблематична.

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

2. Методика решения. В качестве базового рассмотрим шестигранный изопараметрический КЭ с полилинейным законом представления перемещений и т и геометрии 2 т в виде 8

и т = X и (т)Ф к (Х 1. Х 2> Х 3) . к=1

13

Ф к (Х 1. Х 2. Х 3) = 8 П (1 + РгкХг);

8 г=1

8

2т = X 2 (тк) Ф к (Х 1. Х 2. Х 3) . к=1

Здесь и (т) и 2 (т) - узловые перемещения и координаты КЭ в глобальном базисе; ф к - функции формы, построенные на базе одномерных интерполяционных полиномов Лагранжа; ргк = ±1 - естественные координаты узлов элемента (задаются в виде матрицы 3x8 по столбцам).

Ковариантные компоненты тензора деформаций в локальном базисе вычисляем по формуле [5]

£ j= 2 (z m Um i + Z m U,

(2)

' 1 у 2 т . У т . ^ т . гт .у/'

(суммирование по повторяющемуся индексу) где 2 т ■= д 2 т / д Ху - компоненты тензора преобразования координат; и т. 1 = д и т / д х 1.

Деформации элемента представим в виде отрезков ряда Маклорена:

£11 = £11 + е11.2Х 2+ е11.3Х 3+ е11.23Х 2Х 3 ;

е22 = е22 + е22.1Х 1+ е22.3Х 3+ е22.13Х 1Х 3 ;

Е33 = Е33 + Е33.1Х 1+ е33.2Х 2 + е33.12Х 1Х 2 ; (3)

е О — е О + е О ?Х ? ; ^^ — 1Х 1 ;

где £ и , £ 1},а

1; ар . а. в = 1. 2. 3 - «моменты деформаций» в точке ~ Л „ .

т а т 1х 1=х 2=х 3=0

Связь между деформациями {е} и узловыми перемещениями {и} элемента с учетом выражений (2) и (3) представим в матричной форме

{е}=ММ. (4)

где блочная матрица [D] = [ [D] 1 [D] 2 к [D] 8 ] ;

(6x24)

субматрица

[D] к = [{D(1)} k ^(2)} k ^(3)}] k , k = 1.2.....8 ;

(6x3)

вектор-столбец

Лк [2т.1 + (2т.12 + 2т.1 p2k) Х2 + (2т.13 + 2т.1 p3k) Х3 + (2т.123 + 2т.12 p3k + 2т.13 p2k + 2т.1 p2k p3k) Х2 Х3]

{D (m)} k

p2k Р m, 2 + (zm,12 + zm,2 pik) X1 + (zm,23 + zm,2 p^k) X3 + (zm,123 + zm,12 p3k + zm,23 p1k + zm,2 p1k p3k) X1 X3 ]

p3k [zm,3 + (zm,13 + zm,3 p1k) X1 + (zm,23 + zm,3 p2k) X2 + (zm,123 + zm,13 p2k + zm,23 p1k + zm,3 p1k p2k) X1 X2]

zm,1 p2k + zm,2 p1k + (zm,13 p2k + zm,1 p2k p3k + zm,23 p1k + + zm,2 p1k p3k) X3

zm,1 p3k + zm,3 Ak + (zm,12 p3k + zm,1 p2k p3k + zm,23 ftk + + zm,3 p1k p2k ) X2

zm,2 p3k + zm,3 p2k + (zm,12 p3k + zm,2 p1k p3k + zm,13 p2k + + zm,3 p1k p2k) X1

£ 23 = £23 + £23,1X 1

Здесь обозначено

Zm,j Э,

д 2 Z,

x1= x2 =x3 =0

m'aß дx „Эx

a ß

x1=x2 = x3 =0

а, ß = 1,2, 3,

д 2 z

m, 123"

д x 1 д x 2 д x

t= 0 (m ~m , а ß= 0 ,

3= 0.

Для пятигранного КЭ деформации вычисляем по формулам:

е11 = е11 + е11, 3 х3 ; е22 = е22 + е22, 3 х3 ; е33 = е33 + е33 , 1 х1 + е33 , 2 х2 ;

Блочная матрица [Б], входящая в уравнение (4), приобретает структуру типа

[Б] = [[Б] 1 [Б] 2 ...[В] 6 ] ,

(6x18)

где блок [Б] к = [{Б(1)} к{Б(2)} к {Б(3)}] к,к = 1,2,6 ;

(6x3)

субматрица

{(m) } = 8

Zm, 1 Ф k , 1+ (~m , 13 Ф k ,1+ ~m , 1 Ф k, 13) x3

~m , 2 Ф k , 2+ (~m, 23 Ф k, 2+ Zm, 2 Ф k, 23) x3

Zm , 3 Ф k, 3+ (Zm ,13 Ф k , 3+ Zm, 3 Ф k , 13) x1 + (zm, 23 Ф k , 3+ + Zm , 3 Ф k , 23) x2

Zm , 1 Ф k , 2+ Zm, 2 Ф k , 1+ (Zm, 13 Ф k , 2+ Zm , 1 Ф k , 23+ + Zm , 23 Ф k , 1+ Zm, 2 Ф k, 13) x3

Zm, 1 Ф k , 3+ Zm, 3 Ф k , 1

^, 2 ф к , 3+ Zm, 3 ф к , 2

Здесь производные от функций формы:

ф к, а = дф к / д ха I 0 ;

|х1 =х2 = х3 = 0

Ф k, aß =д 2 Ф k / д xaд x ß

x1 = x2 = x3 = 0

Матрицу жесткости [ к ] и вектор-столбец узловых сил {/} элемента формируем на основании принципа возможных перемещений

[ h ] = j[ D ] T [ E ][ D ] dv ;

(3nx3n) v

(5)

{/} = j ([F]T {q} + [D]T[E]{£0 } - [D]T{O0 }) dv +

(3n)

х1= х2 = х3 =0

Для случая, когда оси х т ортогональны, имеем

j [ F ]T { p } ds.

(6)

где блочная матрица [ Е ] = [[Е]1 [Е] 2 ... [Е] п ],

(3x3п)

[^]к = diag [фк фк фк ] (для пятигранного элемента

Фк к); № = {?(1) 2(2) ?(3)}Г , {Р} = {Р(1) Р(2) Р(3)}Т - соответственно векторы объемных и поверхностных сил, заданные в глобальном базисе (и р^ - физические компоненты). Элементы матрицы упругих постоянных [Е] выразим через физические компо-

тензора модулей упругости по формуле

ненты E,

(ijki)

Eijkl =

E(

(i j ki)

■^SiiS j j 8k k Sil

а, в = 1,2,3 . Выражения для функций ф к получаем

на основании формул для функций ф к с помощью

принципа вырождения (объединения соответствующих узлов).

Выражения для элементов матрицы [ Б] для шести- и пятигранных КЭ получены с помощью символьного процессора в среде Мар1е V Я5.

(по индексам в круглых скобках суммирование не производится); gi i, gjj, ... - ковариантные компоненты метрического тензора.

Вычисление интегралов, входящих в выражения (5) и (6), осуществляем численно с использованием квадратурных формул Гаусса.

3. Численное исследование сходимости. В настоящее время накоплен определенный опыт в области апробации конечных элементов. В частности, для статического анализа конструкций в форме пластин и оболочек разработчики БЕА-программ используют ряд тестовых («учебных») задач, решения которых получены аналитически или численно с применением более жестких гипотез, чем в тестируемом алгоритме МКЭ. Исследование сходимости разработанного конечно-элементного алгоритма выполнялось на следующих тестовых примерах.

Пример 1. Квадратная и круглая пластины с относительными размерами а/И=20 и а/И=100 (а - длина ребра / радиус пластины, И - толщина пластины) при жестком и шарнирном закреплении на контуре. Варианты нагрузки - равномерно распределенная и сосредоточенная в центре пластины сила. Варьировалось число элементов по толщине (однослойная и двухслойная модель) и в плане. Вычисления выполнялись при значениях коэффициента Пуассона V = 0 и V = 0,25. Расчетная схема строилась с учетом симметрии геометрии, граничных условий и схемы нагруже-ния. При ансамблировании круглой пластины в зоне

д

Z

m

a

m

£12 = £12 + £12, 3 x3 ; £13 = £13 ; £23 = £23 •

полюса использовались пятигранные КЭ. Результаты, полученные с помощью разработанного «моментно-го» алгоритма МКЭ, сравнивались с аналитическим решением [5] и данными численного расчета, выполненного с помощью коммерческой БЕЛ-программы СОSМОS/Dеsign 8ТЛЯ. Установлено, что разработанный комплекс позволяет достичь приемлемой точности по прогибу в центре пластины (5 тах < 2 %) при

использовании двухслойной схемы и равномерной разбивки пластины в плане на 16х16КЭ для а/И = 20 и 32x32 КЭ для а/И = 100. Количество элементов ансамбля в системе СО8МО8 при такой же точности соответственно составило 2177 и 6396. Было также выявлено, что стандартные изопараметрические элементы абсолютно неконкурентоспособны «момент-ным» элементам на принятых сетках.

Пример 2. Толстостенные квадратная и круглая плиты с относительным размером а/И = 5, жестко защемленные по контуру и нагруженные равномерно распределенной нагрузкой. Конечно-элементная модель, аппроксимирующая 1/4 часть плиты, представляла собой четырехслойную схему с равномерным шагом сетки в плане 10x10 КЭ. Численное решение сравнивалось с данными [6], полученными «методом определяющих состояний». В данном примере результаты разработанного алгоритма и стандартной схемы МКЭ дают практически одинаковую точность по перемещениям в центре (5 тах = 0,1 %). Наибольшая концентрация напряжений наблюдается на лицевых поверхностях в местах заделки и в центре плиты, во внутренних слоях напряжения минимальны.

Пример 3. Цилиндрический свод, нагруженный собственным весом, по торцам опирающийся на не-деформируемые диафрагмы и имеющий свободные продольные края. Этот пример использовался многими разработчиками БЕЛ-программ для апробации новых оболочечных КЭ на предмет учета «жестких» смещений. Исходные данные: средний радиус свода Я = 7,62 м; длина образующей Ь = 15,24 м; толщина И = 7,62-10-2 м; величина проекции внешнего давления на вертикальную ось q = 4,37-103Н/м2; модуль упругости Е = 2,1-1010 Н/м2; коэффициент Пуассона V = 0. По толщине свод моделировался однослойной конечно-элементной схемой. Использовались восьмиузловые объемные КЭ, построенные по «моментной» и стандартной схемам МКЭ. Установлено, что первый тип элемента достаточно точно моделирует поведение свода на сетке 32*32 КЭ. Стандартный (изопарамет-рический) КЭ на такой же сетке приводит к неверному результату, так как неточно описывает смещения как жесткого целого. На рис. 1 и 2 представлены результаты конечно-элементного моделирования для цилиндрического свода (перемещения увеличены в 100 раз).

Z3

141210-

Zi

Zi, м

а б

Рис. 1. Конечно-элементная модель свода: а - до деформации; б - после деформации

Н/м2

Z2 , м

Рис. 2. Картина распределения эквивалентных (по Мизесу) оэкв напряжений в цилиндрическом своде

Кроме рассмотренных примеров были решены также следующие задачи: плоский изгиб консольной балки; разрезное кольцо, нагруженное на свободном конце сосредоточенной силой; толстостенный цилиндр под действием внешнего давления (задача Ламе); цилиндр, нагруженный в центре двумя сосредоточенными диаметральными силами. Во всех тестах разработанный «моментный» элемент обеспечивает высокую точность и монотонную сходимость на относительно разреженных (в сравнении с пакетом COSMOS) сетках.

Литература

1. Brebbia C.A. Finite Element Systems. A Handbook. N. Y.,

1985.

2. Сахаров А. С. Моментная схема конечных элементов (МСКЭ) с учетом жестких смещений // Сопротивление материалов и теория сооружений. 1974. Вып. 24. С. 147— 156.

3. Завьялов Г.Г., Киричевский В.В., Сахаров А.С. Уточненные схемы МКЭ для расчета массивных конструкций // Проблемы прочности. 1978. № 6. С. 76-82.

4. Метод конечных элементов в механике твердых тел / Под ред. А.С. Сахарова и И. Альтенбаха. Киев; Лейпциг; 1982.

5. Седов Л.И. Механика сплошной среды. Т. 2. М., 1994.

6. Лисицин Б.М. Расчет защемленных плит в постановке пространственной теории упругости // Прикладная механика. 1970. Т. 6. Вып. 5. С. 18-23.

Южно-Российский государственный технический университет (НПИ)

17 сентября 2002 г.