ISSN 0236-235X (P)
ISSN 2311-2735 (E)

## Journal influence

Higher Attestation Commission (VAK) - К1 quartile
Russian Science Citation Index (RSCI)

## Next issue

4
Publication date:
16 September 2023

2023
2022

## Programs for simulating temperature fields in cylindrical products

Date of submission article: 10.10.2018
UDC: 658.52
The article was published in issue no. № 2, 2019 [ pp. 313-317 ]
Abstract:The paper describes a temperature field model in a cylindrical product during an asymmetric convective-radiative heat transfer of surfaces with the environment and protecting surfaces (heating elements) of production equip-ment. 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 co-efficient and given degrees of blackness). The paper describes the program development features related to assign-ing functions of differential equations, initial and boundary conditions. There are the main functions of the devel-oped 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 prob-lem 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 distri-bution in a hollow cylinder. When solving a thermal two-dimensional cylindrical problem, the temperature gradi-ents in height are almost absent. Therefore, there is an excellent coincidence of the temperature field with the re-sults 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.
Аннотация:В статье рассмотрена модель температурного поля в изделии цилиндрической формы при несимметричном конвективно-радиационном теплообмене поверхностей изделия с окружающей средой и ограждающими поверхностями (нагревательными элементами) технологического оборудования. Получены соотношения для расчета температурного распределения в изделии с использованием численных конечно-разностных методов. Приведен пример расчета конвективно-радиационного охлаждения цилиндрического тела. Сформулирована возможность решения по-ставленной задачи в среде программирования Matlab. На основе стандартной функции pdepe в среде MatLab разработана программа, позволяющая по заданным теплофизическим характеристикам материала (коэффициентам теплопроводности, температуропроводности) и параметрам конвективно-радиационного теплообмена (коэффициентам конвективной теплоотдачи и приведенным степеням черноты) моделировать температурное поле в изделии. Рассмотрены особенности разработки программы, связанные с заданием функций типа дифференциального уравнения, начальных и граничных условий. Приведены основные функции разработанных программ и результаты расчета температурного распределения. Сделан сравнительный анализ решения задачи в среде Matlab с результатами, полученными с помощью конечно-разностных соотношений, продемонстрировано их хорошее совпадение. Рас-смотрено решение двухмерной задачи теплопроводности в цилиндрическом теле с помощью функции parabolic, являющейся одним из инструментов Matlab для решения задач параболического типа. Произведена триангуляция цилиндрической области с использованием функции in-itmesh. Приведены фрагменты программных кодов функций граничных условий и визуальное изображение двухмерного температурного распределения в полом конечном цилиндре. При решении тепловой двухмерной цилиндрической задачи градиенты температур по высоте практически отсутствуют, в связи с чем наблюдается отличное совпадение температурного поля с результатами, полученными конечно-разностными методами и с помощью функции pdepe. Разработанные методы в перспективе позволяют перейти к решению задач моделирования температурного поля в сортовых стеклоизделиях цилиндрической формы (стакан, бутылка, колба), имеющих, кроме боковой, еще и донную поверхность.
 Authors: Margolis B.I. (borismargolis@yandex.ru) - Tver State Technical University, Tver, Russia, Ph.D, Mansoor Gubran Ali (gubran_ali@mail.ru) - Tver State Technical University, Department "Automation of technological processes", Tver, Russia Keywords: two-dimensional heat conduction problem, cylindrical product, convective-radiant heat transfer, finite difference methods, modeling, initial and boundary conditions, differential equation, finite element method, triangulation, a tempering Page views: 6101 PDF version article Full issue in PDF (6.72Mb)

Font size:       Font:

Математическое моделирование температурных полей в стеклоизделиях играет важнейшую роль при расчете процессов отжига стек- ла [1]. Геометрическая форма изделий может быть различной – плоской, цилиндрической. В таком случае одним из подходов при решении температурных задач является приведение сложного изделия к плоской форме за счет введения характерного или эффективного разме- ра [2]. Однако при этом точность расчетов ниже и не учитывается реальная геометрия тел цилиндрической формы. В данной работе рассматриваются подходы к моделированию тем- пературных полей в полых цилиндрических изделиях [3].

Для получения температурного поля необходимо решать дифференциальное уравнение теплопроводности в цилиндрических координатах:

(1)

где t(r, t) – значение температуры в координате по радиусу r в момент времени t (°C); a – коэффициент температуропроводности материала изделия (м2/с) [4].

Большую роль для точного моделирования температурных распределений в изделиях играет корректный учет начальных и граничных условий, которые в большинстве стадий технологических процессов термообработки имеют сложный характер. Так, начальные условия обычно неравномерны, граничные условия соответствуют несимметричному конвективно-радиационному теплообмену поверхностей изделия с окружающей средой и ограждающими поверхностями (нагревательными элементами) технологического оборудования [5, 6]. В этом случае наиболее удобно применять численные конечно-разностные методы, при использовании которых частные производные в уравнении (1) по времени, координате 2-го и 1-го порядков заменяются на их численные аналоги: [ti+1(j) – ti(j)]/Dt, [ti(j – 1) – 2ti(j) + + ti(j + 1)]/(Dr)2 и [ti(j + 1) – ti(j)]/Dr. Здесь  – номер временной точки расчета;  – номер точки по толщине; Dt – шаг по времени; Dr = (R2 – R1)/N – толщина слоя; R1, R2 – радиусы внутренней и внешней поверхностей полого цилиндра.

Тогда формулы для расчета температурного распределения во внутренних точках полого цилиндрического тела будут следующими:

(2)

При расчете конвективного теплообмена необходимо приравнять тепловые потоки на поверхности, например, r = R1, в соответствии с законами Фурье – q(R1, t) = –l×¶t(R1, t)/¶r и Ньютона –  q(R1, t) = a1[(t(R1, t) – tc1(t)] [4]. Градиент ¶t(R1, t)/¶r = a1[(tc1(t) – t(R1, t)]/l, производные 2-го порядка  на внутренней r = R1 и внешней r = R2 поверхностях:

(3)

где a1, a2 – коэффициенты конвективной теплоотдачи для внутренней и внешней поверхностей изделия (Вт/(м2град.)); tc1(t), tc2(t) – температура среды внутри и снаружи (°C); l – коэффициент теплопроводности материала (Вт/(м·град.)).

При расчете радиационного теплообмена с учетом закона Стефана–Больцмана по аналогии с (3) вторые производные на поверхностях r = R1 и r = R2 находятся по выражениям:

(4)

где es1, es2 – приведенные степени черноты верхней и нижней ограждающих поверхностей технологического оборудования; s = 5,67·10-12 (Вт/(м2K4)) – постоянная Стефана–Больцмана; Tн1(t), Tн2(t) – температуры ограждающих поверхностей над и под изделием (K).

C учетом передачи тепла точкам поверхностей изделия от внутренних соседних точек за счет теплопроводности получим окончательные выражения для расчета температур поверхностей цилиндрического тела в виде:

(5)

Для проверки модели температурного поля (2)–(5) был произведен расчет конвективно-радиационного охлаждения полого бесконечного цилиндра при следующих значениях параметров: a = 0,2 (см2/мин.); l = 0,00838 (Вт/(см∙град.)); R1 = 0,3 (см); R2 = 0,6 (см); εs1 = εs2 = 0,91; N = 14; a1 = a2 = 0,001256 (Вт/см2∙град.); t0(j) = 707 (°C); tс1 = tс2 = 27 (°C); tн1= tн2 = 27 (°C); Dt = 0,03 (с). Результаты расчета температуры в цилиндрическом теле приведены в таблице 1.

Таблица 1

Распределение температуры в точках по толщине неограниченного полого цилиндра

Table 1

Temperature distribution at points along the thickness of an unlimited hollow cylinder

 Температура Время расчета t, с 0 6 10 20 30 60 120 ¥ Внутренняя t(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 среды программирования Matlab [7, 8]. Характерные особенности разработки такой программы приведены в работе [9]. Функции типа дифференциального уравнения pdedif, начальных условий pdebeg и граничных условий pdebound не изменяются. Для учета цилиндрической формы тела необходимо задать константу формы m = 1.

Результаты расчета температурного поля в неограниченном полом цилиндре с использованием стандартных функций Matlab приведены в таблице 2 и на рисунке 1. Разница с результатами, полученными с помощью конечно-разностных соотношений (2)–(5), не превышает 1 °C.

Для повышения точности мо- делирования процесса отжига необходимо учитывать конечные размеры стеклоизделий цилиндрической формы, в связи с чем возникает задача перехода от одномерного температурного поля по толщине изделия к двухмерному распределению, учитывающему градиенты температуры по высоте изделия. В пакете Matlab имеется набор инструментов для решения дифференциального уравнения теплопроводности в пространственных цилиндрических координатах, в частности, функция parabolic [7, 8]. Функция parabolic – решатель многомерных краевых задач методом конечных элементов [10] параболического дифференциального уравнения в частных производных:

которое для случая осесимметричной температурной задачи в цилиндрической системе координат запишется следующим образом:

(6)

где z – координата по высоте изделия; d = l/a, c = l – коэффициенты при производных; 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 R1 R2 R2]') – описание геометрической области цилиндрического тела в прямоугольной системе координат; hmax – максимальный размер граничного элемента. Функция initmesh формирует выходные параметры: p – массив горизонтальных (по высоте) и вертикальных (по радиусу) координат узлов конечно-элементной сетки; e – матрица граничных элементов, содержащая номера начальных и конечных узлов этих элементов и номера граничных сегментов, которым они принадлежат; tr – матрица треугольных конечных элементов, содержащая номера входящих в них вершин (узлов).

Визуальное изображение разбиения цилиндрической области на треугольные элементы, полученное с помощью функции pdeplot(p, e, t), показано на рисунке 2.

Обращение к решателю parabolic имеет синтаксис: uv = parabolic(u0, tauv, b, p, e, tr, c, a1, f, d), где uv – температурное поле в моменты времени tauv, для которых вычисляется решение уравнения (6); u0 – узловое распределение температуры в начальный момент времени t = 0; b = @boundaryFileHolParabolic – имя функции, вычисляющей матрицы описания граничных условий.

Граничные условия, при которых осуществляется решение задач параболического типа, могут быть 1-го (Дирихле) и 3-го (Неймана) рода и описываются уравнениями hu = r и  соответственно. Для рассматриваемой задачи конвективно-радиационного теплообмена граничные условия 1-го рода отсутствуют, поэтому коэффициенты h = r = 0.

Граничные условия 3-го рода Неймана для температурной задачи:

–l×¶t(r, t)/¶r + qt(r, t) = g.                       (7)

В то же время для конвективно-радиационного теплообмена на поверхностях цилиндрического тела (например, внутренней) по законам Ньютона и Стефана–Больцмана [3, 4]:

(8)

В связи с этим с учетом перехода к абсолютным температурам коэффициенты q, g в уравнении (7) для внутренней цилиндрической поверхности записываются в виде   Приведем текст функции граничных условий несимметричного конвективно-радиационного теплообмена конечного полого цилиндрического тела:

function [ q, g, h, r ] = boundaryFileHolParabolic( p, e, u, time )

global alfa2 Tc2 es2 Tn2 sig umy alfa1 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*Tn2^4)*y;

end

end

Результаты расчета двухмерного температурного поля в полом конечном цилиндрическом теле можно наглядно визуализировать с помощью функции Matlab pdeplot. Температурное поле в момент времени t = 2 мин. иллюстрирует рисунок 3.

На рисунке видно отличное совпадение температурного распределения с результатами, полученными конечно-разностными методами и с помощью функции 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. DOI: 10.15827/0236-235X.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 с.

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.

 Permanent link:http://swsys.ru/index.php?page=article&id=4598&lang=&lang=en&like=1 Print versionFull issue in PDF (6.72Mb) The article was published in issue no. № 2, 2019 [ pp. 313-317 ]

Perhaps, you might be interested in the following articles of similar topics: