Научная статья на тему 'Конечные элементы повышенной точности для решения трехмерных задач теории упругости'

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

CC BY
149
29
i Надоели баннеры? Вы всегда можете отключить рекламу.

Аннотация научной статьи по физике, автор научной работы — Гайджуров П. П.

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

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по физике , автор научной работы — Гайджуров П. П.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

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

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

УДК 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

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

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 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.