Численное решение задачи термопластичности с дополнительными параметрами состояния

( Numerical solution of thermal plasticity problem with additional parameters
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Галанин М.П., Гузев М.А., Низкая Т.В.
(M.P.Galanin, M.A.Guzev, T.V.Nizkaya)

ИПМ им. М.В.Келдыша РАН

Москва, 2007
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 05-01-00618)

Аннотация

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

Abstract

This work describes the numerical algorithm for solution of thermoplasticity equations. As a model problem we considered heating and cooling of a plain plate, made of an elastoplastic material with temperature-dependant properties. We assume that the deformations and the deformation velocities are small and the plate is in plain stress state. The algorithm was programmed, some computational examples are provided.

СОДЕРЖАНИЕ

Введение                                                                                                                              3

§ 1. Математическая модель пластического течения с учетом температуры       4

§ 2. Вычислительный алгоритм                                            8

§ 3. Программная реализация и проведение тестовых расчетов                                17

Заключение                                                                                                                         18

Список литературы                                                                                                          20


Введение

 

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

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

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

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

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

Авторы выражают благодарность И.В. Станкевичу и Ю.М. Темису за полезные обсуждения и советы.

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (проект РФФИ № 05 - 01 - 00618).

 

1. Математическая модель пластического течения с учетом температуры

 

1.1. Уравнения движения

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

Тогда уравнения движения среды принимают вид [6]:

                                      (2.1)

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

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

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

При малых перемещениях верно следующее соотношение между тензором полной деформации и перемещениями:

.

Связь между деформациями и напряжениями задается законом Гука:

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

.

Температурная деформация подчиняется закону теплового расширения:

Здесь T0 - температура естественного состояния среды.

Закон изменения  можно записать в следующем виде:

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

Тогда скорость изменения полной деформации можно записать в виде:

                (2.2)

В итоге система уравнений принимает вид:

                                             (2.3)

В дальнейшем мы будем рассматривать изотропную среду, для которой закон Гука имеет вид:

где – модуль сдвига,  - модуль объемного сжатия, ,  - девиаторы  деформаций и напряжений соответственно.

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

 

1.2 Теория пластического течения

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

В теории пластического течения [2] для описания поведения материала вводят функцию упрочнения  

Основные предположения теории пластического течения таковы:

1.     Напряжения в любой точке удовлетворяют условию

2.        При  материал ведет себя упруго, т.е. .

3.        Выполнен принцип максимальной мощности: .

4.        Изменения объема всегда происходят упруго, т.е. .

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

                             (2.4)

где  имеет смысл множителя Лагранжа, а  означает производную от по времени при условии . Условие  называют условием активного нагружения.

Если материал изначально изотропен, то тензорные величины могут входить в функцию  только в виде инвариантов. Изменение объема всегда упруго, поэтому функция упрочнения не зависит от первого инварианта тензора напряжений. Более того, для широкого класса материалов можно считать, что от третьих (кубических) инвариантов напряжений и деформаций эта функция также не зависит.

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

Выпишем явный вид функции упрочнения для некоторых моделей.

 

Идеальный упруго-пластический материал (материал без упрочнения) [2,6]:

где  - предел текучести, зависящий только от температуры.

Материал с изотропным упрочнением [6]:

где  - предел текучести,  - параметр Одквиста.

Материал с трансляционным упрочнением [6]:

где  - тензор микронапряжений.

Модель Кадашевича - Новожилова [7]:

.

Модель с учетом несовместности пластических деформаций [8]:

 

2. Вычислительный алгоритм

 

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

Введем разбиение  и усредним уравнения движения по временному интервалу :

Интегралы будем аппроксимировать следующим образом:

Считая известными все значения на предыдущем временном слое, получим уравнения относительно приращений , :

.

Записав уравнение для невязки   можно указать диапазон значений , при которых она не возрастает:

Граничные условия записываются аналогичным образом:

Связь между напряжениями и деформациями задается соотношениями:

где

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

В классических моделях пластичности удается получить явное выражение для параметра  в зависимости от скорости изменения напряжений и температуры. Тогда при выполнении условий активного нагружения можно выписать нелинейные уравнения относительно перемещений и решать их методом простой итерации (он же метод переменных параметров упругости) или методом Ньютона [4,5]. Алгоритм решения такой задачи методом переменных параметров упругости приведен в 2.1.

В модели с учетом несовместности пластических деформаций [8] в закон изменения дополнительного параметра  входят производные по пространству. Выписать явное выражение для  в этом случае невозможно, для его вычисления придется использовать метод типа проекции на ограничение. Смысл его заключается в следующем: вначале решают задачу без учета ограничений, а затем методом последовательных приближений определяют приращение пластической деформации, обеспечивающее выполнение ограничения. При этом учитывается смещение поверхности нагружения за счет зависимости дополнительных параметров от поля пластических деформаций. Алгоритм метода дополнительных деформаций с проекцией на ограничение приведен в 2.2.

 

2.1. Метод переменных параметров упругости

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

Учитывая, что , а , получим следующее выражение для параметра :

.                                    (2.5)

Тогда скорость пластической деформации при активном нагружении имеет вид:

                                    (2.6)

где  - тензор активных напряжений, а  - параметр упрочнения.

С учетом (2.4) равенство (2.2) можно записать в виде:

После преобразований получим:

Отсюда:

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

Тогда в предположении о монотонности функций  получим

,

где

, ,  -

значения в некоторой промежуточной точке. Здесь и далее индекс k  при  опущен для краткости.

Вычислим приращение пластической деформации за шаг по времени, считая, что тензорные функции , монотонны по  при :

,

где , - значения в некоторой промежуточной точке. Теперь проинтегрируем соотношения (2.2) по t , зафиксировав некоторую среднюю поверхность пластичности, соответствующую значению : 

Отсюда

                                  (2.7)

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

Таким образом, на каждом временном шаге имеем систему уравнений вида:

дополненную граничными условиями:

Тензоры  определены выше. Их средние значения вычисляются по формулам (2.7).

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

На каждой итерации s=0,1,2... метода переменных параметров упругости решаем уравнение следующего вида:

с соответствующими граничными условями.

В качестве нулевого приближения задается  .

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

Критериями прекращения  итераций выбраны:

по  - ,

по s.

Таким образом, алгоритм решения задачи на очередном временном слое выглядит следующим образом:

0.     Начало временного слоя:,

(равенство понимается здесь в смысле присваивания)

1.     Вычисление невязок предыдущего слоя:

,  .

2.     Вычисление правой части:.

3.     Решение задачи в предположении об упругости материала:

где

4.     Вычисление нулевого приближения напряжений: .

5.     Нахождение областей активного нагружения, в которых

выполняются условия  и .

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

7.     Расчет очередной итерации метода переменных параметров упругости:

                 где

8.     Вычисление  и  по текущим значениям ,,.

9.     Пересчет ,, с учетом новых значений  и .

10.  Проверка условия сходимости итераций по p:

если выполнено, переход на 11, если не выполнено – на 8.

11. Проверка условия сходимости итераций по s: 

если выполнено, переход на 0, если не выполнено – на 6.

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

Конечномерное вариационное уравнение на каждой из итераций выглядит следующим образом:

где  - координаты -го узла сетки (общее число узлов – N),  приращение перемещений в -м узле,  - кусочно-линейные базисные функции, удовлетворяющие условиям .

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

 

2.2. Метод дополнительных деформаций

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

Решение на каждой итерации этого метода находят в два этапа [5]. На первом этапе рассматривают упругую (без учета ограничений) задачу и находят нулевое приближение деформаций и напряжений:

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

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

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

 

1.     Начало временного слоя:,

(равенство понимается здесь в смысле присваивания)

2.     Вычисление невязок предыдущего слоя:

,  .

3.     Вычисление правой части:.

4.     Решение задачи в предположении об упругости материала:

    где

5.     Вычисление предиктора: .

6.     Нахождение областей активного нагружения, в которых выполняются условия  и .

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

8.     Пересчет .

9.     Вычисление очередной итерации :

 где

10. Вычисление напряжений:

11. Вычисление очередного приближения пластической деформации из условия , .

12. Проверка условия :

если выполнено, переход на 0, если не выполнено – на 8.

 

3. Программная реализация и проведение тестовых расчетов

 

Алгоритм метода переменных параметров упругости реализован с использованием программного комплекса, разработанного в [9]. Проведены тестовые расчеты для модельного материала, обладающего изотропным упрочнением:

,

Зависимость коэффициентов упругости от температуры взята следующей:

Температурная деформация является изотропной:

Температурное поле задавалось следующим образом:

,

где .

Зависимость температуры пятна от времени:

Траектория движения:

Значения параметров: , , , .

Граничные условия:   

На рисунках приведено температурное поле, поле эффективных напряжений  и поле эффективной пластической деформации в различные моменты времени: t=0.06, t=0.2, t=0.4, t=0.6, t=1 (сверху вниз).

Видно, что после окончания процесса в материале остаются ненулевые напряжения.

 

Заключение

 

Рассмотрена задача о локальном нагреве и охлаждении тонкой пластины из упругопластического материала, свойства которого зависят от температуры. Для описания пластических свойств материала использована модель пластического течения. В работе приведены различные виды поверхностей нагружения и алгоритмы численного решения уравнений модели на основе метода переменных параметров упругости [4] и метода дополнительных деформаций [5]. Метод на основе метода переменных параметров упругости реализован программно с использованием конечных элементов 1-го порядка.

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

 

 

Рис.1. Температурное поле, поле эффективных напряжений и эффективных пластических деформаций (слева направо) в различные моменты времени (сверху вниз: t=0.06, 0.1, 0.2, 0.4, 0.6, 1)

 

4. Список литературы       

1.     Биргер И.А. Остаточные напряжения. М: Машгиз. 1963, 233 с.

2.     Быковцев Г.И., Ивлев Д.Д. Теория пластичности. – Владивосток: Дальнаука. 1998. 528 c.

3.     Клюшников В.Д. Физико-математические основы прочности и пластичности. – М.: Изд-во МГУ, 1994. 189 с.

4.     Ю.М. Темис. Прикладные задачи термопластичности и термоползучести. В энциклопедии «Машиностроение». Т.I-3. Кн. 1. Динамика и прочность машин. Теория механизмов и машин. М.: Машиностроение, 1994. С. 226 - 272.

5.     Термопрочность деталей машин. Под ред. И.А. Биргера и Б.Ф. Шорра. - М. Машиностроение. 1977 г. 455 с.

6.     Зарубин В.С., Кувыркин Г.Н. Математические модели термомеханики. – М.: Физматлит. 2002. 168 c.

7.     Кадашевич Ю.И., Новожилов В.В. Теория пластичности, учитывающая остаточные микронапряжения // ПММ. 1958. Т.22. Вып. 1. С. 78-89.

8.     Гузев М.А. Термомеханическая модель упругопластического материала с дефектами структуры. // Механика твердого тела. 1998. № 4. С. 341 – 360.

9.     Галанин М.П., Гузев М.А., Низкая Т.В. Разработка и реализация вычислительного алгоритма для расчета температурных напряжений, возникающих при нагреве металла, с учетом фазовых переходов. // Препр. Ин-та прикл. матем. им. М.В. Келдыша PAH. - 2005. - № 139. 19 с.