CC BY
75
УДК 658.52
DOI: 10.15827/0236-235X.126.313-317
Дата подачи статьи: 01.10.18 2019. Т. 32. № 2. С. 313-317
Программы моделирования температурных полей в изделиях цилиндрической формы
Б.И. Марголис 1, д.т.н, зав. кафедрой, [email protected] Мансур Губран Али 1, аспирант, [email protected]
1 Тверской государственный технический университет,
кафедра «Автоматизация технологических процессов», г. Тверь, 170026, Россия
В статье рассмотрена модель температурного поля в изделии цилиндрической формы при несимметричном конвективно-радиационном теплообмене поверхностей изделия с окружающей средой и ограждающими поверхностями (нагревательными элементами) технологического оборудования. Получены соотношения для расчета температурного распределения в изделии с использованием численных конечно-разностных методов. Приведен пример расчета конвективно-радиационного охлаждения цилиндрического тела. Сформулирована возможность решения поставленной задачи в среде программирования Matlab.
На основе стандартной функции pdepe в среде MatLab разработана программа, позволяющая по заданным теплофизическим характеристикам материала (коэффициентам теплопроводности, температуропроводности) и параметрам конвективно-радиационного теплообмена (коэффициентам конвективной теплоотдачи и приведенным степеням черноты) моделировать температурное поле в изделии. Рассмотрены особенности разработки программы, связанные с заданием функций типа дифференциального уравнения, начальных и граничных условий. Приведены основные функции разработанных программ и результаты расчета температурного распределения.
Сделан сравнительный анализ решения задачи в среде Matlab с результатами, полученными с помощью конечно-разностных соотношений, продемонстрировано их хорошее совпадение. Рассмотрено решение двухмерной задачи теплопроводности в цилиндрическом теле с помощью функции parabolic, являющейся одним из инструментов Matlab для решения задач параболического типа. Произведена триангуляция цилиндрической области с использованием функции initmesh. Приведены фрагменты программных кодов функций граничных условий и визуальное изображение двухмерного температурного распределения в полом конечном цилиндре. При решении тепловой двухмерной цилиндрической задачи градиенты температур по высоте практически отсутствуют, в связи с чем наблюдается отличное совпадение температурного поля с результатами, полученными конечно-разностными методами и с помощью функции pdepe.
Разработанные методы в перспективе позволяют перейти к решению задач моделирования температурного поля в сортовых стеклоизделиях цилиндрической формы (стакан, бутылка, колба), имеющих, кроме боковой, еще и донную поверхность.
Ключевые слова: температурное поле, изделие цилиндрической формы, конвективно-радиационный теплообмен, конечно-разностные методы, моделирование, начальные и граничные условия, дифференциальное уравнение, метод конечных элементов, триангуляция, двухмерная задача теплопровод-
Математическое моделирование температурных полей в стеклоизделиях играет важнейшую роль при расчете процессов отжига стекла [1]. Геометрическая форма изделий может быть различной - плоской, цилиндрической. В таком случае одним из подходов при решении температурных задач является приведение сложного изделия к плоской форме за счет введения характерного или эффективного размера [2]. Однако при этом точность расчетов ниже и не учитывается реальная геометрия тел цилиндрической формы. В данной работе рассматриваются подходы к моделированию тем-
пературных полей в полых цилиндрических изделиях [3].
Для получения температурного поля необходимо решать дифференциальное уравнение теплопроводности в цилиндрических координатах:
dt(r,т) дт
d2t(r,т) 1 dt(r,т)
dr
dr
(1)
где ^г, т) - значение температуры в координате по радиусу г в момент времени т (°С); a - коэффициент температуропроводности материала изделия (м2/с) [4].
ности.
= a
Большую роль для точного моделирования температурных распределений в изделиях играет корректный учет начальных и граничных условий, которые в большинстве стадий технологических процессов термообработки имеют сложный характер. Так, начальные условия обычно неравномерны, граничные условия соответствуют несимметричному конвективно-радиационному теплообмену поверхностей изделия с окружающей средой и ограждающими поверхностями (нагревательными элементами) технологического оборудования [5, 6]. В этом случае наиболее удобно применять численные конечно-разностные методы, при использовании которых частные производные в уравнении (1) по времени, координате 2-го и 1-го порядков заменяются на их численные аналоги: [и+^) - иа)]/кг, [^ - 1) - 2^) + + + 1)]/(Дг)2 и [и(/ + 1) - иа)]/Дг. Здесь , = 0, т - номер временной точки расчета; у = 0, N - номер точки по толщине; Дт - шаг по
времени; Дг = (Л2 - Я{)М - толщина слоя; R1, Я2 - радиусы внутренней и внешней поверхностей полого цилиндра.
Тогда формулы для расчета температурного распределения во внутренних точках полого цилиндрического тела будут следующими:
у)=, у)+дт- а (дг )2 х
X {[*(у -1) - 2*, (у) +1 (у +1)] + [^ (у) -1 (у)]/у}, (2)
у = 1, N-1.
При расчете конвективного теплообмена необходимо приравнять тепловые потоки на поверхности, например, г = R1, в соответствии с законами Фурье - q(Rl, т) = -Х-д^1, т)/дг и Ньютона - q(Rl, т) = а,1[(^1, т) - ^(т)] [4]. Градиент д^1, т)/дг = а1[(^(т) - t(Rl, т)]/Х, производные 2-го порядка О (^) =д2 г(^, т)/дг2 = = 2а [^ (т) - * (0)]/(ХДг) на внутренней г = Rl и внешней г = R2 поверхностях: Ос(Я2) = д2 г (Я2, т)/дг2 =
= 2а2 [гс2 (т) - г, (N)]/(ХДг),
где а1, а2 - коэффициенты конвективной теплоотдачи для внутренней и внешней поверхностей изделия (Вт/(м2град.)); ^(т), ^(т) - температура среды внутри и снаружи (°С); X - коэффициент теплопроводности материала (Вт/(м-град.)).
При расчете радиационного теплообмена с учетом закона Стефана-Больцмана по аналогии с (3) вторые производные на поверхностях
r = R\ и r = R2 находятся по выражениям:
Q (R) = 28лст[т>)- T4(0)]/(XAr);
Q (R2) = 28,2а[Гя42(т)-T/(N)]/(XAr),
где 8s\, 8s2 - приведенные степени черноты верхней и нижней ограждающих поверхностей технологического оборудования; ст = 5,6710-12 (Вт/(м2К4)) - постоянная Стефана-Больцмана; Ти(х), Тй(х) - температуры ограждающих поверхностей над и под изделием (K).
C учетом передачи тепла точкам поверхностей изделия от внутренних соседних точек за счет теплопроводности получим окончательные выражения для расчета температур поверхностей цилиндрического тела в виде: Y+1 (0) = t. (0) + Ах- a х
х [ (1) -1 (0) )/Ar - (2/ Ar + 1/ R)+ QCi (R) + QK (R)]; 1 (5)
t.+1 (N) = t, (N) + Ах- a х
х [ (N -1) -1 t (N) )/Ar - (2/Ar +1R) + Qc (R) + QK (R)] ■
Для проверки модели температурного поля (2)-(5) был произведен расчет конвективно-радиационного охлаждения полого бесконечного цилиндра при следующих значениях параметров: a = 0,2 (см2/мин.); X = 0,00838 (Вт/(смтрад.)); R\ = 0,3 (см); R2 = 0,6 (см); = £.й = 0,91; N = 14; a\ = а,2 = 0,001256 (Вт/см2трад.); to(j) = 707 (°C); tc\ = tC2 = 27 (°C); ki= ^ = 27 (°C); Ах = 0,03 (с). Результаты расчета температуры в цилиндрическом теле приведены в таблице 1.
Таблица 1
Распределение температуры в точках по толщине неограниченного полого цилиндра
Table 1
Temperature distribution at points along the thickness of an unlimited hollow cylinder
Температура Время расчета х, с
0 6 10 20 30 60 120 œ
Внутренняя i(0) 707,0 619,6 584,2 511,8 456,4 344,8 226,8 27,0
Средняя t(N/2) 707,0 648,7 608,8 529,7 470,2 352,6 230,5 27,0
Внешняя t(N) 707,0 609,6 574,6 504,4 450,5 341,4 225,1 27,0
Решение поставленной задачи также может быть получено с помощью стандартной функции pdepe среды программирования Mat-lab [7, 8]. Характерные особенности разработки такой программы приведены в работе [9]. Функции типа дифференциального уравнения pdedif, начальных условий pdebeg и граничных условий pdebound не изменяются. Для учета
цилиндрической формы тела необходимо задать константу формы m = 1.
Результаты расчета температурного поля в неограниченном полом цилиндре с использованием стандартных функций Matlab приведены в таблице 2 и на рисунке 1. Разница с результатами, полученными с помощью конечно-разностных соотношений (2)-(5), не превышает 1 °C.
Для повышения точности моделирования процесса отжига необходимо учитывать конечные размеры стеклоизделий цилиндрической формы, в связи с чем возникает задача перехода от одномерного температурного поля по толщине изделия к двухмерному распределению, учитывающему градиенты температуры по высоте изделия. В пакете Matlab имеется набор инструментов для решения дифференциального уравнения теплопроводности в пространственных цилиндрических координатах, в частности, функция parabolic [7, 8]. Функция parabolic - решатель многомерных краевых задач методом конечных элементов [10] параболического дифференциального уравнения в частных производных:
дu(r, z, х) , ч
d--V-(cVu(r, z, х)J + au(r, z, х) = f,
дх
которое для случая осесимметричной температурной задачи в цилиндрической системе координат запишется следующим образом:
dt (r, z, х)
Рис. 1. Температурное поле в неограниченном полом цилиндре Fig. 1. Temperature field in an unlimited hollow cylinder
Таблица 2
Распределение температуры в неограниченном полом цилиндре
Table 2
Temperature distribution in an unlimited hollow cylinder
Температура В ремя расчета x, с
0 6 10 20 30 60 120 œ
Внутренняя t(0) 707,0 620,7 585,4 512,9 457,4 345,6 227,3 27,0
Средняя t(N/2) 707,0 649,2 609,4 530,3 470,8 353,2 230,9 27,0
Внешняя t(N) 707,0 608,3 573,5 503,8 450,2 341,4 225,2 27,0
d
dx
d 2i (r, z, x)
1 dt(r, z, x) д t(r, z, x) +------ + ■
dr
dz1
(6)
dr1
+ alt (r z х) = f,
где z - координата по высоте изделия; d = X/a, c = X - коэффициенты при производных; a1 = 0, f = 0 - отсутствующие множитель при температуре и внешний тепловой поток.
Для решения цилиндрической задачи с использованием parabolic необходимо произвести триангуляцию выбранной геометрической области с помощью функции initmesh. Обращение к этой функции в Matlab выглядит следую-
щим образом: [p, e, tr] = initmesh(g, 'Hmax', hmax), где g = decsg([3 4 - H/2 H/2 H/2 - H/2 R1 R1R2 R2]') - описание геометрической области цилиндрического тела в прямоугольной системе координат; hmax - максимальный размер граничного элемента. Функция initmesh формирует выходные параметры: p - массив горизонтальных (по высоте) и вертикальных (по радиусу) координат узлов конечно-элементной сетки; e - матрица граничных элементов, содержащая номера начальных и конечных узлов этих элементов и номера граничных сегментов, которым они принадлежат; tr - матрица треугольных конечных элементов, содержащая номера входящих в них вершин (узлов).
Визуальное изображение разбиения цилиндрической области на треугольные элементы, полученное с помощью функции pdeplot(p, e, t), показано на рисунке 2.
Обращение к решателю parabolic имеет синтаксис: uv = parabolic(u0, tauv, b, p, e, tr, c,
Рис. 2. Триангуляция полой конечной цилиндрической области
Fig. 2. Triangulation of a hollow finite tube domain
a\, f, d), где uv - температурное поле в моменты времени tauv, для которых вычисляется решение уравнения (6); u0 - узловое распределение температуры в начальный момент времени х = 0; b = @boundaryFileHolParabolic - имя функции, вычисляющей матрицы описания граничных условий.
Граничные условия, при которых осуществляется решение задач параболического типа, могут быть 1-го (Дирихле) и 3-го (Неймана) рода и описываются уравнениями hn = г и п ■ (сУ и) + qu = g соответственно. Для рассматриваемой задачи конвективно-радиационного теплообмена граничные условия 1-го рода отсутствуют, поэтому коэффициенты h = r = 0.
Граничные условия 3-го рода Неймана для температурной задачи:
-X-dt(r, х)/Эг + qt(r, х) = g. (7)
В то же время для конвективно-радиационного теплообмена на поверхностях цилиндрического тела (например, внутренней) по законам Ньютона и Стефана-Больцмана [3, 4]: - X -дТ(Я,,х)/дг = a, T(х) -T(R,,х)] +
Г 4 4 1 (8)
+ 8,1ст[Ti (х) - Т4(Ri, х) ] ■
В связи с этим с учетом перехода к абсолютным температурам коэффициенты q, g в уравнении (7) для внутренней цилиндрической поверхности записываются в виде q = a Т (R, х) + + 8,1 стТ4 (Ri, х); g = aT, (х) + 8,^ (х). Приведем текст функции граничных условий несимметричного конвективно-радиационного теплообмена конечного полого цилиндрического тела:
function [ q, g, h, r ] = boundaryFileHolParabolic(p, e, u, time)
global alfa2 Tc2 es2 Tn2 sig umy alfal Tc1 es1 Tn1
N = 1; ne = size(e,2); q = zeros(N^2, ne); g = zeros(N, ne);
h = zeros(N^2, 2*ne); r = zeros(N, 2*ne);
xy1 = p(:,e(1,:)); xy2 = p(:,e(2,:));
xyMidEdge = (xy1+xy2)/2; y = xyMidEdge(2, i);
T1 = umy(e(1,:)); T2 = umy(e(2,:)); TMidEdge = (T1+T2)/2; for i=1:ne switch(e(5,i)) case 1
q(1,i) = (alfa1+es1*sig*TMidEdge*3)*y;
g(1,i) = (alfa1 *Tc1+es1 *sig*Tn1^4)*y; case 3
q(1,i) = (alfa2+es2 *sig *TMidEdge^3)*y; g(1,i) = (alfa2 *Tc2+es2 *sig *Tn2n4)*y; end end
Результаты расчета двухмерного температурного поля в полом конечном цилиндрическом теле можно наглядно визуализировать с помощью функции Matlab pdeplot. Температурное поле в момент времени х = 2 мин. иллюстрирует рисунок 3.
Температурное поле в конечный момент времени (120 с)
О-1-1-1-1-1-1
-□ 03 -D 02 -0.01 0 0 01 0.02 0.03
Рис. 3. Температурное поле в цилиндрическом теле в конечный момент расчета
Fig. 3. A temperature field in a cylindrical body at the final moment of calculation
На рисунке видно отличное совпадение температурного распределения с результатами, полученными конечно-разностными методами и с помощью функции pdepe. Разработанные методы в перспективе позволяют перейти к решению задач моделирования температурного поля в сортовых стеклоизделиях цилиндрической формы (стакан, бутылка, колба), имеющих, кроме боковой, еще и донную поверхность. Предложенные модели температурного поля в дальнейшем можно использовать для идентификации условий теплообмена и оптимизации режимов термообработки стеклоизде-лий цилиндрической формы.
Литература
1. Гулоян Ю.А. Физико-химические основы технологии стекла. Владимир: Транзит-Икс, 2008. 736 с.
2. Мазурин О.В., Белоусов Ю.Л. Отжиг и закалка стекла. М.: Изд-во МИСИ и БТИСМ, 1984. 114 с.
3. Рубанов В.Г., Величко Д.В., Луценко О.В. Математическая модель динамики температурного поля
стеклоизделий сложной конфигурации при их отжиге // Стекло и керамика. 2018. № 5. С. 3-8.
4. Формалеев В.Ф. Теплопроводность анизотропных тел. Ч. 1. Аналитические методы решения задач. М.: Физматлит, 2014. 349 с.
5. Мазурин О.В., Лалыкин Н.В. Математическая модель процесса отжига листового стекла // Стекло и керамика. 1984. № 1. С. 13-15.
6. Марголис Б.И. Программа идентификации условий теплообмена для изделий плоской формы // Программные продукты и системы. 2017. № 1. С. 148-151. БО!: 10.15827/0236-235Х.117.148-151.
7. Дьяконов В.П. MATLAB: Полный самоучитель. М.: ДМК Пресс, 2014. 768 с.
8. Сирота А.А. Методы и алгоритмы анализа данных и их моделирование в MATLAB. СПб: БХВ-Петербург, 2016. 384 с.
9. Марголис Б.И. Программы моделирования температурных полей в изделиях плоской формы // Программные продукты и системы. 2016. № 2. С. 124-127. DOI: 10.15827/0236-235X. 114.124-127.
10. Крайнов А.Ю., Миньков Л.Л. Численные методы решения задач тепло- и массопереноса. Томск: STT, 2016. 92 с.
Software & Systems
DOI: 10.15827/0236-235X.126.313-317
Received 01.10.18 2019, vol. 32, no. 2, pp. 313-317
Programs for simulating temperature fields in cylindrical products
B.I. Margolis l, Dr.Sc. (Engineering), Head of Chair, [email protected] Mansoor Gubran Ali l, Postgraduate Student, [email protected]
1 Tver State Technical University, Department "Automation of technological processes", Tver, 170026, Russian Federation
Abstract. The paper describes a temperature field model in a cylindrical product during an asymmetric convec-tive-radiative heat transfer of surfaces with the environment and protecting surfaces (heating elements) of production equipment. The equations for calculating temperature distribution in the product using numerical finite-difference methods. There is a calculation example of convective-radiative cooling of a cylindrical body. It is stated that the described problem is possible to solve in Matlab.
Based on the standard MatLab function pdepe, there is a developed program that allows simulating a product temperature field according to specified on thermo-physical characteristics of the material (thermal conductivity, temperature conductivity) and parameters of the convective-radiative heat transfer (a convective heat transfer coefficient and given degrees of blackness). The paper describes the program development features related to assigning functions of differential equations, initial and boundary conditions. There are the main functions of the developed programs, as well as the temperature distribution calculation results.
The authors present a comparative analysis of the problem solution in Matlab with the results obtained using finite-difference relations. They also demonstrate their good agreement. A two-dimensional heat conduction problem in a cylindrical body is solved using a parabolic function, which is one of Matlab tools to solve parabolic type problems. The cylindrical domain has been triangulated using the initmesh function. There are some fragments of the program codes of boundary conditions functions and a visual image of a two-dimensional temperature distribution in a hollow cylinder. When solving a thermal two-dimensional cylindrical problem, the temperature gradients in height are almost absent. Therefore, there is an excellent coincidence of the temperature field with the results obtained by finite-difference methods and using the function pdepe.
The developed methods in the long term allow solving problems of modeling temperature field in high-quality cylindrical glass products (a glass, a bottle, a flask) that have a side face and a back surface.
Keywords: temperature field, cylindrical product, convective-radiative heat transfer, finite difference methods, modeling, initial and boundary conditions, differential equation, finite element method, triangulation, two-dimensional heat conduction problem.
References
1. Guloyan Yu.A. Physical and Chemical Basics of Glass Technology. Vladimir, Tranzit-Iks, 2008, 736 p.
2. Mazurin O.V., Belousov Yu.L. Annealing and Tempering Glass. Moscow, MISI i BTISM Publ., 1984, 114 p.
3. Rubanov V.G., Velichko D.V., Lucenko O.V. A mathematical model of the temperature field dynamics of complex configuration glassware during annealing. Glass and Ceramics. Moscow, Ladya Publ., 2018, no. 5, pp. 3-8 (in Russ.).
4. Formaleev V.F. Thermal Conductivity of Anisotropic Bodies. Part 1. Analytical Methods for Solving Problems. Moscow, Fizmatlit Publ., 2014, 349 p.
5. Mazurin O.V., Lalykin N.V. A mathematical model of sheet glass annealing process. Glass and Ceramics. Moscow, Stroyizdat Publ., 1984, no. 1, pp. 13-15 (in Russ.).
6. Margolis B.I. The program of identification of heat transfer conditions for flat products. Software & Systems. 2017, no. 1, pp. 148-151 (in Russ.). DOI: 10.15827/0236-235X.117.148-151.
7. Dyakonov V.P. MATLAB: Full Self-Teaching Guide. Moscow, DMK Press, 2014, 768 p.
8. Sirota A.A. Methods and Algorithms for Data Analysis and Modeling in MATLAB. St. Petersburg, BHV-Peterburg Publ., 2016, 384 p.
9. Margolis B.I. Programs for modeling temperature fields in flat-shaped products. Software & Systems. 2016, no. 2, pp. 124-127 (in Russ.). DOI: 10.15827/0236-235X.114.124-127.
10. Kraynov A.Yu., Minkov L.L. Numerical Methods for Solving Problems of Heat and Mass Transfer. Tomsk, STT Publ., 2016, 92 p.
