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

Bookmark

Next issue

4
Publication date:
16 September 2020
-->

Modeling programs and temperature field identification in sort glassware

Date of submission article: 2019-04-08
UDC: 658.52
The article was published in issue no. № 1, 2020 [ pp. 119-124 ]
Abstract:The paper proposes the structure and the main automated system modules for calculating the annealing modes of glass products. There is the problem of modeling temperature fields for cylindrical products with asymmetric con-vection-radiation heat exchange with the environment and enclosing surfaces of technological equipment. The au-thor formulated the possibility of solving the two-dimensional temperature problem in the Matlab programming environment. There was the parabolic function in the performance of the temperature field-modeling problem in cylindrical glass products. This is one of the Matlab tools for solving differential equations of parabolic type with complex boundary conditions. There are the programs based on the standard fmincon function in Matlab that allow to identify the coeffi-cients of artificial convective heat exchange between the air and the product surface in each of the zones by using the determinate thermal material characteristics and the temperature-time regime parameters in the annealing fur-nace. There is a parameter identification example of the convective-radiation heat transfer in the sort glassware (glasses) annealing furnace based on the temperature field-modeling program in Matlab. The author considered the program development features associated with the need to take into account the changing initial and boundary conditions at each temperature-time stage annealing of the product. There are the program codes of the main program, the modeling temperature field functions the CilinderHolParabolic1 and the boundary conditions boundaryFileHolParabolic1. The author studied the program results and demonstrated a good match of the convective heat transfer coefficients with their physically reasonable values in artificial convec-tive heat exchange in annealing furnaces of glassware. The attained condition identification results of the convec-tive-radiation heat exchange with the environment and enclosing surfaces we can use to optimize the heat treat-ment of cylindrical products.
Аннотация:В статье рассмотрены структура и основные модули автоматизированной системы расчета режимов отжига стеклоизделий. Приведена постановка задачи моделирования температурных полей для изделий цилиндрической формы при несимметричном конвективно-радиационном теплообмене с окружающей средой и ограждающими поверхностями технологического оборудования. Сформулирована возможность решения поставленной двухмерной температурной задачи в среде программирования Matlab. При решении задачи моделирования температурного поля в стеклоизделиях цилиндрической формы использована функция parabolic, являющаяся одним из инструментов Matlab для решения дифференциальных уравнений параболического типа со сложными граничными условиями. На основе стандартной функции fmincon в среде Matlab разработаны программы, позволяю-щие по заданным теплофизическим характеристикам материала и параметрам температурно-временного режима в печи отжига идентифицировать коэффициенты искусственного конвективного теплообмена между воздухом и поверхностью изделия в каждой из зон. Приведен пример идентификации параметров конвективно-радиационного теплообмена в печи отжига сортовых стеклоизделий (стаканов) на основе программы моделирования температурного поля в среде Matlab. Рассмотрены особенности разработки программы, связанные с необходимостью учета изменяющихся начальных и граничных условий на каждом из этапов температурно-временного режима отжига изделия. Приведены программные коды основной программы, функций моделирования температурного поля CilinderHolParabolic1 и граничных условий boundaryFileHolParabolic1. Произведен анализ результатов работы программы и продемонстрировано хорошее совпадение полученных коэффициентов конвективной теплоотдачи с их физически обоснованными значениями при искусственном конвективном теплообмене в печах отжига стеклоизделий. Полученные результаты идентификации условий конвективно-радиационного теплообмена с окружающей средой и ограждающими поверхностями можно использовать для оптимизации режимов термообработки изделий цилиндрической формы.
Authors: Margolis B.I. (borismargolis@yandex.ru) - Tver State Technical University, Tver, Russia, Ph.D, G.A. Mansour (gubran_ali@mail.ru) - Tver State Technical University, Tver, Russia
Keywords: triangulation, identification of heat transfer conditions, cylindrical product, convective-radiant heat transfer, initial and boundary conditions, modeling, a tempering, temperature-time mode, two-dimensional heat conduction problem
Page views: 1033
PDF version article
Full issue in PDF (4.91Mb)

Font size:       Font:

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

Рассмотрим структуру и основные модули автоматизированной системы расчета режимов отжига на основе анализа необходимых для ее работы данных. Для моделирования полей температур и напряжений в сортовых стеклоизде- лиях необходимо знание теплофизических (теплопроводность, температуропроводность, плотность, удельная теплоемкость) и механических (модуль Юнга, коэффициент Пуассона, фотоупругая постоянная) свойств стекла при любой температуре. Наиболее удобно определять значения этих свойств по известному химическому составу с использованием аддитивного метода, метода Аппена и хорошо известных линейных зависимостей [2]. В связи с этим в структуре системы должны присутствовать модули ввода химического состава (Состав) и расчета по нему вышеуказанных свойств (Свойства).

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

Для расчета полей напряжений в стеклоизделиях нужно иметь значения постоянных структурной релаксации. Необходимо определить постоянные, которые обеспечивают наименьшее отклонение рассчитанных по релаксационной модели значений характерного свойства стекла от его экспериментальных значений. Наиболее удобно в качестве таких экспериментальных данных использовать зависимость относительного удлинения стекла от температуры. Поиск релаксационных постоянных стекла осуществляется с помощью стандартного метода Matlab fmincon многомерного поиска минимума функции нескольких переменных [4] (модули Дилатометрия и Релаксация).

Тепловая история стеклоизделия начинается с процесса формования, после окончания которого стекло не сразу попадает в печь отжига, так как его доставка занимает некоторое время [1]. В связи с недостаточно высокой температурой на входе в печь возникает необходимость нагрева стеклоизделия до температуры отжига и расчета полей температур и напряжений также и для доотжигового периода. Различие условий теплообмена при нахождении стеклоизделия на участке от формующей машины до печи отжига и в ней делает необходимым наличие в автоматизированной системе расчета модулей Доотжиговый период и Период отжига. Составляющие систему модули приведены на рисунке 1.

В данной работе рассмотрена часть системы, связанная с моделированием тепловой истории стеклоизделий при их отжиге. Подходы к моделированию температурных полей в полых цилиндрических изделиях рассмотрены в [5]. Продемонстрированы возможности и перспективность использования метода parabolic математической среды программиро- вания Matlab [4, 6] для решения задач тепло- проводности в цилиндрических телах при сложных несимметричных условиях конвективно-радиационного теплообмена. Распространим рассмотренные в [5] методы на сортовые стеклоизделия, имеющие, кроме боковой, еще и донную поверхность. В качестве примера рассмотрим цилиндрический стакан.

Для решения цилиндрической задачи с использованием  parabolic необходимо произвести триангуляцию выбранной геометрической области с помощью функции initmesh. Для цилиндрического тела с полой боковой и сплошной донной поверхностями описание занимаемой геометрической области в прямоугольной системе координат выглядит следующим образом: reg1 = [3 4 Hb H H Hb R1 R1 R2 R2]'; reg2 = [3 4 0 Hb Hb 0 0 0 R2 R2]'; reg = [reg1 reg2]; g = decsg(reg), где R1, R2 – внутренний и внешний радиусы полой боковой поверхности; H, Hb – высота стакана и толщина его донной поверхности [3].

 

Рис. 2. Триангуляция цилиндрического стака-на

Fig. 2. Cylindric glass triangulation technique
Для рассмотренных в [5] геометрических размеров и Hb = 0,3 (см) визуальное изображение разбиения цилиндрической области стакана на треугольные элементы, полученное с помощью функции pdeplot(p, e, t), показано на рисунке 2.

На этом же рисунке показаны номера граничных поверхностей геометрической области цилиндрического тела: 1 – внутренняя боковая, 2 – торцевая, 4 – внешняя донная, 5 – внутренняя донная, 7 – наружная донная, 8 – наружная боковая. Обозначения 3 и 6 не являются физическими границами тела и соответствуют внутренним точкам донной поверхности на оси стакана и точкам на пересечении донной и боковой поверхностей. Для решения задачи моделирования температурного поля в цилиндрическом стакане по аналогии с работой [5] используем функцию Matlab parabolic  с измененным файлом граничных условий  b=boundaryFileHolParabolic1(p,e,umy,0). Приведем код функции граничных условий:

function [ q, g, h, r ] = boundaryFileHolParabolic1(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, 5}  % внутренние боковая и донная поверхности

            q(1,i) = (alfa1+es1*sig*TMidEdge^3)*y;

            g(1,i) = (alfa1*Tc1+es1*sig*Tn1^4)*y;

         case {2}% торцевая поверхность

           q(1,i) = (alfa2+es1*sig*TMidEdge(i)^3)*y(i);

           g(1,i) = (alfa2*Tc2+es1*sig*Tn1^4)*y(i);

         case {7, 8}% торцевая, наружные донная и боковая поверхности

           q(1,i) = (alfa2+0.5*(es1+es2)*sig*TMidEdge(i)^3)*y(i);

           g(1,i) = (alfa2*Tc2+0.5*sig*(es2*Tn2^4+es1*Tn1^4))*y(i);

         case 4% внешняя донная поверхность

            q(1,i) = (alfa2+es2*sig*TMidEdge^3)*y;

            g(1,i) = (alfa2*Tc2+es2*sig*Tn2^4)*y;

      end

   end

Результаты расчета двухмерного темпе- ратурного поля в стакане в момент времени t = 2 мин. для указанных в [5] условий приведены на рисунке 3, полученном с помощью функции Matlab pdeplot. По сравнению с приведенными в [5] результатами основное отличие заключается в наличии градиента температур по высоте изделия, наиболее сильно выраженное в верхней части боковой поверхности. Это связано с конвективно-радиационным охлаждением с торцевой поверхности стакана. Можно сделать вывод о том, что разработанные методы успешно используются для реше- ния задач моделирования температурных полей в сортовых стеклоизделиях цилиндрической формы.

При описании процесса отжига сортового стекла после создания математической модели температурного поля возникает  необходимость  проверки совпадения результатов расчетов с экспериментально измеренными температурами характерных точек изделия (для стакана это обычно температуры средних точек наружной боковой и внутренней донной поверхностей). Для решения этой задачи необходимо идентифицировать параметры конвективного (коэффициенты теплоотдачи для внутренней и внешней поверхностей изделия (a1, a2 (Вт/(м2град.))) и (или) радиационного (степени черноты верхней и нижней ограждающих поверхностей технологического оборудования es1, es2) теплообмена [7]. Выбор определяемых параметров зависит от типа оборудования в соответствии с преобладающими механизмами внешнего теплопереноса. Различают конвективные, радиационные и конвективно-радиационные печи [1]. В данной работе приведен пример идентификации параметров конвективно-радиационного теплообмена в печи отжига сортовых стеклоизделий (стаканов) на основе вышеприведенной программы моделирования температурного поля с файлом граничных условий boundaryFileHolParabolic1.

При расчете необходимо учитывать изменяющиеся начальные и граничные условия на каждом из этапов температурно-временного режима отжига изделия при его перемещении в следующую зону печи отжига [8, 9]. Приведем функцию для моделирования температурного поля в цилиндрических изделиях:

     function THol=CilinderHolParabolic1(tau, c, d, f, apar, p, e, tmesh, b);

     global lambd a sig Kel alfa1 alfa2 es1 es2 a1 b1 a2 b2 d3 a3 b3 a4 b4 Tc1 Tc2 Tn1 Tn2 umy

     tfin=tau(end);

     dt=tau(2)-tau(1); Nt=round(tfin/dt);

     tv=0:dt:dt;

     u0=d3; % initial temperature

     uKel=d3.*ones(size(p,2),1); umy=uKel;

     Tc1=a1; Tc2=a2; Tn1=a3; Tn2=a4;

     [ q, g, h, r ] = boundaryFileHolParabolic1( p, e, umy, 0);

        for i=1:Nt

           ti=i*dt;

           Tc1=a1+b1*ti;

           Tc2=a2+b2*ti;

           Tn1=a3+b3*ti;

           Tn2=a4+b4*ti;

           uv = parabolic1(u0, tv, b, p,e,tmesh,c,apar,f, d);

           u0=uv(:,2);

           uKel=[uKel u0];

           umy=u0;

           [ q, g, h, r ] = boundaryFileHolParabolic1( p, e, umy, 0);

        end;

     THol=uKel-Kel;

Величины a1, a2, a3, a4 (°C) и b1, b2, b3, b4 (°C/мин.) – начальные температуры и скорости изменения температур среды внутри и снаружи изделия и печи  под и над транспортирующей лентой для каждой из зон отжига. Параметры обращения к решателю parabolic1 среды Matlab подробно рассмотрены в работе [5]. Температура в программе выражена в градусах Кельвина для удобства записи радиационных тепловых потоков в файле граничных условий.

Вспомогательной при идентификации условий теплообмена является функция Iden1, определяющая отклонение рассчитанной по модели температуры середины наружной боковой поверхности стакана от экспериментально полученной tl_iden:

     function y=Iden1(par,m,tau)

     global a1 b1 a2 b2 a es1 es2 sig Kel lambd alfa1 alfa2 a3 b3 a4 b4 tbeg x tl_iden t

     global c d f apar p e tmesh b hiden

     alfa1=par(1); alfa2=par(2); es1=par(3); es2=par(4);

    t=CilinderHolParabolic1(tau, c, d, f, apar, p, e, tmesh, b);

    tfin=t(:,end);

     y=abs(tfin(hiden)-tl_iden);

Исходные данные для решения задачи идентификации, описывающие условия в цехе отжига сортового стекла на Спировском стеклозаводе «Индустрия», приведены в таблице 1.

Скорость движения изделия в печи составляет V = 0,27 (м/мин.), поэтому по данным таблицы 1 можно рассчитать время нахождения в зоне отжига и скорости изменения температур среды и печи на каждом этапе температурно-временного режима. Рассматриваются равные степени черноты es1 = es2 и естественный конвективный теплообмен для внутренних и искусственный конвективный теплообмен для наружных поверхностей стеклоизделия. Данные таблицы 1 в среде Matlab удобно читать из файла Excel. При расчете использованы значения параметров: a = 0,4351·10-6 (м2/с), l = 0,9024 (Вт/(м∙град.)).

Таблица 1

Исходные данные по зонам печи отжига

Table 1

Feed data zoned at annealing furnaces

Зона отжига

Координата печи, м

Температура сре­ды a1= a2, °C

Температура печи a3= a4, °C

Температура поверхности tl_iden, °C

1

2,4

545

540

545

2

4,8

525

520

526

3

7,2

490

480

491

4

9,6

455

425

456

5

12,0

345

325

346

6

14,4

300

280

302

7

16,8

250

230

251

8

19,2

200

180

201

9

21,6

150

130

151

10

24,0

100

80

101

Для идентификации параметров теплообмена в среде Matlab необходимо использовать функцию fmincon, позволяющую эффективно решать задачи нелинейной оптимизации функции нескольких переменных с ограничениями типа неравенств и равенств. Приведем фрагмент основной программы с использованием функции fmincon:

global Nzones Ntau yst es1 es2 alfa1 alfa2 tl_exp tau_end tau_zone tbeg lambd tl_iden

global Vyst Kel d3 a f apar p e tmesh b H hin hmid hiden a1 b1 a2 b2 a3 b3 a4 b4 m c d t

      Tbeg=tbeg;

     [klam,kroc]=TempDepend(tbeg, temp, alfa); Bio_nature=3*4.1868*(R2-R1)/(lambd*klam);

     par_all=[]; TolX=0.02; nonlcon=[]; optset=[];

    par_prev=[alfa1; alfa2; es1; es2];  Aeq=[]; beq=[]; A=[]; bneq=[];

        for i=1:Nzones

          tau_fin=tau_zone(i)*60;

          tau = linspace(0,tau_fin,Ntau);

          a1=yst(1,i)-Vyst(1,i)*tau_zone(i)+Kel; a2=yst(2,i)- -Vyst(2,i)*tau_zone(i)+Kel;

          a3=yst(3,i)-Vyst(3,i)*tau_zone(i)+Kel; a4=yst(4,i)- -Vyst(4,i)*tau_zone(i)+Kel;

          b1=Vyst(1,i)/60; b2=Vyst(2,i)/60; b3=Vyst(3,i)/60; b4=Vyst(4,i)/60; d3=Tbeg+Kel;

          [klam,kroc]=TempDepend(Tbeg, temp, alfa);

          c = sprintf('%g*y', klam*lambd);

          d = sprintf('%g*y', kroc*lambd/a);

          alfa1=Bio_nature*klam*lambd/(R2-R1); alfa2=10*alfa1;

          par0=[alfa1; alfa2; par_prev(3); par_prev(4)];

          tl_iden=tl_exp(i); lb=[alfa1/keps; 6; par0(3)/keps; par0(4)/keps];

          ub=[alfa1*keps; 1000; par0(3)*keps; par0(4)*keps];

         optset=optimset('Display','iter','TolX',TolX)

         [par,fval]=fmincon(@Iden1,par0,A,bneq,Aeq,beq,lb,ub,nonlcon,optset,m,tau);

          par_all=[par_all; par'];

        end;

     Tbeg=t(:,end); MeshOut(i, p, e, tmesh, H, Tbeg, Sfig, tau_end, iden_yes2);

     SurfOut(i, hin, hmid, hiden, t, x1, tau, tau_end,' период отжига');

Рассчитанное температурное распределение показано на рисунке 4, идентифицированные значения параметров теплообмена и температур поверхности изделия приведены в таблице 2.

Таблица 2

Результаты идентификации условий теплообмена

Table 2

The Identification results of the heat-exchange conditions

Зона отжига

Экспериментальная температура tl_iden, °C

Рассчитанная температура tfin(hiden), °C

Коэффициент теплоотдачи α2

1

545

545,0

103,3

2

526

526,2

103,4

3

491

492,2

102,3

4

456

457,4

100,2

5

346

346,0

969,1

6

302

304,0

91,1

7

251

251,0

471,6

8

201

201,0

479,5

9

151

151,0

468,5

10

101

101,0

465,3

Из таблицы 2 видно, что отклонения рассчитанных температур середины наружной бо- ковой поверхности стакана от ее экспериментальных значений составляют десятые доли градуса. Полученные коэффициенты конвективной теплоотдачи a2 согласуются с физическими представлениями об искусственном конвективном теплообмене в печах отжига стеклоизделий.

Полученные результаты идентификации условий конвективно-радиационного теплообмена с окружающей средой и ограждающими поверхностями можно использовать для оптимизации режимов термообработки изделий цилиндрической формы с целью сокращения длительности технологического процесса или энергетических затрат на него [10].

Наличие в Matlab наряду с parabolic функции fmincon, решающей задачи оптимизации функции нескольких переменных с ограничениями, позволяет эффективно использовать ее для оптимизации режимов термообработки стеклоизделий цилиндрической формы.

Литература

1.    Гулоян Ю.А. Физико-химические основы технологии стекла. Владимир: Транзит-Икс, 2008. 736 с.

2.    Матвеев М.А., Матвеев Г.М., Френкель Б.Н. Расчеты по химии и технологии стекла. М., 1972. 240 с.

3.    Мазурин О.В., Белоусов Ю.Л. Отжиг и закалка стекла. М., 1984. 114 с.

4.    Дьяконов В.П. MATLAB: Полный самоучитель. М.: ДМК Пресс, 2014. 768 с.

5.    Марголис Б.И., Мансур Г.А. Программы моделирования температурных полей в изделиях цилиндрической формы // Программные продукты и системы. 2019. № 2. С. 313–317. DOI: 10.15827/0236-235X.126.313-317.

6.    Сирота А.А. Методы и алгоритмы анализа данных и их моделирование в MATLAB. СПб: БХВ-Петербург, 2016. 384 с.

7.    Марголис Б.И. Программа идентификации условий теплообмена для изделий плоской формы // Программные продукты и системы. 2017. № 1. С. 148–151. DOI: 10.15827/0236-235X.117.148-151.

8.    Формалеев В.Ф. Теплопроводность анизотропных тел. Ч. 1: Аналитические методы решения задач. М.: Физматлит, 2014. 349 с.

9.      Крайнов А.Ю., Миньков Л.Л. Численные методы решения задач тепло- и массопереноса. Томск, 2016. 92 с.

10. Мазурин О.В., Лалыкин Н.В. Математическая модель процесса отжига листового стекла // Стекло и керамика. 1984. № 1. С. 13–15.

References

  1. Guloyan Yu.A. Physical-Chemical Bases of Glass Technology. Vladimir, 2008, 736 p. (in Russ.).
  2. Matveev M.A., Matveev G.M., Frenkel B.N. Calculations on Chemistry and Glass Technology. Moscow, 1972, 240 p. (in Russ.).
  3. Mazurin O.V., Belousov Ju.L. Annealing and Tempering of Glass. Moscow, 1984, 114 p. (in Russ.).
  4. Dyakonov V.P. MATLAB: Complete Tutorial. Moscow, DMK Press, 2014, 768 p. (in Russ.).
  5. Margolis B.I., Mansoor G.A. Programs for simulating temperature fields in cylindrical products. Software & Systems. 2019, no. 2, pp. 313–317. DOI: 10.15827/0236-235X.126.313-317 (in Russ.).
  6. Sirota A.A. Methods and Data Analysis Algorithms and Their Simulation in MATLAB. Saint-Petersburg, BHV-Petersburg Publ., 2016, 384 p. (in Russ.).
  7. Margolis B.I. Program for heat conditions identification in flat products. Software & Systems. 2017, no. 1,
    pp. 148–151. DOI: 10.15827/0236-235X.117.148-151 (in Russ.).
  8. Formaleev V.F. The thermal conductivity of anisotropic solids. Pt. 1. Analytical Approaches of Solving Problems. Moscow, Fizmatlit, 2014, 349 p. (in Russ.).
  9. Kraynov A.Yu., Minkov L.L. A Numerical Method for the Calculation of Heat and Mass Transfer. Tomsk, 2016, 92 p. (in Russ.).
  10. 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.).

Permanent link:
http://swsys.ru/index.php?page=article&id=4685&lang=&lang=en&like=1
Print version
Full issue in PDF (4.91Mb)
The article was published in issue no. № 1, 2020 [ pp. 119-124 ]

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