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

УЧЕНЫЕ ЗАПИСКИ НАГИ

Том VIII 1 977 М3

УДК 518:512.831 539:377

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

А. Б. Кудряшов

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

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

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

1. При наличии начальных деформаций функционал потенциальной энергии для одного элемента имеет вид

П = -±-| а*(е - е0)гіа - | и./оГа - [и-рйЪ, (1)

2 и 2

с0 (&Охх’ “ОУУ’ £Огг' 2е0жу, 2в0хг, 2е0уг)

■—шестимерные векторы, составленные из компонентов тензоров напряжений, деформаций и начальных деформаций соответственно („*“ означает знак транспонирования)*,

/> Р> 11 ~ трехмерные векторы объемных и поверхностных сил и перемещений;

2 — объем, занимаемый конечным элементом;

£ — его поверхность.

Напряжения и деформации в (1) связаны между собой законом Гука

' а — О (в — е0), (2)

где О —матрица связи, размером 6X6.

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

а = 017; (3)

и = Ыи, (4)

где и — вектор обобщенных узловых перемещений; й и N—соот-

ветствующие матрицы связи.

Подстановка (2) — (4) в (1) дает после ряда преобразований

П = и* й* ОйийО, — | ^*0*0ео<ЙЗ +

а а

+ -1- [ ео Ое0 ая — $и* ЛГ* /<Ш — | и* ЛГ* р 41. (5)

2 9 2

Минимизация функционала по и приводит к системе линейных

алгебраических уравнений для конечного элемента

Ки = Р + Р0, (6)

где

К = | Я* вй (19 (7)

—матрица жесткости конечного элемента;

р=^ы*/ай+^р(И. (8)

2 а

— вектор обобщенных внешних узловых сил;

Р0 = $0* Ое0сг9 (9)

— вектор обобщенных внешних узловых сил, обусловленных начальной деформацией.

* В частном случае изотропной температурной деформации вектор начальной деформации имеет вид

£о = (аГ, аТ, аТ, 0, 0, 0),

где Т — температура в точке, а — коэффициент линейного температурного расширения.

Для изопараметрических конечных элементов в объем, занимаемый элементом, вводится криволинейная система координат тг), С, каждая из которых меняется на элементе от — 1 до 1. Переход к интегрированию в этих координатах дает

1 1 1

Г0= [ | ^ й* Ое0 \ de^.J \ йЫг\с1ы, (10)

-1 -1 -1

где

дх ду дг

Ж Ж Ж

дх ду дг дт\ <5г) дт\

дх ду дг

Ж Ж Ж

— матрица Якоби.

Интеграл (10) вычисляется, как правило, численно.

2. С помощью описанного подхода автором были разработаны алгоритмы формирования матриц жесткостей, вектора правой части и тензора напряжений для задач с температурными деформациями. Указанные алгоритмы были реализованы в комплексе программ „Система-4“, составленного автором совместно со Сниса-ренко Т. В., Чубанём В. Д. и Шевченко Ю. А.

При этом были рассмотрены следующие типы изопараметрических конечных элементов.

Трехмерные элементы сплошной среды [5] (фиг. 1 ,а). Каждый из них имеет форму криволинейного параллелепипеда с узлами в вершинах и с одинаковым числом (от 0 до 2) промежуточных узлов* на каждом ребре элемента. Закон аппроксимации геометрической формы элемента имеет вид

т

Я (5. С) —£/?,#,(?. 7), С), (11)

/~1

—*■

где /?($, С) —радиус-вектор точек %, % С элемента в общей системе координат (общей системой координат мы называем ту, в

которой задается геометрия элемента); Я1 — радиус-вектор г-го узла элемента**; N^1, -ц, С) — функция формы г-го узла элемента; /и—-число узлов элемента.

Аппроксимация поля перемещений принимается в виде

т

«(?, *1, С) = Х“/^(«. С), (12)

г=1

где и{%, -г), С)'—перемещение точек %, 7], С элемента; щ — перемещение г-го узла элемента.

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

** Выражения для трех-, двух- и одномерных функций формы приведены

в [Б].

Локально-двумерные элементы (фиг. 1,6). Каждый из этих элементов имеет форму криволинейного четырехугольника переменной толщины с узлами в вершинах и с одинаковым числом (от

О до 2) промежуточных узлов на каждой стороне элемента. Ап-

проксимация геометрии каждого из таких элементов осуществляется по формуле

т / \

яв ъ о=2 (я< + к‘) в ^ (13>

где к1 — толщина элемента в 1-м узле;

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

J, 7

а — нольточечный трехмерный элемент, б — двухточечный локально-двумерный элемент, в — одноточечный локально-одномерный элемент

Фиг. 1. Изопараметрические конечные элементы

Местная система координат в данном случае определяется соотношениями:

- а* / дЯ • > < | \ / •/ » / (< >, / УД / дЯ

1 — д5 / ; / = к! X *'; Л’ = V' х / <*т,

где = /? (£, Г|) задает положение срединной поверхности элемента.

Элемент пластины, работающей в своей плоскости, рассмотрен в [5]. В качестве неизвестных в узлах элемента принимаются два перемещения в плоскости пластины в общей системе координат. Аппроксимация поля перемещений производится по закону

т

“(Е> Ъ С) = 2и^(£. (15)

1=1

Элемент безмоментной оболочки. В качестве неизвестных в узлах элемента принимаются три перемещения в общей системе координат. Аппроксимация поля перемещений производится по закону (15).

Ш

Элемент моментной оболочки рассмотрен в работе [6]. В качестве неизвестных в узлах элемента принимаются три перемеще-

—» —

ния и два угла поворота относительно ортов V и )' в местной системе координат, определяемой соотношениями (14). Аппроксимация перемещений производится по закону

. т ~ ~

И (5, ч, С) = £

к = 1

ик + ('

}кУх'к + ІьУу'к)

\ .

Ц). (16)

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

Н, я с) = 2

і = і

4-

-г] А,-

і'і “Ь

С*/

(17)

где е1 — эксцентриситет в направлении /' в узле г;

/гг — толщина элемента в направлении у' в узле г;

Ьь — толщина элемента в направлении к' в узле г, и 1гг — орты местной системы координат в узле г, которые опре-

деляются следующим образом:

V =

с»5

у = к'Х і'; к' = к.

(18)

При этом предполагается, что £?=/?(£) задает положение оси

->

стержня и он расположен в плоскости, нормальной к вектору

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

Ъ С) = £и,.ЛЛ(!). (19)

Элемент балки, работающей на растяжение — сжатие и изгиб в своей плоскости. В качестве неизвестных в узлах элемента принимаются два перемещения и угол поворота относительно вектора № в местной системе координат, определяемой соотношениями (18). Аппроксимация перемещений производится по закону

и,

к= 1

Иь +

ГгЬ

(20)

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

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

где Т(%, т], С) — температура в точке £, -ц, С элемента; 7^ —значение ее в г'-м узле; Л^ —функция формы г'-го узла одно-, двух- или трехмерного элемента.

При реализации описанных выше элементов на ЭЦВМ зависимость характеристик среды от температуры учитывалась согласно соотношениям:

где Е — модуль упругости материала; а-—коэффициент линейного температурного расширения; V — коэффициент Пуассона.

3. В качестве примеров применения комплекса программ „Система-^ к расчетам конструкций с температурными деформациями рассмотрим следующие задачи.

Расчет температурных напряжений в толстой трубе конечной длины при осесимметричном нагреве. В [7] приведено аналитическое решение этой задачи для произвольного распределения температур. В расчете был принят закон изменения температуры 7’(г) = &г2. Для этого закона имеем

Были приняты следующие значения входящих в (23) и (24) параметров:

Расчет производился для кольцевого сектора трубы с углом 90° с постановкой соответствующих граничных условий. Использовались трехмерные элементы с двумя промежуточными узлами. Число элементов в расчете — 16, по 4 в радиальном и окружном на* правлениях. Общее число переменных составило 780. Получено практически точное совпадение численного и аналитического решений (фиг. 2).

Расчет температурных напряжений в конической оболочке. В [8] приведено аналитическое решение задачи термоупругости для конуса с параметрами

которое получено в предположениях гипотезы Кирхгофа — Лява.

Закон изменения температуры вдоль образующей конуса принимался 7’(5)=7"252. В окружном направлении температура не изменялась. Расчеты проводились при значениях параметров

/г=1 см; а = 10~6 1 /град; £=107 кГ/см2; Т^=\ град/см2.

т

(21)

£=£„ + £, Т+Е2 Т*; а = а0 + 04 Т а2 Т2\

4 = 7о + *1 Т + *2

(22)

(24)

(23)

а = 12-10 6 1/град; Е = 2-Ю4 кГ/мм2;

£ = 0,05; '^ = 0,3; а = 40 мм; й— 100 мм.

52/Л=10; = 0,4682 тс; 51 = 0,2 52,

8—Ученые записки № 3

11$

О — расчет методом конечного элемента,

----- аналитическое решение

Фиг. 2. Температурные напряжения в толстостенной трубе

На фиг. 3 и 4 приведены результаты расчетов по методу конечного элемента и аналитическое решение. Для расчета методом конечного элемента из конуса вырезался сектор с углом раствора 15° с постановкой соответствующих граничных условий. Было проведено два расчета с числом элементов вдоль образующей 6 и 12 (в окружном направлении в обоих расчетах брался один элемент) с общим числом переменных 260 и 500 соответственно. В расчетах использовались двухточечные элементы. Результаты обоих расчетов по методу конечного элемента практически совпали между собой.

Наблюдаемое различие в перемещениях, полученных аналитическим путем и по методу конечного элемента, следует отнести, по-видимому, к тому, что исследуемая оболочка является недостаточно тонкой (А/($2 — 5^= 1/8) и аналитическое решение, проведенное в рамках гипотезы Кирхгофа—Лява, является для данного конуса неточным*.

Температурные напряжения в цилиндрической оболочке, подкрепленной продольными ребрами. Рассмотрим замкнутую подкрепленную оболочку, состоящую из 2 т цилиндрических панелей и промежуточных продольных ребер. Оболочка имеет постоянную толщину /г, радиус /? и длину I. Е и V — модуль упругости и коэффициент Пуассона материала оболочки; ^0, Е0 и а0 — площадь поперечного сечения, модуль упругости и коэффициент линейного

* При построении элементов моментной оболочки, использовавшихся в этом расчете, принималась более общая гипотеза Рейсснера [9].

------ расчет методом конечного элемента,

О — аналитическое решение

Фиг. 3. Перемещения конической оболочки при неравномерном нагреве

о-

№

О

----- расчет методом конечного элемента,

О — аналитическое решение для внутренней поверхности,

X — аналитическое решение для наружной поверхности

Фиг. 4. Напряжения в конической оболочке .

при неравномерном нагреве

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

Результаты, нанесенные сплошными линиями, получены в [91 аналитическим путем при упрощающих предположениях (теория Лурье — Власова). Штриховыми линиями изображены результаты расчета по методу конечного элемента. Приведенные величины 5, ЛГ0 и л; связаны со сдвигающими усилиями в оболочке вдоль

К

0,6

0,4

0,2

О 0,1 0,2 0,3 0,4 %

------- метод конечного элемента,

-------в предположениях теории Лурье — Власова

Фиг. 5. Расчет подкрепленной цилиндрической оболочки

ребра 5, сжимающими усилиями в ребре Ы0 и координатой х соотношениями

5 =----------------------—— ; М0 =-^-; = —.

2а0Т0Е0Р(! 0 а0Т„Е0Р0 I

Решения проведены для числа /га = 1 для значений параметров ~гп = (IIЯ) |/=20 и % = ЕН1Ы Е0 4.

Имеет место хорошее согласование результатов на всех участках цилиндра, кроме его торцов, где теоретическое решение дает логарифмическую особенность для сдвигающих сил 5. Разумеется, расчет по методу конечного элемента не может дать „бесконечность" для 5 при л: = 0,5*, и вблизи этой точки на расстояниях

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

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

* Конечпоэлементным расчетом было получено значение 5 = 38,3 при *=0,5. _

** Теоретическое значение сжимающих усилий в оболочке М, = О при *=0,5.

Расчетная часть разбивалась на 135 одноточечных элементов мо-ментной оболочки (15 — вдоль ребра и 9 в окружном направлении). В окружном направлении разбивка на элементы была более густой у ребер, чем на периферии. Продольное ребро было разбито на 15 стержневых одноточечных элементов. Расчет производился при значениях параметров

= 25 см; Л-—1 см; / = 100 см; Е = 6283, 1852 кГ/см2,

V = 0,3; Г0 = 5 см2; Е0 = 104 кГ/см2; ос0 = 10~6 1/град; Г=100°.

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

ЛИТЕРАТУРА

1. Аргирис Д ж. Энергетические теоремы и расчет конструкций. Сб. .Современные методы расчета сложных статически неопределимых систем* под ред. Филина А. П., Л., Судпромгиз, 1961.

2. Галлахер Р., Падлог П., БейлондрДж. Анализ напряжений в конструкциях сложной формы, подверженных нагреву. „Ракетная техника и космонавтика” (русский пер.), № 5, 1962.

3. Иванов Ю. И. К расчету температурных деформаций и напряжений методом конечного элемента. Труды ЦАГИ, вып. 1728, 1975.

4. Т а b а г г о к В., Н о а V. S. Thermal stress analysis of plates and shells by hybrid finite element method. „CANCAM 73, Montreal 1973“, 1973.

5. Зенкевич О. Метод конечных элементов в технике. М., .Мир*, 1975.

6. Ahmad S., Irons В. М., Z i е n k i е vv i с h О. М. Analysis of thick and thin shell structures by curved finite elements. International J. Num. Meth. in Engng., 1970, N 4.

7. Боли Б., Уэйнер Дж. Теория температурных напряжений. М., „Мир", 1964.

8. Коваленко А. Д. Основы термоупругости. Киев, „Наукова думка“, 1970,

9. Б а л а б у х Л. И., Шаповалов Л. А. Исследование температурных напряжений в цилиндрической оболочке, подкрепленной продольными ребрами. В сб. „Расчеты на прочность", вып. 12, М., „Машиностроение", 1966.

Рукопись поступила 20jIV 1976 г.