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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м V 1974

№ 6

УДК 512.831

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

А. Б. Кудряшов, В. Д. Чубань, Ю. А. Шевченко

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

В работе [1] были предложены новые принципы построения матриц жесткостей для конечных элементов (получивших название изопараметрических), что дало мощный толчок применению метода конечного элемента к решению трехмерных задач математической теории упругости [2], теории толстых оболочек [3], линейной теории теплопроводности [4] и многих других.

В согласии с этими принципами декартовы координаты произвольной точки объема 2, занимаемого конечным элементом, могут быть представлены в виде:

Я = £/?гЛЛ($, 7,, С), (1)

1=1

где /? = /?(£, у\, С) — радиус-вектор произвольной точки объема й;

Н1 — радиус-вектор г'-го узла элемента (г = 1, 2и — число узлов элемента); -Л/г(£, т), С) — интерполяционный полином /-го узла (1=1, 2,

Выражение (1) фактически означает введение криволинейной системы координат г\, С в области 2. При этом в качестве обла-

* Вид интерполяционных полиномов приводится в приложении.

IЧ. I п I. I ^ К1 •

Заметим, что выражение (1) можно рассматривать также как отображение куба со на криволинейную область 2.

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

Фиг. 1. Разбиение криволинейного параллелепипеда на элементы и типы агрегатов, применяемых в расчетах

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

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

Применяя несложные алгоритмы, с помощью выражения (1) можно получить радиус-векторы любого узла любого элемента для данного агрегата.

В зависимости от числа базисных точек авторами рассматривались три типа агрегатов, показанных на фиг. 1, б, в, г соответственно.

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

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

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

---------по формуле Журавского при Л/6 = 1;

---------по формуле Журавского при К/Ь = 3:

О — результаты расчета

Фиг. 3. Распределение касательных напряжений, возникающих в поперечном сечении консольно заделанной балки при ее косом изгибе

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

Внешний вид конструкции и разбиение ее на элементы представлены на фиг. 2, а.

На фиг. 3 представлены характер нагружения балки и результаты расчетов. Для сравнения приведены кривые распределения касательных усилий, полученных по формуле Журавского.

5— Ученые записки ЦАГИ № 6

65

Несовпадение результатов расчета и теоретической кривой при •й/й = 1 говорит о том, что формула Журавского (верная, вообще говоря, при ЩЬ~^> 1, [5]) плохо „работает" в данном случае.

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

узлами на каждом ребре. Учи' ' тывая осевую симметрию за-

дачи, решение проводилось для четверти цилиндра при соответствующих граничных условиях. ^ г & т, град

\Р= 6,6'6к г/ммг-/; ^50мм, г = 38 мм-, V = 6,3 ■,

£=7-10*яг/мм2

результа ты расчета

076

08В

20

10

і * /, ^б^-іа^т, А = 2-103£т/с< г,-5,0 см г = 3,8см у 'см2-и-грао

9гі г о результат і расчета.

г/г, 10

0,88

г/Гі 1,0

Фиг. 4. Распределение радиальных оЛ и тангенциальных и* компонент тензора напряжений по радиусу в бесконечной цилиндрической оболочке под внутренним давлением

Фиг. 5. Стационарное распределение поля температур по радиусу в бесконечной цилиндрической оболочке, находящейся под воздействием внешнего теплового потока

Для задания геометрии агрегата и его разбиения на элементы (фиг. 2, б) применялся агрегат с 32 базисными точками (см. фиг. 1, г).

На фиг. 4 представлены результаты расчета и для сравнения приведены теоретические кривые распределения радиальных и тангенциальных компонентов тензора напряжений.

Аналогично была решена задача о нагреве бесконечной цилиндрической оболочки постоянным внешним тепловым потоком для стационарного случая. На внутренней поверхности оболочки поддерживалась постоянная температура £ = 0.

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

Задача определения концентрации напряжений в цилиндрической оболочке, подкрепленной ребром жесткости. На фиг. 2, в изображены внешний вид конструкции и ее разбиение на элементы с одним промежуточным узлом на каждом ребре *.

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

* Расчеты производились при числе разбиений конструкции вдоль оси х вдвое большем, чем на фиг. 2, в.

всех узлах ребра жесткости, лежащих на его торце, были зафиксированы перемещения их. Все узлы на гранях оболочки, перпендикулярные к оси х, были закреплены, причем узлы, лежащие на грани, не примыкающей к ребру жесткости, получали одинаковые перемещения их— 10~3, т. е. нагружение конструкции проводилось кинематически.

----------- вдоль линии 1 по оболочке;

--------— вдоль линии / по ребру жесткости;

-----------вдоль линии 2

■Фиг. 6. Распределение компонент ах и ог тензора напряжений по оси х

----------вдоль линии 1\

—---------вдоль линии 2

Фиг. 7. Распределение компоненты ах тензора напряжений вдоль оси х в гладкой оболочке

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

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

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

Все расчеты с помощью метода конечного элемента, а также выполнение чертежей конструкции графопостроителем производились с помощью комплекса программ для ЭЦВМ БЭСМ-6 „СИСТЕ-МА-4“, составленного А. Б. Кудряшовым, Т. В. Снисаренко,.

В. Д. Чубань и Ю. А. Шевченко.

ПРИЛОЖЕНИЕ

Интерполяционные полиномы.

Пусть ?)г, С/ — криволинейные координаты г-го узла элемента.

Введем обозначения:

^0 = %, ^0 = ЧП-И Со = СС,-.

Тогда интерполяционные полиномы Л/^(£, •/), С) можно представить в виде:

а) для элемента с 8 узлами (см. фиг. 1, б)

ад ч, с)= 4-0 +^о)(1-Но)(1+С0);

б) для элемента с 20 узлами (см. фиг. 1, в)

N\ (£, "Ч, С) — -у (1 + Ео) 0 + т1о)(1 + С0)(Е0 + + Со — 2)

для узлов, лежащих в вершинах куба;

ад -П, 0=4- (1+&о)(1 + 7|0)(1-С2)

для узлов, лежащих в плоскости С = 0;

ад % <э=4-(1 +чо)(1 +с0)(1-?2)

для узлов, лежащих в плоскости 1 = 0;

ад ъ с)=4(1 + ^)(1+с0)(1-^)

для узлов, лежащих в плоскости 7] = 0;

68

в) для элемента с 32 узлами (см. фиг. 1, г)

ад Ч, С)~-^-(1 +Е0)(1 +Чо)(1 + бо) [—19+9(5* + ча + С*)]

для узлов, лежащих в вершинах куба;

А/ (5, Ч, С) = ^ (1 + 10) (1 - Чо) (1 - С2) (1 + 9 С0)

г I 1

для узлов, лежащих в плоскостях С=+-д-;

ад Ч, С) = -^-(1+Чо)(1+Со)(1-52)(1+9«о)

для узлов, лежащих в плоскостях Е=+-д-;

ад Ч, С) = -^(1 +У(1 +С0)(1-Ч2)(1 + 9^0)

, 1

для узлов, лежащих в плоскостях '1,1 = + -з--

Приведенные три семейства полиномов обладают тем важным свойством, что

ад. V Су) = а„,

где 8^ — символ Кронекера.

ЛИТЕРАТУРА

1. Irons В. М. Engineering applications of numerical integration in stiffness methods. AIAA J., 1966, vol. 4, No 11.

2. ErgatoudisJ., Irons В. М., ZienkiewiczO. M. Curved isoparametric quadrilateral elements for finite element analysis, Int. J. Solids Structures, 1968, No 4.

3. Ahmad S., Irons В. М., Zi.enkiewicz О. C. Analysis of thick and thin shell structures by curved finite elements. International J. Num. Meth. in Engng., 1970, No 3.

4. ZienkiewiczO. C., Parekh C. J. Transient field problems: two — dimentional and three — dimentional analysis by isoparametrical finite elements. International J. Num. Meth. in Engng., 1970, № 2.

5. Беляев H. М. Сопротивление материалов. М., Физматгиз,

1959.

6. Ч у б а н ь В. Д. Об одном экономичном алгоритме решения систем алгебраических уравнений метода конечного элемента. Журн. вычисл. матем. и матем. физ., 13, 3, 1973.

Рукопись поступила 8jX 1973 г. Переработанный вариант поступил 19\Н 1974 г.