CC BY
31
УДК 535:621
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕПЛО- И ЭЛЕКТРОПЕРЕНОСА ПРИ ВОЗДЕЙСТВИИ СИЛЬНОТОЧНОГО ИМПУЛЬСА НА ЭЛЕКТРОД
Р. В. Арутюнян
Арутюнян Роберт Владимирович, кандидат физико-математических наук, доцент кафедры математического анализа, Московский технический университет связи и информатики, rob57@mail.ru
В статье исследовано влияние нелинейностей теплофизических параметров и фазовых переходов плавления и испарения на электрические и тепловые процессы при нагреве металлического электрода сильноточным импульсом. Сформулирована математическая модель, а также разработаны конечно-разностный метод и программы для ЭВМ, позволяющие эффективно осуществлять компьютерное моделирование тепло- и электрофизических процессов при воздействии сильноточного импульса на металлические электроды. Расчет полей осуществляется на основе сквозного энтальпийного метода. Осуществлена серия расчетов для информативного случая железа. Установлено значительное влияние нелинейностей теплофизических параметров, фазовых переходов плавления и испарения, вида краевых условий на значения температурного и электрического полей.
Ключевые слова: металлический электрод, сильноточный импульс, электрическое поле, температурное поле, расчет, сквозной энтальпийный метод.
DOI: 10.18500/1816-9791 -2016-16-2-138-144
1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
В теории электрических контактов, плазмотронов, электросварки и т.д является актуальной задача расчета электрического и температурного полей с учетом нелинейностей электро- и теплофизических свойств материала, фазовых переходов и других факторов [1-4]. В предлагаемой статье результаты отмеченных исследований развиваются на основе сквозного «энтальпийного» метода, позволяющего эффективно учитывать нелинейности теплофизических характеристик материала, фазовые переходы (плавления, испарения), радиационное и конвективное охлаждения поверхности материала.
Расчет тепло- и электрофизических процессов осуществляется на основе нестационарной пространственной осесимметричной модели. Предполагается, что электрод занимает полупространство г > 0, г — пространственная координата (аппликата) в цилиндрической системе координат. Ток величиной 10 проводится через токовое пятно радиусом Я0 на поверхности электрода (г = 0).
1.1. Задача Стефана
При моделировании рассматриваются три этапа теплового процесса: нагрев материала до температуры плавления (твердая фаза); нагрев расплава и дальнейшее проплавление твердой части материала (жидкая фаза); начало интенсивного испарения и кипения материала (фаза испарения и кипения).
В статье для целей компьютерного моделирования применялся численный метод сквозного счета [5,6], основанный на преобразовании многофазной задачи Стефана к «энтальпийному» виду с сосредоточенной теплоемкостью. Краевая задача для уравнения теплопроводности имеет вид
д^ = 1 д дг г д
д Л
— = г = 0, г > 0, г> 0; qs = ст(«4 - «0) + к(« - «о) - 5исп(«), (2)
« = «0, г> 0, г> 0, г = 0; « ^ и0, а/г2 + г2 ^ го, ^ = 0 (г = 0), (3)
дг
и
где ^ — энтальпия: ф(«) = / р(г>)е(г>)^г> + р1 «плгпл 1(« — «пл), 1(«) — единичная функция, и — тем-
ио
пература, «пл — температура плавления, гпл — скрытая теплота плавления, «0 — начальное значение температуры тела («0 < «пл), р = р(«) — массовая плотность, р1 («) — плотность твердой фазы,
дЛ
дг
д 2л
+ + ^,
г > 0, г > 0,
г > 0,
(1)
Р. В. Арутюнян. Математическое моделирование тепло- и электропереноса
соответственно р2(и) — плотность расплава, с = с(и) — удельная теплоемкость материала электрода,
и
Л — тепловой потенциал: Л(и) = / А(-) А(-и) — теплопроводность, Ь — время, г — полярный
ио
радиус, ст — постоянная Стефана - Больцмана, к — коэффициент конвективного теплообмена, Жо — энергия, а Ьр — время токового импульса, ду — плотность джоулевых источников: ду = 7, где 7 = 7(и) — электропроводность, р — электрический потенциал, дисп — удельные потери на испарение с поверхности материала: дисп = р2(и)гисп-исп, где гисп — скрытая теплота, а -исп — скорость фронта испарения:
р (и.)
- исп =
М
2пДи. Р2(ив)'
где р — давление насыщенного пара при температуре поверхности и. = и(г, 0, Ь), М — масса моля и Л — газовая постоянная. Давление насыщенного пара: р(и.) = ро ехр (^Нит — ^ТгО, где ро — давле ние газа окружающей тело среды, и
Ни Нив
температура кипения при давлении р0. Толщина испаренного
слоя 2исп(Ь) определялась как интеграл: 2исп(Ь) = / -исп¿Ь.
о
1.2. Модель электрического поля
Электрический потенциал является решением краевой задачи [1-4,7]:
1А /^г7 (и)
г дг V дг
д / + а; 7(и)
= 0, г > о,
2 > 0.
(4)
В исследовании рассматривались два вида краевых условий:
10,
I о
ПН2
г ^ Д, 2 = 0, г > До, ; = 0;
др
р = Ро, г ^ До, ; = 0; — = 0,
д2
г > До,
2 = 0.
(5)
(6)
В (5) предполагается, что ток в токовом пятне распределен равномерно. Соответствующим примером является воздействие электрической дуги или молнии на электрод. В (6) в токовом пятне электрический потенциал принимает постоянное значение, что соответствует задачам теории электрических контактов (как известно, линии тока перпендикулярны контактной площадке) [1-4]. Условие
ду дг
на бесконечности: р ^ 0 при Vг2 + 22 ^ го. На оси симметрии (г = 0) : ду =0.
2. КОНЕЧНО-РАЗНОСТНЫЙ метод решения
Дифференциальные соотношения краевой задачи приближенно заменяются системой конечно-разностных уравнений:
Qi.-j Qi
1
ЛЩ — А^1 1 — Л^ \
г+1-г------^ | +
к
г,г+1
к
г,г+1
+ -
к
' ЛР+1 _ ЛР+1 ЛР+1 _ ЛР+1 \ Лг,0 + 1 Лг,3 Лг,3 Лг,0-1 \ , в+1
+ дУ,г,3'
2,0
к
2,0 + 1
к2
г = 1,...,пг — 1, з = 1,..., п2 — 1, р = 0,1,...,
где т — шаг по времени, Ьр = рт, р — номер временного слоя; кг^ — шаг сетки по координате г, г — номер узла сетки координаты г, к2,0- — шаг сетки по координате 2, 3 — номер узла сетки пространственной координаты 2, кг^ = (кг,^ + кг,^+1 )/2, к2,0- = (к2,0- + к2,0+1 )/2, ир,0- — сеточная функция, иР0 ~ и(г^, 20, Ьр), г = 0,..., Пг, 3 = 0,..., Пг, р = 0,1,..., Qf,0• = Q(up,0•), Лр,0- = Л(ир,0-). Начальное условие запишется в виде: и°,0- = ио, г = 0,..., пг, з = 0,..., п2. Краевое условие на поверхности тела с первым порядком точности:
р р р ) и^,1 — и^,о Чо У
А(и^о -— = д.(и^о), г = 0, ...,п, р = 1, 2,
¿2,1
г
ггкг,г
1
Изв. Сарат. ун-та. Нов. сер. Сер. Математика. Механика. Информатика. 2016. Т. 16, вып. 2 со вторым порядком точности:
Р V 2, 1 + ' 2,2) р '2,1 р
= _-_____-_01 _ ___II _ ___-_____-____
г,° (2Лг>1 + '2,2)^2,2 гД (2Лг>1 + '2,2)Лг>2 г'2 (2Лг>1 + '2,2) А«0)'
Условие на оси симметрии с первым порядком точности: Ор 3- = Ор 3-; со вторым порядком точности (? = 1,...,пг - 1, р = 0,1,-..):
р _ ('г, 1 + 'г, 2)2 'г, 2 '2 ,1 'г,2
2'г,1 + 'г,2 Ор,3 2'г,1 + 'г,2 О 2,3
Условия на границах сеточной области, соответствующие условиям на бесконечности: О 3 = 0, ? = 1,..., п2 - 1, Ор,Пг =0, г = 0,..., Пг, р = 0,1, —
Система конечно-разностных уравнений, аппроксимирующая краевую задачу для электрического потенциала:
1 | ^ ^Р . ^г+13 - ^г,з _ ^^Р. ^г,з - ^г-1,з | + Гг'г,^ у » 'г,г+1 * ^ 'г,г /
+ ^ Лр 3 _ <3 ^Р <3 _ = 0
г = 1,..., пг _ 1, ? = 1,...,п2 _ 1, р = 0,1,...; 7Р • 1 + 7Р •
С = , = 7 «,-),
оР ,ЛР ,ЛР
/ = ^Р (Е2 + Е2 Ч Е . . = ^+1,7_^г-1 ,3 Е . . = ^+1,7_, 3
V, г , 3 = 7г , 3 (,Ег , г , 3 + Е2 , г , 3) , Ег , г , 3 = 2' . ' Е , г, 3 = 2'
Т,г ^'"2,3
Конечно-разностная аппроксимация краевых условий на границах области: Для случая (5):
Р <1 _ <0 /1°/М2), Гг ^ Я°,
'2,1 [0, Гг > л°,
Для случая (6): ^р,° = при Гг ^ Д°, ^р,° = , г = 0,..., Пг, р = 0,1,... . Условие на оси симметрии: ^р,3 = ^р,3, ? = 1,..., п2 _ 1, р = 0,1,... . Сеточные уравнения, соответствующие условию на бесконечности:
=0, г = 0, ...,пг, ^ ,3 =0, ? = 1,...,п2 _ 1 р = 0, 1,....
На каждом временном слое система конечно-разностных уравнений решалась методом Зейделя. Вначале осуществлялся расчет электрического поля, что позволяло рассчитать плотность джо-улевых источников. На следующем этапе для каждого внутреннего узла сетки система решалась методом Ньютона относительно сеточного значения энтальпии Ор+1 при известном Ор3, г = 1, ...,пг _ 1, ? = 1,...,п2 _ 1. По этим значениям вычислялась температура в тех же узлах сетки: -ир+1 = О-1 (Ор,+1), р = 0,1,..., где О-1 — функция, обратная к функции О(о), вычислявшаяся также методом Ньютона как решение нелинейного уравнения относительно вида: О(у) = г. Величина -ир+1 в конечно-разностной производной краевого условия при г = 0 выражалась через остальную часть выражения и также находилась в процессе итераций по методу Зейделя.
Контроль погрешности численного решения осуществлялся по значениям теплового и токового баланса (порядка 5%).
3. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ
В качестве материала рассмотрено железо, которое характеризуется нелинейностью теплофизиче-ских характеристик [8,9]. Размер расчетной области, длительность и величина токового импульса были выбраны типичными для электрических контактов: £г = 10-3 м, Х2 = 10-3 м, = 2.25 ■ 10-3 с, 1° = 500 А. Площадь и радиус токового пятна: 5° = 10-7 м2, = д/5°/п.
Р. В. Лрутюнян. Математическое моделирование тепло- п электропереноса_
Параметры конечно-разностного метода выбирались в соответствии с требованиями точности: Пг = 30, иг = 70, т = 10-6.
Область в районе токопроводящего пятна с размерами Ьг = 5 ■ 10-4 м, Ьг = 5 ■ 10-4 м покрывалась сеткой с меньшим шагом при количестве узлов: пг0 = 20, п20 = 50.
3.1. Результаты для первого варианта краевых условий (5)
На рис. 1 представлены серии графиков температурного поля в зависимости от пространственных координат (1 — зависимости от г, г фиксированное, 2 — зависимости от г, г фиксированное, шаг по времени 2.5 ■ 10-4 с). Из общего множества графиков на рисунках здесь и далее представлен каждый третий.
На рис. 2 представлены графики (в относительных единицах) зависимости от времени координат фронтов фазовых переходов (1 — среднее значение в токовом пятне, 2 — в начале координат, 3 — график максимальной температуры, 4 — гпл(£), 5 — гисп(£)). Наибольшие значения координат фронтов фазовых переходов: хпл(£г) = 80 мкм, хисп(£г) = 3.4 мкм. Скорости движения фронта испарения: "Уисп(¿г) = 0.02 м/с.
и, К Ъ 3500
икип
2800
2100
ипл
1400 700
41
и 9 ■1 2
\1 \1
г, о.е 1.0
0.2
0.4 0.6 0.8 г, г мм
0.5 1.0 1.5 2.0
Т, мс
Рис. 1. Графики профилей температурного поля
Рис. 2. Графики температур и фазовых переходов
На рис. 3 представлены графики, описывающие характер электрического поля в электроде. В работе исследовались также потери на радиационное излучение и конвективное охлаждение, которые составляют порядка 1% и являются пренебрежимыми. На рис. 3 и 4 1 — зависимости от г, г фиксированное, 2 — зависимости от г, г фиксированное.
/, мВ 85
68
51
34
17
/мВ1
803
643 1 2 482
0.20 0.40 0.60 0.80 г г, мм
321 161
0.20 0.40 0.60 0.80 г г, мм
б
Рис. 3. Зависимости потенциала на 2-м шаге (а) и на 2000-м шаге (б) по времени
1
0
а
0.20 0.40 0.60 0.80 г, мм
аб Рис. 4. Плотность тока на 2-м шаге (а) и на 2000-м шаге (б) по времени
Как следует из рис. 4, плотность тока относительно слабо зависит от времени.
3.2. Результаты для второго варианта краевых условий (6)
На рис. 5 представлены серии графиков температурного поля в зависимости от пространственных координат (с шагом по времени 2.5 ■ 10-4 с).
и,к
2585
2068
ипл
1551
1034
517
и,к 3500
2800
2100
ипл '
0.20 0.40 0.60 0.80 к г, мм
аб Рис. 5. Профили температурного поля на 200-м шаге (а) и на 2000-м шаге (б) по времени
2, о.е. 1.0
2.3 Т, мс
Рис. 6. Графики: 1 — электрического потенциала в токовом пятне, 2 — тока, 3 —максимальной температуры, 4 — ¿пл(Ь), 5 — ¿исп(Ь)
На рис. 5 графики серии 1 выражают зависимости от г при фиксированном г, 2 — зависимости от г при фиксированном г.
При относительно малых напряжениях за несколько миллисекунд происходит установление стационарного режима, при больших напряжениях происходит интенсивное кипение и испарение материала.
Наибольшие координаты плавления и испарения (рис. 6): хпл(£г) = 80мкм, хисп(£г) = = 1.7 мкм. Скорость движения фронта испарения: г>исп(£г) = 0.012 м/с.
Графики, описывающие характер электрического поля в электроде, представлены на рис. 7 и 8.
Р. В. Арутюнян. Математическое моделирование тепло- и электропереноса
/мВ 600
480 360 240 120
/мВ 600
480 360 240 120
0 0.20 0.40 0.60 0.80 г, г, мм 0 0.20 0.40 0.60 0.80 г, г, мм
а б
Рис. 7. Зависимости потенциала на 2-м шаге (а) и на 2000-м шаге (б) по времени: 1 — г фиксировано,
2 — г фиксировано
Л/мм2* 57017
45614
34210 А1
22807
11403
0 0.20 0.40 0.60 0.80 г г, мм
а
Л/мм2* 6933
5547 ■■
4160 |
2773
1387 -
0 0.20 0.40 0.60 0.80 г г, мм
б
Рис. 8. Плотность тока на 2-м шаге (а) и на 2000-м шаге (б) по времени: 1 — г фиксировано,
2 — г фиксировано
В отличие от варианта (5) плотность тока на рассматриваемом интервале времени изменяется на порядок. Наибольшие значения плотности тока и температуры достигаются на краях проводящей ток площадки. Это согласуется с известным фактом теории контактов, что сваривание контактов часто происходит по краям контактного пятна.
ЗАКЛЮЧЕНИЕ
Сформулирована математическая модель, а также разработаны конечно-разностный метод и программы для ЭВМ, позволяющие эффективно осуществлять моделирование тепло- и электрофизических процессов при воздействии сильноточного импульса на электроды. Результаты работы могут применяться в практике исследования и проектирования электрических аппаратов и других электротехнических устройств.
Библиографический список
1. Таев И. С. Электрические контакты и дугогаси-тельные устройства аппаратов низкого напряжения. М. : Энергия, 1973.
2. Ульрих Т. А. Математическое моделирование процесса контактной точечной сварки : автореф. дис. ... канд. техн. наук. Пермь, 2000. 15 с.
3. Абрамов Н. Р., Кужекин И. П., Ларионов В. П. Характеристики проплавления стенок металлических объектов при воздействии на них молнии // Электричество. 1986. № 11. С. 22- 27.
4. Борисенко П. А., Павлейно О. М., Павлей-но М. А. Методы численного решения нелинейных
нестационарных термоэлектромеханических контактных задач // Современные проблемы электрофизики и электрогидродинамики жидкостей : сб. тр. IX Междунар. науч. конф. СПб., 2009. С. 287291.
5. Самарский А. А., Моисеенко Б. Д. Экономичная схема сквозного счета для многомерной задачи Стефана // Журн. вычисл. матем. и матем. физ. 1965. Т. 5, № 5. С. 816-827.
6. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М. : Едиториал УРСС, 2014.
7. Weisenfels C., Wriggers P. Numerical modeling of electrical contacts // Computational Mechanics. 2010. Vol. 46, iss. 2. P. 301-314.
8. Теплофизические свойства расплавов. Справ.-информ. интернет-портал ebibl/umkd/Mamina, 2009. URL: http://files.lib.sfu-kras.ru/ebibl/umkd/ Mamina/u_lectures.pdf (дата обращения : 11.11. 2015).
9. Теплоемкость железа. Справ.-информ. интернет-портал «Лаборатория крупного слитка steelcast. ru», 2009. URL: http://steelcast.ru/iron_heat_ capacity (дата обращения : 11.11.2015).
Simulation of the Temperature and Electric Fields by High-current Pulse to the Electrode
R. V. Arutyunyan
Robert V. Arutyunyan, Moscow Technical University of Communications and Informatics, 8a, Aviamotornaya st., 111024, Moscow, Russia, 111024, rob57@mail.ru
The paper investigates the influence of nonlinearities thermophysical properties and phase transitions of melting and evaporation on the electrical and thermal processes at heating of a metallic electrode of high-current pulse. We formulate a mathematical model and develop a finite-difference method and a computer program, allowing effectively to carry out computer modeling of thermal and physical processes when exposed to high current pulse to metal electrodes. The article describes the results of the calculation of the fields based on cross-cutting enthalpy method. A series of informative calculation is implemented for the case of iron. The study found the following effects - a significant influence of the nonlinearity of the thermophysical parameters, phase transitions melting and evaporation, types of boundary conditions on the values of temperature and electric fields.
Key words: metal electrode, a high current pulse electric field, temperature field, calculation, cross-cutting enthalpy method.
References
1. Taev I. S. Elektricheskie kontakty i dugog-asitel'nye ustroistva apparatov nizkogo napri-azheniia [Electrical contacts and arcing devices low-voltage devices]. Moscow, Energiia, 1973 (in Russian).
2. Ul'rikh T. A. Matematicheskoe modelirovanie protsessa kontaktnoi tochechnoi svarki [Mathematical modeling of the resistance spot welding process]: Dis. ... kand. tekhn. nauk. Perm', 2000. 15 p. (in Russian).
3. Abramov N. R., Kuzhekin I. P., Larionov V. P. Characteristics of penetration of the walls of metal objects when exposed to lightning. Elektrichestvo [Electricity], 1986, no. 11, pp. 22-27 (in Russian).
4. Borisenko P. A., Pavleino O. M., Pavleino M. A. Metody chislennogo resheniia nelineinykh nestat-sionarnykh termoelektromekhanicheskikh kontakt-nykh zadach [Methods for the numerical solution of nonlinear transient termoelektromehanicheskih contact problems]. Sovremennye problemy elek-trofiziki i elektrogidrodinamiki zhidkostei [Modern Problems of Electrophysics and Electrohydro-
dynamics Liquids] : Proc. IX Intern. Sci. Conf. S.-Petersburg, 2009, pp. 287-291 (in Russian).
5. Samarskii A. A., Moiseyenko B. D. An economic continuous calculation scheme for the Stefan multidimensional problem. U.S.S.R. Comput. Math. Math. Phys, 1965, vol. 5, no. 5, pp. 43-58. DOI: 10.1016/0041-5553(65)90004-2.
6. Samarskii A. A., Vabishchevich P. N. Vychis-litel'naia teploperedacha [Computational Heat Transfer]. Moscow, Editorial URSS, 2014 (in Russian).
7. Weisenfels C., Wriggers P. Numerical modeling of electrical contacts. Computational Mechanics, 2010, vol. 46, iss. 2, pp. 301-314.
8. Thermal properties of melts. Reference and information portal ebibl/umkd/Mamina, 2009. Available at: http://files.lib.sfu-kras.ru/ebibl/umkd/ Mamina/u_lectures.pdf (accessed : 11.11.2015).
9. The specific heat of iron. Reference and information portal "Laboratory of large ingots steel-cast.ru", 2009. Available at: http://steelcast.ru/ iron_heat_capacity (accessed : 11.11.2015).
