CC BY
17
Пентегов И. В. Pentegov I. V.
доктор технических наук, профессор, ведущий научный сотрудник отдела электротермии, Институт электросварки им. Е.О. Патона Национальной академии наук Украины, г. Киев, Украина
Рымар С. В. Rymar S. V.
доктор технических наук, старший научный сотрудник, ведущий научный сотрудник отдела электротермии, Институт электросварки им. Е.О. Патона Национальной академии наук Украины, г. Киев, Украина
УДК 536.12:621.78.01
ПРАКТИЧЕСКОЕ ПРИМЕНЕНИЕ МЕТОДА ТЕПЛОВЫХ ИСТОЧНИКОВ ПРИ АНАЛИЗЕ ТЕПЛОВЫХ ПРОЦЕССОВ В ЭЛЕКТРОТЕХНИЧЕСКИХ СИСТЕМАХ
Цель. Разработка алгоритмов интегрирования в модернизированном методе тепловых источников, обеспечивающих ускорение счета, а также сравнение времени, затрачиваемого на тепловые расчеты, в различных задачах.
Методология. Метод основывается на научных положениях теоретической теплофизики: теории распространения тепла в сплошных средах для определения распределения тепловых потоков в массивных телах, математических методах и способах интегрирования.
Результаты. Разработаны алгоритмы интегрирования в модернизированном методе тепловых источников, обеспечивающие ускорение счета. На примере пакета MathCAD показано, что скорость расчетов в предложенном методе, по сравнению с традиционными методами расчетов, увеличивается от нескольких десятков до нескольких тысяч раз в зависимости от сложности функции удельной мощности. Алгоритмы были проверены при решении ряда задач электротермии и показали высокую точность совпадения результатов расчета с опытными данными, а также с данными расчета прикладных компьютерных программ ведущих мировых фирм, использующих метод конечных элементов. Разработанные алгоритмы интегрирования в модернизированном методе тепловых источников позволяют решать задачи высокочастотного индукционного нагрева объемных тел и особенно эффективно тел с простыми геометрическими формами.
Оригинальность. Ускорение расчетов тепловых процессов с использованием разработанных алгоритмов интегрирования в модернизированном методе тепловых источников столь существенно, что позволяет рассчитывать распространение тепла в нагреваемом изделии в режиме реального времени.
Практическая ценность. Ускорение расчетов приводит к появлению нового качества метода — выбор режимов нагрева изделий может быть автоматизирован и встроен в технологическую линию по их производству.
Ключевые слова: тепловые расчеты, тепловая мощность, распространение тепла, модернизированный метод тепловых источников, алгоритмы интегрирования, повышение скорости расчетов, высокочастотный индукционный нагрев.
Electrical facilmes and systems
PRACTICAL APPLICATION OF THE HEAT SOURCES METHOD TO THE ANALYSIS OF THERMAL PROCESSES IN ELECTROTECHNICAL SYSTEMS
Purpose. Development of integration algorithms, ensuring accelerated computations in the framework of the modernized method of heat sources, and comparison of times spent on thermal calculations in various problems.
Methodology. The method is based on scientific principles of theoretical thermal physics — theories of heat propagation in continuous medium, used for determination of distribution of thermal flows in massive bodies, and also on mathematical methods and integration procedures.
Results. Algorithms of integration providing acceleration of computations are developed for the modernized method of heat sources. By the example of calculations in the MathCAD software package it is shown that the speed of calculations with proposed method, in comparison with traditional methods of calculations, increases by a factor of several tens to several thousand, depending on the complexity of the power density function. The proposed algorithms have been verified in the solutions of some problems of electrothermics and have shown high accuracy of obtained results in comparison with empirical data, and also with the data calculated in specialized software packages, engineered by the leading global developers and utilizing the method of finite elements. The developed algorithms of integration in the modernized method of heat sources allow solving problems of high-frequency induction heating of volumetric bodies and are especially effective in the case of simple geometrical forms.
Originality. The speed-up of calculations of thermal processes with the use of the developed integration algorithms in the modernized method of heat sources is so significant, that it allows conducting a real-time computation of the heat propagation in heated body.
Practical value. The acceleration of calculations leads to a new quality of the method — the choice of the products heating regimes can be automated and built into their manufacturing line.
Key words: thermal calculations, heat rating, heat propagation, modernized method of thermal sources, integration algorithms, acceleration of computations, high-frequency induction heating.
При анализе распределения тепловых потоков в задачах высокочастотного индукционного нагрева в электротехнических системах часто используется аналитический метод расчета тепловых источников Фурье [1], развитый Н. Н. Рыкалиным [2, 3], оперирующий конечной величиной импульса энергии. В этом методе в точке возникновения импульса энергии возникает бесконечно высокая температура, что является существенным недостатком метода, вносящим погрешность при расчете температур, особенно вблизи источника нагрева.
В работе [4] предложен новый способ суммирования радиус-векторов в методе тепловых источников, что позволило избежать бесконечно высоких температур в точке возникновения импульса энергии и тем самым существенно повысить точность расчетов, сделав существенный шаг в развитии аналитических методов расчета тепловых процессов. Этот модернизированный метод базиру-
ется на интегральном подходе при решении задачи распространения тепла. Однако взятие интегралов в компьютерных математических программных пакетах довольно медленный процесс, значительно снижающий скорость вычислений.
Существует ряд задач, для которых сокращение времени счета очень актуально. Это задачи тепловых расчетов при закалке и отпуске разнотипных изделий, таких как трубы, рельсы, профили, фланцы, а также сварка и гибка изделий в процессе серийного производства конструкций, нагрев электродов при различных видах сварки, определяющий режим сварки, и другие технологические процессы. Все эти расчеты желательно проводить в реальном времени, так как переналадку режимов нагрева необходимо осуществлять быстро при каждой смене типа и номенклатуры обрабатываемого изделия. На примере пакета MathCAD [5] будет показано, что при использовании предложенного в
настоящей статье способа интегрирования в модернизированном методе в ряде случаев расчетное время уменьшается в тысячи раз. Это приводит к появлению нового качества модернизированного метода — выбор режимов нагрева изделий может быть автоматизирован и встроен в технологическую линию по их производству.
Целью и задачей статьи является разработка алгоритмов интегрирования в модернизированном методе тепловых источников [4], обеспечивающих ускорение компьютерных расчетов, а также сравнение времени, затрачиваемого на тепловые расчеты в различных задачах. Необходимо дать рекомендации по использованию встроенных функций «erf» в пакетах математических программ, позволяющие отказаться от использования явных функций интегрирования и резко сократить время счета.
Исследование времени счета при различных алгоритмах вычисления кратных интегралов
Распределение приращений температур в неограниченном трехмерном пространстве при источниках тепла, равномерно распределенных внутри прямоугольного параллелепипеда со сторонами X, Y, Z и началом координат в центре параллелепипеда, при новом методе суммирования радиус-векторов описывается уравнением (17) работы [4], справедливого для случая отсутствия теплоот-вода вверх от полупространства:
xt2 (4=4
D Р 1 Г - 0 1 &1(x,y,z,t) = р" з|—Ц-l о
су(4тш)2 o(i-t)2_J/2
I е 4a^dy0 I е Aa(,-^dz0dx
-Г/2 -Z/2
(1)
точке с координатами (х0,у0,20) выделяется энергия dQ (эта энергия равна dQ = руётёУ, где ёУ — элементарный объем ёУ = ёх ^ёу в области распределенных источников тепла [4]); (1 - т) — время распространения тепла, связанного с выделившейся энергией dQ; (х, у, z) — координаты точки тела. В расчете обычно принимают средние значения тепло-физических параметров за время нагрева.
Истинная температура в точке (х, у, z) в момент времени 1 будет равна
Т1(Х, У, 2, г) = в(х, У, 2, г) + т.. (2)
Здесь Т0 — температура окружающей среды (начальные температуры всех точек нагреваемого изделия также равны Т0).
Здесь и далее номера цифровых индексов при параметрах температур в и Т соответствуют номеру определяющей в формулы.
Если же в уравнении (1) не разделять внутренние интегралы, то для приращений температур получим классическую формулу
четырехкратного интеграла: р
®з(*>У>2>0 =
су(4тш)2
* Z/2 Г/2 Х/2 (*>-*)2+<У<>-У)2+(го-4
WW
е
4и((-т)
dx0dy0dz0dz, (3)
где pV — удельная мощность, выделяемая в элементе объема, Вт/м3. Сначала будем полагать pV = const, так что эта величина может быть вынесена за знак интеграла; с — удельная теплоемкость материала, Дж/(кг-Х); у — плотность материала, кг/м3; а = X(cy) — коэффициент температуропроводности, м2/с; X — коэффициент теплопроводности, Вт/(м^); t — текущее время, отсчитываемое от начала процесса нагрева; т — момент времени, отсчитываемый от начала процесса нагрева, когда в
О -г/2-Г/2-ЛГ/2
где (х - х) + (у - у) + (2 - 2) — квадрат расстояния от точечного источника тепла в точке с координатами (ха у, г) до точки тела с координатами (х, у, 2), в которой определяется приращение температуры. Истинная температура в этой точке при использовании формулы (3) будет равна
Т(Х, У, 2, г) = в(х, У, 2, г) + т. (4)
Сравним время вычислений по формулам
(1) и (3).
Для хронометрирования вычислений воспользуемся встроенной в MathCAD функцией текущего времени (в секундах) йте(ф), где ф — любой идентификатор [5]. Эта функция требует двукратного вычисления: до и после расчетного фрагмента.
Значения параметров, использованных при хронометрировании
Рассмотрим решение гипотетической задачи поверхностного нагрева нижнего неограниченного полупространства (немагнитная сталь) прямоугольным высокочастотным индуктором с площадью Х2 и длиной Z, с цен-
тром в начале координат, расположенного над полупространством (см. рисунок 1), на котором У/2 — глубина проникновения тока в массивное тело, при отсутствии потерь тепла с поверхности тела в верхнее полупространство на излучение и конвекцию. При этом верхняя половина параллелепипеда с размерами X, У/2, I должна рассматриваться как «зеркальное отражение», которое «запирает» для прохода тепловых потоков верхнюю границу нижнего полупространства и обеспечивает требуемые граничные условия [4]. Размеры высокочастотного индуктора примем равными: ширина настила тока X = 0,1 м, длина I = 0,5 м. Частота тока индуктора / = 1,6 кГц, угловая частота ю = 2п/ = 104 рад/с. Средние значения теплофизических параметров стали X = 47,14 Вт/(мК), с = У^а) = 572,5 Дж/(кгК); а = 1,05 10 5 м2/с, g = 7860 кг/м3. Среднее значение удельного электрического сопротивления стали р = 6,45 10—7 Омм. Глубина проникновения магнитного поля в сталь Д = У/2= = ^рДюцо) = 0,01 м (здесь уже учтено зеркальное отражение Д (У = 2Д) вверх от плоскости _у0 = 0 для запирания этой границы для потока тепла вверх [4]); магнитная постоянная ц 0 = 4п10-7 Гн/м.
Удельная поверхностная мощность тепловло-жения постоянна и равна р = 600 кВт/м2, эквивалентная, равномерно распределенная по глубине проникновения удельная мощность объемного тепловложения ру= 2 р5 /У = 60 МВт/м3.
Активная мощность индуктора Р = р8Х1= = ргХУ1/2= 30 кВт.
Начальная температура Т0 = 293 К. Время нагрева 1п = 90 мин = 5400 с.
Хронометрирование
Д =172
Рисунок 1. Геометрия области с источниками тепла под индуктором (сечение по = 0, индуктор не показан)
Для хронометрирования использовался персональный компьютер с процессором Intel Core i5-3450.
Расчетный фрагмент представляет собой температуры массива точек
х,. = х - +-
х„„„ - х„
i); У] Уш
+ —{у
>
\ шах лип /' * J * шш \
где х . = -0,2; x = 0,2; y . = -0,4; y = 0;
mm y 5 max y 5 s mm 5 5 s max 5
i = 1, ..., N; j = 1, ..., N.
Вычисления проводились для сечения z = 0 при N = 80 по следующей схеме. Расчетный блок для интеграла (1) в пакете MathCAD имеет вид
t0 =time(0);
(5)
Здесь t0 — время начала формирования матрицы температуры F1 ; t1 — время конца формирования матрицы температуры; At1 — время счета матрицы температуры по формуле (1).
Аналогично, расчетный блок для интеграла (3) будет иметь вид
t2 = time(2); /F3,J.=©3(xI.)y,,0,i„);
(6)
Ai3 = г3 -i2.
Здесь t2 — время начала формирования матрицы температуры F3 ; t3 — время конца формирования матрицы температуры; At3 — время счета матрицы температуры F3, Величины At1 и At3 зависят от типа компьютера, но отношение этих величин меняется слабо и его можно использовать для сравнения быстроты счета по формулам (1) и (3). Вычисления по формулам (5) и (6) показывают, что это отношение лежит в интервале At J At х = 50 - 60. (7)
Таким образом, один только переход к новому методу суммирования радиус-векторов в методе тепловых источников позволяет сократить расчетное время в 50-60 раз по сравнению с классическим методом.
Достоинства нового метода этим не ограничиваются. При разделении внутреннего интеграла в уравнении (3) на произведение трех интегралов в уравнении (1) открывается возможность еще во много раз увеличить скорость вычислений.
Электротехнические комплексы и системы
Для случая, когда ру зависит от времени и не является функцией координат, запишем уравнение (1)в виде
I xji
4а(1-т)
, -Жо
Jm(t-i)
т/2 _iZoZZ)L г/2 (го"г)2
Г 6 4>~Z) jy [
J JnaU-i) J
4u(i-x)
=dz„dx
(8)
_у/2 х/Й^Г70 ^тш^-У где ру (т) — любая функция времени.
Внутренние интегралы в уравнении (8) безразмерны и, как показано в работе [4], выражаются через встроенные в МаЛСДС интегралы вероятности:
Х/2
4a(i-x)
-Х/2
I A=erf
^na(t-T)
-erf
х-
X
2 yla(t-x)
II2
4о(!-Т)
-7/2
^na(t-i)
dy0 = erf
Z/2
-erf
kri
4o(i-x)
2ja(t-x)
-Z/2
e_^
^Tia(t-x) °
-erf
= erf
X+2 2 sla(t-t)
= Erf(jc,X,a,i,T);
У " 2 Va(f-t)
= Erf(^,r)e,i,T); Z
z + —
(9)
(10)
2 yja(t-z)
2ja(t-%)
= Erf(z,Z,e,i,T).
(11)
тело, при отсутствии отвода тепла вверх от полупространства (см. рисунок 1). Активная мощность, передаваемая от индуктора в нижнее полупространство в момент времени т, при этом будет равна ру (т)ХУ1/2.
Вычисление этой функции не связано с взятием четырехкратного интеграла и время вычисления чрезвычайно мало.
Расчетный блок для интеграла (12) имеет вид
Г4 =йше(4);
t5 =time(5);
(13)
Ain=t5-
В дальнейшем будем применять краткие обозначения Erf вместо разностей интегралов вероятности, и везде, где это возможно, переходить от интегрирования к встроенным функциям erf. Произведя замены в уравнении (8), получим приращение температуры в точке (x, y, z):
®n(x,y,z,t) = -^^pr(x)-Vxi{x,X,a,t,%) x 0
x ^ri{y,Y,a,t,x)-m(z,Z,a,t,i)dx. (12) Этот интеграл позволяет рассчитывать температуру в процессе нагрева массивного тела (неограниченного полупространства) сверху плоским высокочастотным индуктором, у которого Z — длина индуктора, X — ширина настила тока под индуктором, Y/2 — глубина проникновения тока в массивное
Здесь t4 — время начала формирования матрицы температуры F12 ; t — время конца формирования матрицы температуры; At — время счета матрицы температуры по формуле (12).
Вычисления показывают, что отношения времен счета интегралов (1) и (3) ко времени счета интеграла (12) лежат ориентировочно в интервалах:
Atj/At12 = 26-34; At3/At12 = 1300-2040. (14) Таким образом, переход к встроенным функциям erf позволяет сократить расчетное время более чем на 3 порядка по сравнению с классическим методом, см. уравнение (3). Это качественный скачок.
Если удельная мощность pV , выделяемая в элементе объема, зависит от одной из координат, например pv = pV (х0,т), то интеграл (12) может быть записан в виде
()2
t Х/2
О-Х/2
Jna(t-x)
i£c0 x
х Erf(.y,7,a,í,т)•Erf(z,Z,a)í)т)í/т. (15) Здесь внутренний интеграл в общем случае не может быть выражен через интегралы вероятности, и время счета Д1 по сравнению с Д возрастает в зависимости от сложности функции ру (х0,т) приблизительно на порядок.
Если удельная мощность ру , выделяемая в элементе объема, зависит от двух координат, например ру = ру (х0,у0,т), то интеграл (15) может быть записан в виде
(\-х)1+(угу)2
t 7/2 Х/2
Mw-O^f J J
РГ(Х0,У0,Ф 4fl(("T>
О-7/2-Х/2
m(t-x)
-dx0dy{ x
х Ш{г,г,а,г,т)с1т: . (16)
Здесь время счета Дtl6 по сравнению с Дtl5 возрастает в зависимости от сложности функции ру (х0у0,т) приблизительно в два раза.
И, наконец, если удельная мощность ру , выделяемая в элементе объёма, зависит от трёх координат, ру = ру (х0у0,г0,т), то интеграл (16) преобразуется в
©17 (*>7.2.0 =
I г/2 7/2 ЛГ/2
су! 111
О -Z/2 -У/2 -Х/2
pF(x0, y0,z0,x)e 4°<"> (4 mf\t-xf11
-dx0dy0dz0di.
(17)
Это аналог формулы (3) с тем лишь различием, что здесь время счета At17 по сравнению с At3 возрастает в зависимости от сложности функцииpV (х0, y0, z0, т) приблизительно в 4-5 раз, и по сравнению с At получаем: At17/At3 = 4-5; At17/At12 = 5000-10000. (18)
Заметим, что при каждом вычислении At нужно делать полный пересчет программы.
Аналогичный результат наблюдается при пределах интегрирования, являющихся функциями других параметров, например, если Z = fX) и т.д.
Таким образом, разбиение многократного интеграла на независимые интегралы сокращает время счета в 50-60 раз.
При использовании в модернизированном методе с независимыми интегралами встроенной функции «erf» пакетов математических компьютерных программ время счета, по сравнению с многократным интегралом, сокращается в 1300-2000 раз.
Если удельная мощность, выделяемая в элементе объема, зависит от трех координат, то перейти к встроенным функциям не удается, и время счета по сравнению со временем счета с использованием встроенных функций «erf» возрастает в зависимости от сложности функции удельной мощности приблизительно в 10000 раз.
Наиболее существенное уменьшение времени счета происходит при разбиении области с источниками тепла на элементарные параллелепипеды с равномерным распределением источников тепла внутри каждого параллелепипеда, отдельного расчета тепловых процессов для каждого параллелепипеда и использования принципа наложения темпе-
ратурных полей, что допустимо при усредненных и постоянных значениях теплофизи-ческих параметров.
Во всех этих случаях для ускорения счета необходимо использовать метод суммирования (наложения) тепловых вкладов от элементарных параллелепипедов, на которые приближенно разбивается объем с источниками тепла. Этому методу будут посвящены дальнейшие исследования.
Авторы выражают благодарность кандидату физико-математических наук Всеволоду Игоревичу Пентегову за полезные замечания, сделанные им при просмотре статьи.
Выводы
1. Разработаны алгоритмы интегрирования в модернизированном методе тепловых источников, обеспечивающие уменьшение времени счета. На примере пакета МаШСАО показано, что в ряде случаев расчетное время уменьшается в тысячи раз, по сравнению с традиционным интегрированием многократного интеграла.
2. Ускорение расчетов тепловых процессов с использованием предложенных алгоритмов интегрирования в модернизированном методе тепловых источников столь существенно, что позволяет рассчитывать распространение тепла в нагреваемом изделии в режиме реального времени. Это приводит к появлению нового качества метода — выбор режимов нагрева изделий может быть автоматизирован и встроен в технологическую линию по их производству.
3. Предложенные алгоритмы интегрирования в модернизированном методе тепловых источников позволяют решать задачи высокочастотного индукционного нагрева объемных тел, и особенно эффективно тел с простыми геометрическими формами.
4. Разработанные алгоритмы были проверены при решении ряда задач электротермии и показали высокую точность совпадения результатов расчета с опытными данными, а также с данными расчета прикладных компьютерных программ ведущих мировых фирм, использующих метод конечных элементов.
Список литературы
1. Fourier J.B. Théorie analytique de la chaleur. Paris: Chez Firmin Didot, père et fils, 1822. 639 p.
2. Рыкалин Н.Н. Расчеты тепловых процессов при сварке. М.: Машгиз, 1951. 297 с.
3. Фролов В.В., Винокуров В.А., Волченко В.Н. и др. Теоретические основы сварки / Под. ред. В.В. Фролова. М.: Высшая школа, 1970. 592 с.
4. Пентегов И.В. К теории метода тепловых источников, используемого при анализе тепловых процессов в электротехнических системах // Электротехнические и информационные комплексы и системы. 2014. Т. 10. № 3. С. 5-5.
5. Кирьянов Д.В. MathCAD 14. СПб.: БХВ-Петербург, 2007. 704 с.
References
1. Fourier J.B. Théorie analytique de la chaleur. Paris: Chez Firmin Didot, père et fils, 1822. 639 p.
2. Rykalin N.N. Raschety teplovykh protsessov pri svarke. M.: Mashgiz, 1951. 297 s.
3. Frolov V.V., Vinokurov V.A., Volchenko V.N. i dr. Teoreticheskie osnovy svarki / Pod. red. V.V. Frolova. M.: Vysshaja shkola, 1970. 592 s.
4. Pentegov I.V. K teorii metoda teplovykh istochnikov, ispol'zuemogo pri analize teplovykh protsessov v elektrotekhnicheskikh siste-makh // Elektrotekhnicheskie i informatsionnye kompleksy i sistemy. 2014. T. 10. № 3. S. 5-15.
5. Kir'ianov D.V. MathCAD 14. SPb.: BKhV-Peterburg, 2007. 704 s.
