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

Количество теплоты, попавшее в элемент через другие грани, находится аналогично. Кроме этого, в (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 теор}, что нашло свое отражение при расчете погрешности.