Программные продукты и системы

Математическое моделирование температурных полей в стеклоизделиях играет важнейшую роль при расчете процессов отжига стек- ла [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].

Для рассмотренных в [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

Для идентификации параметров теплообмена в среде 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

Из таблицы 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.).