Содержание.

Введение.

§ 1. Математическая модель

§ 2. Численная реализация

            2.1 Описание алгоритма

            2.2 Результаты тестовых расчетов

§ 3. Анализ данных вычислительного эксперимента

Литература.

Введение

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

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

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

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

§1. Математическая модель

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

Уравнения равновесия в эйлеровых переменных записываются в обычном виде:

                                              (1.1)

Тензор напряжений зависит от внутренней энергии U (на единицу массы, далее используем и энтропию S на единицу массы) упругого тела:

                                                 (1.2)

Для описания состояния в качестве первого приближения возьмем модель:

               (1.3)

гдe ,  - удельная теплота плавления,  - удельная объемная теплоемкость при постоянной деформации. Типичный вид функции c_γ(T) приведен ниже на рис. 1.1.

Рис. 1.1

Выражение (1.2) дает:

                   (1.4)

В (1.3), (1.4) λ(T) и μ(T) – коэффициенты Ламе, β = (3 λ + 2 μ) α_T, α_T (T)- коэффициент линейного теплового расширения,  - коэффициент Пуассона. Характерные зависимости этих величин от температуры приведены на рисунках 1.2 и 1.3.

                   Рис. 1.2.                                                    Рис.1.3.

Тензор напряжений можно представить в виде:

,

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

Рис. 1.4

Усреднив полученные соотношения по толщине пластины и предполагая, что  (случай плосконапряженного состояния), получим следующее выражение для тензора напряжений:

Здесь и далее все индексы пробегают значения от 1 до 2.

Уравнение энергии для данного случая имеет вид:

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

Далее получим уравнение теплопроводности в следующем виде:

.             (1.8)

Оно является нелинейным относительно температуры T. Функции  являются разрывными функциями температуры и при численном решении будут заменены кусочно – линейными (см. рис. 1.2 и 1.3).

Данные о параметрах материалов и протекающих при сварке процессах взяты из работ [1 – 4].

§ 2. Численная реализация

2.1 Описание алгоритма

Итак, в случае малых перемещений имеем следующую систему уравнений:

         ,

         .

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

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

,

где  на каждом шаге ищется с помощью метода Ньютона:

,

Условия прекращения итераций:

     ,    .

Функция  имеет разрыв при , причем скачок . При численном решении ее можно заменить некоторой непрерывной функцией [6], в простейшем случае кусочно-линейной.

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

2.2 Результаты тестовых расчетов

Данный алгоритм реализован программно на языке C++ с применением библиотеки разреженных матриц SparseLib++. Для построения сеток использован встроенный алгоритм пакета MatLab 7.0. Расчеты проводились на трех сетках, которые отличаются друг от друга максимальным размером ребра .

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

                                            (2.1)

Задача решается в безразмерных переменных.

Во всех тестовых расчетах область представляет собой квадрат: . Далее используются следующие обозначения (см. рис. 2.1):

Рис. 2.1. Расчетная область

Для сравнения численного решения с точным используются следующие величины:

где индексы ex и h указывают на точное и приближенное решениие соответственно.

Тест 1. Свободное тепловое расширение

Задано однородное поле температуры:. Решается только задача механики с граничными условиями следующего вида:

На нижней границе заданы нулевые перемещения по y, на левой – нулевые перемещения по x. Остальная часть границы свободна.

Точное решение задачи:

.

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

Таблица 2.1. Погрешность решения задачи 1 на сетке с

Тест 2. Задача Стефана

Решается уравнение теплопроводности без источников тепла, но с учетом фазового перехода:

,

где

, .

На левой и правой границах задан перепад температур, остальная часть теплоизолирована:

.

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

.

При численном решении задачи вместо разрывной функции  используется кусочно-линейная функция:

,

где а  - параметр, отвечающий за ширину «зоны плавления». В данном расчете .

Рис. 2.2. ,                 Рис. 2.3. ,

На рис. 2.2, 2.3 показано положение «зоны плавления» для численного решения при t = 3 (установившийся режим) и различных отношениях . Черной линией показано положение фронта для точного решения.

Тест 3. Линейная задача термоупругости

Решается система уравнений:

,

где

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

Точное решение:

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

Таблица 2.2. Погрешность численного решения при

Таблица 2.3. Погрешность численного решения при

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

Тест 4. Задачи чистого сдвига и растяжения при наличии «жидкой фазы»

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

Рассмотрим задачу (2.2) с тепловым источником следующего вида:

и граничными условиями

,                

На верхней границе заданы ненулевые перемещения по x, нижняя граница закреплена. Перемещения по y запрещены на всей границе. Температура поддерживается постоянной и равной  на верхней и нижней границах, боковая часть границы теплоизолирована.

Параметры материала:

;

Индекс 1 соответствует твердому телу, индекс 2 – расплаву. Чтобы исключить влияние дополнительных факторов, все параметры, кроме m, оставлены непрерывными.

На рисунках 2.4-2.6 приведено поле перемещений в различные моменты времени. Также на этих рисунках показано положение фронта плавления (линия уровня T=1). Приведены результаты расчетов на двух сетках. Левый рисунок соответствует сетке с  правый – с .

Из этих рисунков видно, что до начала процесса плавления и на его начальной стадии поле перемещений линейно по y. При полном проплавлении пластины перемещения в верхней части равны заданным на границе, а в нижней части близки к нулю (см. рис 2.5, 2.6). Это означает, что верхняя часть пластины фактически теряет связь с нижней и движется как свободное твердое тело.

Рис. 2.4 t=0.02

Рис. 2.5 t=0.04

Рис. 2.6 t=0.06

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

                  .

Тепловой источник также выберем «диагональным»:

.

Результаты решения этой задачи представлены на рисунках 2.7-2.10.

                   Рис. 2.7 t=0.02                                 Рис. 2.8 t=0.04

                   Рис. 2.9 t=0.08                                 Рис. 2.10 t=0.1

Эти рисунки показывают, что качественно картина остается той же, что и в случае сдвига в направлении x. Искривленная форма фронта плавления объясняется условиями теплоизоляции, заданными на границе расчетной области.

Рассмотрим теперь задачу растяжения. Для этого зададим на границе области условия следующего вида:

          .

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

Как и ранее, нагрев осуществляется источником вида

.

t=0.04                                               t=0.06

Рис. 2.11

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

§ 3 Анализ данных вычислительного эксперимента

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

.

В качестве граничных условий зададим на боковых границах нулевые перемещения по x, а на верхней и нижней – нулевые перемещения по y:

.

Рассмотрим случай прямолинейного движения источника:

.

Параметры материала:

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

На рисунке, соответствующем моменту времени t=0.02 видно, что интенсивность напряжений достигает максимума на некотором расстоянии от центра пятна. В последующие моменты времени эта картина несколько искажается (максимум напряжений смещается к границе), что связано с условиями закрепления.

t=0.02

t=0.04

t=0.08

Рис. 3.1 Поле перемещений и распределение температуры (слева) и интенсивности напряжений (справа)

Литература

1. Теоретические основы сварки. Под ред. В.В. Фролова. М.: Высшая школа, 1970.

2. Сварка в машиностроении: Справочник. В 4-х т.т. Редкол.: Г.А. Николаев (пред.) и др. - М.: Машиностроенние, 1979. Т. 3. Под ред. В.А. Винокурова. 1979. 567с.

3. Физические процессы в металлах при сварке. Т.1. Элементы физики металлов и процесс кристаллизации. Н.Н. Прохоров. М.: Металлургия, 1968. 695 с.

4. В.Н. Волченко, В.М. Ямпольский, В.А. Винокуров и др. Теория сварочных процессов. – М.: Высшая школа, 1988. 559 с.

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

6. А.А. Самарский, В.Д, Моисеенко. Экономичная схема сквозного счета для многомерной задачи Стефана // ЖВМ и МФ. – 1965, т. 5, с. 816-827.