Решение двумерных задач термоупругости и термоупруго-пластичности методом локальных вариаций Текст научной статьи по специальности «Математика»

МЕХАНИКА

J

УДК 539.434+681.3.06

И. В. Станкевич, Ле Шонг Тунг

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

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

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

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

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

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

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

Математическая формулировка задачи. Рассмотрим в двумерной области Г (с границей дГ = С и С2) квазистатическую задачу термоупругости для обобщенного деформированного состояния при некотором заданном распределении температурной деформации еТ, г,] = 1, 2, и фиксированных кинематических и силовых условиях нагру ения.

Пусть оси х1 и х2 декартовой системы координат лежат в плоскости поперечного сечения тела. Для линейно-упругого изотропного материала с учетом температурной деформации напря ения определя тся по

формулам [2, 3, 6]

аи = 1 [(1 " + ^е22 " (1 + ^^

^22 = 1 ^ [(1 - ^)е22 + ^е11 - (1 + (1) т12 =

где ^ — коэффициент Пуассона, О — модуль сдвига материала. Напряжения должны удовлетворять уравнениям равновесия

+ ¡1 = 0,

(2)

+ ¡2 = 0,

дап дт12

+

дх1 дх2

да22 дт12

+

дх2 дх1

здесь ¡1, ¡2 — проекции объемной силы на оси х1, х2 соответственно. Используя соотношения Коши

е = д^1 е = ди е = 1 / д«1 + дмЛ (3)

дх/ дх2' 2 \ дх2 дх1 у

и учитывая формулы (1), из уравнений равновесия получим

д „г п д

dxj J дхг

IJG

uj,j 1 - 2^

=

г дХг

G(1 +

1 - 2^

, i,j = 1, 2.

(4)

Здесь и далее запятой перед нижним индексом обозначена производная по соответствующей координате, а по повторяющимся индексам г, ] проводится суммирование от 1 до 2.

редполо им, что на участке границы С1 задана с компонентами рг, г = 1, 2, распределенная нагрузка

ик,к - (1+ ^ р . . 2^-1 - 2^-П + (иг,з + и,,г)п = о, г, к = 1, 2, (5)

где пг и щ — компоненты внешней нормали, а на участке границы С2 области ^ заданы перемещения (кинематические граничные условия)

и = и0, г = 1, 2. (6)

Возможности аналитического решения краевой квазистатической задачи термоупругости (4)-(6), (1) и (3) для двумерной области сложной формы при произвольном температурном поле и зависимости физико-механических свойств от температуры и координат ограничены. Поэтому приходится использовать численные метода решения, например, основанные на вариационной постановке.

Решение упругой задачи. Функционал для эквивалентной вариационной формулировки данной задачи применительно к линейно-упругому изотропному материалу имеет вид [1,2]

[и] = !К,к - 3ет]2 + 0[и2к к - ик,кщ,г] +

Р С 1

+ "^[и1,2 + и2,1 ]2|— £ргщйС, г, к = 1, 2. (7)

С

Функционал (7) на допустимых перемещениях имеет минимальное значение. Для минимизации функционала (7) в работе [2] был использован метод локальных вариаций [5, 7].

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

Для этого рассмотрим произвольный конечный элемент (рис. 1) сетки, построенной в области ^ .О бозначим перемещение по оси х1 через и = и1, а перемещение по оси х2 через V = и2. Кроме того, в каждом элементе обозначим узловые перемещения через иг, V, и, V], ик, vk в узлах г, ], к соответственно, а пары координат узлов через (хг, уг), (х], У]X (хк, Ук).

Перемещения в направлении осей х1 и х2 в каждом элементе определяются по приближенным выражениям

и = Ки + N и] + N ик, (8)

V = NгVг + К] V] + Кк Vk, (9)

где Кг, К], Кк — функции формы конечного элемента [4].

Формула для определения компонентов деформации в матричной форме имеет вид

[4]

(£11 £22 £12)т = В(и V и, V, щ ьк)т; (10)

здесь В — матрица градиентов, полученная соответствующим дифференцировани-а] J ем матрицы функций формы N [4].

Таким образом, функционал (7) можно Рис. 1. Перемещения узлов г

г приближенно заменить суммой:

конечноэлементнои модели ^

N

J = £ J(•>,

(11)

e=1

где 7(е) — функционал, определенный на элементе с номером е.

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

перемещении u

(0) v(0)

во всех узлах сетки (но с учетом кинематических граничных условий) и некоторые начальные значения шагов варьирования узловых перемещений Ни > 0, Н,0 > 0. С помощью формул (8)-(10) вычисляется значение функционала (11).

В каждом элементе варьируется сначала значение перемещения и в узле I (рис. 1) с шагом Ни. Из трех возможных значений щ, и + Ни и щ — Ни выбирается то значение, при котором функционал принимает минимальное значение. Подобные вариации проводятся до тех пор, пока при замене и на щ + Ни (или щ — Ни) не произойдет уменьшение значения функционала (11). После этого аналогично варьируется значение перемещения и, в ]-м и затем ик в к-м узлах. Аналогично варьируются значения перемещений V, V,, ьк с шагом Н,0. Процесс варьирования на данном этапе заканчивается тогда, когда разность численных значений функционала (11), полученных в двух последовательных приближениях, станет меньше заданного положительного достаточно малого числа т.е. будет выполняться условие

|j(n) - J(n-1)| <

(12)

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

шах(Ни, ) < Н*, (13)

где h* — заданное достаточно малое число.

Как отмечалось ранее, компонента вектора деформации определяются по формулам (10), а компоненты вектора напряжений вычисляются с помощью равенства

= D

£22 | - £ £12

T

(14)

где О — матрица упругих характеристик, ет — вектор температурной деформации. Матрица О и вектор ет для плоского напряженного и плоского деформированного состояний определяются по следующим формулам:

— для плоского напря енного состояния

D =

E

1 -

1 0

1 0

0 0 1 -2 )

£T = aAT I 1

(15)

— для плоского деформированного состояния

1

' 1 1

E(1 -

D =

1

(1+ ^)(1 - 2^)

1 - ^ 0

1 - ^ 1

0

0 0

1 - 2^ 2(1 - /

(16)

1

£T = (1 + p,)aAT I 1

0

Здесь а — коэффициент линейного расширения материала, зависящий в общем случае от температуры и координат; ДТ — разность средней температуры элемента и некоторой его начальной температуры.

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

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

метод локальных вариации в сочетании с методом переменных параметров упругости [3].

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

£11 = Е*[ап + Н*а22^

^22 = Е* [^22 + (17)

_ Т12

^12 =

2G*'

где E*, G*, — переменные параметры упругости, определяемые по формулам

аи 11 — 2^ аи

Е*___ н*__2_3Е £и _ аи (18)

= 1 + 1 - ои , Н _ 1 - Ои , _ 3£и' ^ '

3Е £и 3Е £и

Здесь £и, аи — интенсивности деформации и напряжений соответственно.

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

В исходном, нулевом приближении считается, что параметры Е*, Н*, С* равны параметрам упругости Е, н, С. Затем решается задача термоупругости, т.е. минимизируется функционал (7) методом локальных вариаций по приведенному ранее алгоритму. В результате определя т-ся значения перемещений. По этим значениям с помощью формул (10) и (14) вычисляются компоненты векторов деформации {£ ц } и напряжений {ац} и соответствующие им значения интенсивности £и и аи [3].

На рис. 2 в координатах £и, аи напряженно-деформированное состояние некоторого конечного элемента отмечено точкой 0, лежащей на луче, тангенс угла наклона которого пропорционален величине 3С [3].

По заданной диаграмме деформирования аи(£и) определяется интенсивность напряжений а*0. Это значение соответствует точке А пересечения диаграммы аи(£и) с вертикальной прямой, проходящей через точку £и0 _ £и (см. рис. 2). С помощью найденных значений и

еи0 по формулам (18) определяются новые значения параметров Е*, и С*, которые будут различными в различных элементах. Таким образом, возникает задача определения напряжений в неоднородном теле. Далее вновь методом локальных вариаций решается задача термоупругости, но уже с новыми параметрами. В результате решения находятся перемещения для первого приближения. Процесс коррекции параметров упругости продолжается до тех пор, пока результаты двух последовательных приближений не станут достаточно близки.

Результаты численного решения. При численном решении задачи упругопластичности в качестве объекта рассматривалась цилиндрическая труба с радиусами Ях = 150 мм и Я2 = 156 мм. Начальные значения узловых перемещений соответственно в радиальном и окрук-ном направлениях были приняты равными и(0) = 0,1 мм, г>(0) = 0 мм, а начальные значения шагов варьирования соответственно в радиальном и окружном направлениях были приняты равными ки = 0,001 мм, к, = 0,0001 мм. Процесс варьирования завершался, когда шаги варьи-

Рис. 2. Диаграмма деформирования

мм.

рования становились меньше к* = 10-10

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

Численные исследования проводились для трех вариантов нагру-ения.

В первом случае расчет проводился для трубы, находящейся под воздействием внутреннего давления, которое изменялось в диапазоне рг = 2,1... 2,5 МПа. Распределения радиальных и окружных напряжений по толщине трубы представлены на рис. 4. Из анализа результатов следует, что при увели-

г ' г Рис. 3. Конечноэлементная модель ци-

чении значения внутреннего дав- линдрической трубы

150 151 152 153 154 155 г, мм

-2,5

Ьг,мпа

150 151 15Z 153 154 155 г, ММ Ö

Рис. 4. Напряжения при действии внутреннего давления:

а — радиальное напряжение; б — окружное напряжение; 1 — р = 2,5 МПа; 2 р = 2,4 МПа; 3 — р = 2,3МПа; 4 — р = 2,2 МПа; 5 — р = 2,1 МПа

152 153 154 г, мм

Рис. 5. Напряжения при температурном воздействии:

а — радиальное напряжение; б — окружное напряжение

ления граница зоны пластичности г(р) перемещается от внутреннего радиуса г(р) = при рг & 2,1 МПа к внешнему радиусу г(р) = Л2 при рг & 2,5 МПа (рис. 4, б).

Во втором случае рассматривалось напряженно-деформированное состояние трубы под влиянием температурного поля. При этом температура распределялась по радиусу трубы по логарифмическому закону:

T (r) = Ть

1п(Д2/г) 1n(R2/Ri)

где Ть = 150 °С — температура на внутренней поверхности цилиндра. Графики изменения радиальных и окружных напряжений в зависимости от радиуса показаны на рис. 5.

В третьем случае труба была нагружена внутреннем давлением (рг = 2 МПа) и температурным полем (Ть = 150 °С). Результаты решения показаны на рис. 6.

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

о

-0,5 -1 -1,5 -2 -2,5 -J

bp,

Рис. 6. Hапряжения при термосиловом воздействии:

а — радиальное напряжение; б — окружное напряжение

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

СПИСОК ЛИТЕРАТУРЫ

1. Балабух Л. И., Колесников К. С., Зарубин В. С. Основы строительной механики ракет. - М.: Высшая школа, 1969. - 494 с.

2. Зарубин В. С. Прикладные задачи термопрочности элементов конструкции. - М.: Машиностроение, 1985. -290 с.

3. Малинин Н. Н. Прикладная теория пластичности и ползучести. - М.: Машиностроение, 1975. - 395 с.

4. Сегерлинд Л. Применение метода конечных элементов. - М.: Мир, 1979. -391 с.

5. Черноусько Ф. Л., Баничук Н. В. Вариационные задачи механики и управления. - М.: Наука, 1973. - 237 с.

6. Тимошенко С. П. Теория упругости. - М.: Наука, 1979. - 560 с.

7. Станкевич И. В., Ле Шонг Тунг. Определение напряженно-деформированного состояния конструкции методом локальных вариаций // Научные материалы Первой межд. конф. "Аэрокосмические технологии", посвященной 90-летию со дня рождения акад. В.Н. Челомея. Москва-Реутов, 24-25 мая 2004. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. - C. 125-129.

Статья поступила в редакцию 23.12.2004

Игорь Васильевич Станкевич родился в 1949 г., окончил в 1976 г. МВТУ им. Н.Э. Баумана и в 1988 г. МГУ им. М.В. Ломоносова. Д-р техн. наук, профессор кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Автор более 50 научных работ в области прикладной термомеханики.

I.V. Stankevich (b. 1949) graduated from the Bauman Moscow Higher Technical School in 1976 and Moscow State University n. a. M.V. Lomonosov in 1988. D. Sc. (Eng.), professor of "Applied Mathematics" department of the Bauman Moscow State Technical University. Author of over 50 publications in the field of applied thermal mechanics.