Математическое моделирование процесса контактной односторонней сварки

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ  ПРОЦЕССА КОНТАКТНОЙ ОДНОСТОРОННЕЙ СВАРКИ

 

1 Разработка математической  модели процесса сварки

 

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

С целью некоторого упрощения модели был принят ряд допущений:

1. Для уменьшения объема  вычислений электрод заменяем  его частью высотой d_э (рис. 1.). На поверхности Z = d₂ + d₁ + d_э, которая содержит дно канала водяного охлаждения, полагаем наличие одинакового потенциала, считая его неизменным или возрастающим от нуля до заданной величины (при имитации модуляции переднего фронта импульса тока) по времени прохождения тока, а также постоянную температуру, равную температуре воды, используемой для охлаждения.

 

 

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

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

4. Учитывая условие  относительной скоротечности процесса  сварки, считаем количество теплоты,  излучаемое с поверхности деталей  и электродов, незначительным. Границы  деталь - воздух и электрод –  воздух справедливо считать адиабатическими стенками; для теплового поля граничными условиями на этих поверхностях будут условия 2-го рода, т. е. тепловой поток через эти границы постоянен и равен нулю.

Исходя из принятых допущений, рассмотрим область решения уравнений (рис. 4.1). Имеются две детали с размерами L ´ H₁ ´ d₁ (верхняя) и L ´ H₂ ´ d₂ (нижняя). Верхняя деталь контактирует с токоподводящим электродом диаметром D_э, высотой h_э и углом наклона рабочей части a_э. Центры контактов электрод – деталь и деталь - деталь привязаны к осям координат меняющимися размерами l_э, и Н₄. Диаметры контактов равны d_э. Размеры каждой детали варьируемые; параметр l_э равен половине расстояния между электродов. Смещение одной детали относительно другой фиксируется значением Н₃. Плоскость ZoY является плоскостью симметрии Электрического и температурного полей. Каждая деталь обладает своими электротеплофизическими свойствами, изменяющимися или неизменными в зависимости от температуры.

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

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

Электрическое и тепловое поле в области решения можно  описать соответственно стационарным и нестационарным уравнениями электрического потенциала и теплопроводности /109, 110, 111/, которые записываем в прямоугольной трёхмерной системе координат:

 

                            (1)

 

где:   – удельная электропроводность материала (См/м),

j_{(x, y, z)} – электрический потенциал точки поля (В),

x, y, z – её координаты (м).

 

                           (2)

 

где: Т – температура (К),

 t – время нагрева (сек),

 а_д – коэффициент температуропроводности материала (м²/с),

w_{(x, y, z, t)} – приращение температуры за счёт внутреннего источника тепла, определяемое по формуле:

 

                                                     (3)

 

где:   j – плотность тока (А/м²),

с – удельная теплоёмкость материала (Дж/(кг К)),

g - плотность материала (кг/м³).

Коэффициент температуропроводности материала равен:

 

                                                                   (4)

 

где:   l – теплопроводность (Вт/м К),

 с – удельная теплоёмкость материала (Дж/(кг К)),

 g - плотность материала (кг/м³).

Удельная электропроводность материала является величиной, обратной удельному электросопротивлению:

 

                                                      (5)

 

где: r– удельное электросопротивление материала (Ом м).

В зависимости от температуры удельное электросопротивление, теплопроводность, удельная теплоёмкость и плотность изменяются (рис. 2 и 3).

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

                                                     (6)

 

где: – вектор плотности тока (А/м²),

       – вектор напряженности электрического поля (В/м).

В свою очередь:

                                     (7)

и

                                     (8)

где: – единичные орты в прямоугольной системе координат,

     – проекция вектора напряженности на оси координат.

 

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

 

                                      (9)

Плотность тока равна:

                                         (10)

 

Начальные и граничные условия сведены в таблицу 1. Начальным условием для решения системы уравнений будет распределение потенциала при Т = Т₀. Для границ деталь – воздух и электрод – воздух первая производная потенциала и температуры по нормальным составляющим равны нулю (граничное условие 2-го рода) /112/. На плоскости симметрии деталей граничное условие для электрического поля j = 0, а на дне водоохлаждаемого канала электрода (при Z = d₂ + d₁ + d_э) - j = U / 2, где U – напряжение на электродах машины.

 

 

Таблица 1 – Начальные и граничные условия 

 

граница	поле
	электрическое	температурное
деталь-воздух
	j = 0
электрод-воздух
Z = d2 + d1 + dэ		T = T0 = Tводы
Z = d2	cz = 0	lz = 0
При t = t0	Распределение потенциала	T = T0

 

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

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

 

2 Особенности численного  решения модели

 

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

 

2.1 Разбивка области  на элементарные ячейки

 

Конечные элементы были выбраны в форме прямоугольного параллелепипеда, линейные размеры  которого зависят от местоположения элемента (рис.4.4). Это сделано с целью уменьшения количества элементов. Высокая точность расчетов требуется в зоне больших температурных градиентов, расположенной под электродом (рис. 5). Для удобства эта зона была выбрана в форме прямоугольного параллелепипеда размерами . В этой зоне размеры элементов минимальны и постоянны. Вне этой зоны линейные размеры элементов увеличиваются независимо друг от друга в геометрической прогрессии по мере удаления (по соответствующей координате) от зоны. Таким образом, все трехмерное пространство оказывается разбитым на элементы, каждому из которых в качестве номера присваивается трехмерный целочисленный индекс (i, j, к) относительно начала координат. По этому индексу определяется не только местоположение элемента в пространстве, но и его размеры и естественным образом индексы элементов, с которым он граничит. Область, как часть пространства, состоит из попавших в нее элементов. При этом элементы, частично выходящие за границы верхней и нижней детали, обрезаются границей, электрод же для удобства расчетов заменяется набором элементов, чьи центры попали в электрод (рис. 4). В обеих деталях увеличение элементов происходит по осям Х и У, а в электроде – по оси Z. Шаг прогрессии увеличения элементов задается вместе с другими исходными данными; в позволяющем большинстве случаев в расчетах использовался шаг 1.12. Зона наибольшей точности (зона под электродом) разбивалась на элементы таким образом, что их число по осям Х и У составляло 60. По оси Z число элементов в верхней детали равно 10; элементы нижней детали имеют такой же размер (по оси Z) и их число пропорционально толщине детали.

Вся область разбита  на три естественные части: нижнюю пластину, верхнюю пластину и электрод. Эти части имеют две области контакта относительно небольшой площади – круг диаметром d_э (малый диаметр электрода) между электродом и верхней пластиной и такой же круг между пластинами. Поэтому расчеты на каждом шаге проводились в два этапа. На первом этапе – в каждой из трех частей как в самостоятельной области; на втором – рассчитывалось влияние контакта на каждую часть. Каждая часть имеет свое дополнительное начало координат и свои дополнительные индексы элементов относительно него, что удобно при распределении машинной памяти и переборе элементов внутри части. При определении размеров элемента и индекса соседнего элемента из другой части в месте контакта используется основной индекс элементов.

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

Элемент имеет шесть граней, а  каждая грань принадлежит двум ячейкам, поэтому общее число граней в три раза больше, чем элементов. Грань элемента имеет 4 индекса: первые три от ячейки, и номер переменной, координатная ось которой ортогональна грани. Последний индекс принимает только три значения: 0,1 и 2. Два элемента, разделенные гранью, отличаются одним из трех индексов на единицу. Присвоим грани меньший из этих индексов. Так, между элементами (i-1,j,k) и (i,j,k) лежит грань (i-1,j,k,0), а между элементами (i,j,k) и (i,j+1,k) – грань (i,j,k,1).

 

2.2 Разностные уравнения

 

а. Уравнение потенциала

 

 

Рассмотрим элемент (i, j, k). Из элемента (i-1, j, k) в (i, j, k) течет ток, величина которого пропорциональна разности потенциалов в элементах, площади их совместной грани и обратно пропорциональна расстоянию между центрами элементов и удельному сопротивлению в центре грани.

 

                             (11)

 

где: - площадь грани,

      - расстояние между центрами элементов,

    - удельное электросопротивление в центре грани.

 

 

Коэффициент

                         (12)

 

характеризует электропроводность грани (i-1,j,k,0). Теперь формула (11) имеет вид

 

                           (13)

 

Аналогично находятся коэффициенты А и токи I, входящие в (i,j,k) через остальные пять граней. Сумма токов входящих в элемент, равна нулю; суммируя токи (13) по всем граням ячейки, приравнивая сумму нулю и выражая j_i,j,k, получаем разностное уравнение для потенциала:

 

     (14)

Здесь:

 

                            (15)

 

б. Уравнение плотности  тока

Плотность тока определяли как модуль векторной суммы составляющих плотностей тока:

         (16)

 

в. Уравнение теплопроводности

Количество теплоты, попавшее в элемент (i,j,k) из соседнего элемента (i-1, j, k) (т.е. через грань (i-1,j,k,0)) будет равно:

 

                (17)

 

Здесь: - температура в элементе (i,j,k) в начальный момент временного отрезка t,

        -коэффициент теплопроводности в центре грани.

Количество теплоты, попавшее в элемент через другие грани, находится аналогично. Кроме этого, в (i, j, k) выделяется теплота за счет прохождения электрического тока:

 

                            (18)

 

где   -плотность тока,

-удельное электросопротивление,

        -объем элемента.

В результате элемент  нагреется до температуры  . Количество теплоты, необходимое для этого, равно:

 

,                       (19)

где - плотность,

        - удельная теплоемкость.

Получаем уравнение

.      (20)

 

Подставляя в него (17), (18), (19) и выражая получаем:

 

         (21)

 

Здесь ,

 

         ,

 

         .

 

Коэффициенты В₃ …В₆ находятся аналогично.

Разностное уравнение (21) является явным, его отличают простота и экономичность вычислений. Однако, серьезным недостатком является его условная устойчивость. Поэтому на практике пользуются неявными разностными схемами, обладающими – безусловной устойчивостью /114/. Неявный вариант уравнения (21) выводится аналогично. Для этого в (17) используют температуру не в начальный, а в конечный момент временного отрезка t, т. е. . Неявное разностное уравнение теплопроводности имеет вид:

                                                      (22)

 

2.3 Решение уравнений  и организация вычислений

 

Уравнение (14) решаем методом итераций с помощью явной трехслойной схемы с постоянными параметрами:

 

,                              (23)

где - решение уравнения (14) для значения в левой части,

      с – номер итерации.

Формула (23) может быть использована при с>1. При с=1 необходимо знать и . В качестве были взяты начальные данные, а получено из начальных данных с помощью (14), т. е. = . Итерационный процесс останавливается при достижении условия

                                                 (24)

во всей области. Здесь e - заданная точность расчета.

Решение уравнения (22) прямыми методами неэкономично, поэтому при расчете использовали локально-одномерную схему, где вместо (22) решалось три одномерных уравнения:

 

      (25)

       (26)

      (27)

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

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

Начальными данными  для уравнения теплопроводности служит постоянная во всей области  температура 293 ⁰К.

 

 

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

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

(Т_пл-η; Т_пл+η), на котором разрыв заменяется линейной функцией. Число η выбирается таким образом, чтобы два смежных элемента, между которыми проходит граница раздела фаз, имели температуру из отрезка (Т_пл-η; Т_пл+η); мы использовали экспериментально подобранную η = 100⁰.

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

С учетом этого в системе  единиц СИ зависимости свойств от температуры примут вид:

 

, (21)

 

,                      (22)

 

, (23)

 

,                          (24)

 

, (25)

 

,                        (26)

 

,  (27)

 

 

.                      (28)

 

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

Расчет поля потенциала начинается с вычисления массива А по формулам (15), который остается неизменным во время всего итерационного цикла. Далее по формулам (14), (23) и (24) осуществляется итерационный процесс. Учитывая явный характер формул (14) и (23) вычисления проводятся во внутренних элементах каждой части области, затем в элементах контактов, и затем сносятся в граничные элементы из ближайших внутренних, кроме некоторых граничных элементов, где поле постоянно.

После нахождения потенциала по формуле (16) находится плотность тока во всех внутренних элементах области.

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

Общий алгоритм программы представлен на рис. 7.

Программа реализована  на языке C++ Builder 5.0. для Windows 95,98, Windows NT. Объем всего пакета около 5 мегабайт.

 

1. Ввод исходных данных.
2. Вычисление начальных  данных.		2.1. Разбивка области  на части и элементы.
		2.2. Определение принадлежности элемента электроду.
		2.3. Определение принадлежности  элемента контакту.
		2.4. Выделение массивов  под сеточные функции.
		2.5. Задание начальных  значений электрического потенциала.
3. Расчет потенциала		3.1. Расчет коэффициента А.
		3.2. Вычисление потенциала  в нижней, верхней пластине, в  электроде, в зонах контактов  и на границе.
		3.3. Проверка итерационного  процесса на завершение. При его  невыполнении переход к пункту 3.2.
4. Расчет плотности  тока во всех элементах.
5. Расчет температуры.	5.1. Расчет коэффициента В.
	5.2. Вычисление температуропроводности  всех элементов.
	5.3. Вычисление температуры  всех элементов.
6. Расчет изменения  электротеплофизических свойств  материала во всех элементах.
7. Проверка условия  окончания расчета по времени. При его невыполнении переход к пункту 3.
8. Запись в массивы  последних расчетных данных, их  сохранение в виде файлов с  расширением .txt. Вывод на экран.
9. Выход.

Рис. 7 – Общий алгоритм программы

 

 

3 Оценка адекватности  математической модели

 

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

Таким образом, суммарную ошибку можно определить по формуле:

 

,                                        (29)

 

где: S_мат – ошибка неточности численного решения,

S_y – ошибка неточности исходных данных,

S_упр – ошибка упрощения модели.

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

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

Для оценки математической точности использовали уравнение:

 

,                            (30)

 

где ∆у – ошибка,

Р – номер итерации.

Ошибка вычисления температуры  составила ∆Т_мат = 41⁰. При пересчете в выходные параметры модели – номинальный диаметр ядра (d _яном) и проплавление каждой из деталей (h_я1 и h_я2) – получили:

 

∆ d _{яном мат} = ∆Т_мат /grad_x Т = 0,33 мм

∆ h_{я1 мат} = ∆Т_мат /grad_z Т = 0,1 мм

∆ h_{я2 мат} = ∆Т_мат /-grad_z Т =0,16 мм.

 

Здесь grad_x Т, grad_z Т и -grad_z Т – проекции градиентов температур на границах литого ядра на оси электрода и в плоскости номинального диаметра ядра по соответствующим осям, которые определены при моделировании в момент выключения тока.

Основную часть суммарной  ошибки обычно составляет ошибка, связанная  с упрощением модели. В нашем случае, так как математическая модель создается  впервые, а не упрощается, необходима экспериментальная проверка ее адекватности. Для этого провели натурный эксперимент – сварку без подкладки образцов из Ст3 шириной 30 мм, толщиной 1+3 мм, расстояние между электродами 40 мм. Использовали следующие режимы сварки: I_св = 10 кА, (U_ээ /2 = 0,8 В), t_св = 0,3 с (15 периодов), F_св = 2500 Н, d_э = 7,5 мм.

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

Значимость расхождения  оценивали сопоставлением дисперсии  адекватности математической точности модели и точности эксперимента по параметру, имеющему наибольшее расхождение:

 

,                            (31)

 

где h_я2iэ – экспериментальное значение проплавления нижней детали в i – ом       эксперименте,

n – число измерений,

S_н – экспериментальное среднеквадратичное отклонение размеров литой зоны.

 

Таблица 2 – Значение теоретических и экспериментальных значений

 размеров литой  зоны

Параметр	Теоретическое	Среднее экспериментальное
dяном, мм	8,9	7,69	13,6
hя1, мм	0,9	0,83	8,4
hя2, мм	0,45	0,68	34

 

Значение F – критерия Фишера, равно 3,0284 и по данным /116/ соответствует доверительной вероятности (Р ≥ 0,95) того, что расхождение результатов моделирования с экспериментальными данными незначимо, а сделанные допущения приемлемы.

Необходимо отметить, что величина проплавления верхней  детали при моделировании определяется расстоянием от контакта электрод –  деталь до дна водоохлаждаемого канала, и изменением этого расстояния в  ту или иную сторону можно значительно повлиять на величину h_{я1 теор}.

Что касается проплавления нижней детали, то при моделировании  эта величина непостоянна (см. рис. 8), что объясняется отсутствием контактных сопротивлений в первую очередь. В табл. 2 при сопоставлении результатов использовано среднее значение h_{я2 теор}, что нашло свое отражение при расчете погрешности.