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

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

Для моделирования температурного поля используется стандартная функция pdepe Matlab [8, 9]. Функция типа дифференциального уравнения pdedif в силу исследования той же геометрической формы останется неизменной [6], тогда как функции начальных pdebeg и граничных условий pdebound изменятся следующим образом:

 function [c,f,s]=pdedif(x,tau,t,dtdx)

 global a

 c = 1/a;

 f = dtdx;

 s = 0;

 function t0=pdebeg(xv)

 global tbeg x

 [xm,num]=min(abs(xv-x));

 t0=tbeg(num);% tbeg;

 function [pl,ql,pr,qr]=pdebound(xl,tl,xr,tr,tau)

 global alfa1 alfa2 es1 es2 sig Kel lambd a1 a2 a3 a4 b1 b2 b3 b4 rs1 rs2

 tc1=a1+b1*tau;

 tc2=a2+b2*tau;

 tn1=a3+b3*tau;

 tn2=a4+b4*tau;

 pl = alfa1/lambd*(tc1-tl)+es1*sig/lambd*(rs1*(tn1+ +Kel)^4-(tl+Kel)^4);

 ql = 1;

 pr = -alfa2/lambd*(tc2-tr)-es2*sig/lambd*(rs2*(tn2+ +Kel)^4-(tr+Kel)^4);

 qr = 1;

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

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

 function y=Iden(par,m,tau)

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

 alfa1=par(1); alfa2=par(2); rs1=par(3); rs2=par(4);

 options=odeset('RelTol',1e-4);

 sol = pdepe(m,@pdedif,@pdebeg,@pdebound,x,tau, options);

 t = sol(:,:,1);

 tfin=t(end,:);

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

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

Таблица 1

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

Table 1

Initial data by annealing furnace zones

Зона отжига	Координата печи, м	Температура среды a1= a2, °C	Температура печи a3= a4, °C	Температура поверхности tl_iden, °C
1	1,8	479	540	477
2	3,6	510	560	496
3	5,7	518	560	513
4	7,5	537	550	531
5	10,5	539	550	532
6	13,5	541	550	533
7	15,9	532	540	526
8	18,2	507	530	500
9	20,7	480	530	473
10	24,3	453	520	446
11	26,2	439	439	430
12	29,9	423	423	413
13	31,6	388	388	371
14	36,4	366	366	352
15	38,8	346	346	329
16	43,3	325	325	306
17	48,4	301	301	284
18	50,8	285	285	260

Скорость движения изделия в печи составляет V=1,2 (м/мин.), поэтому по данным таблицы можно рассчитать время нахождения в зоне отжига и скорости изменения температур среды и печи на каждом этапе температурно-временного режима. Для упрощения расчетов рассмотрена симметричная задача конвективно-радиационного теплообмена, то есть a1=a2 и es1=es2. Данные таблицы 1 в среде Matlab удобно читать из файла Excel. При расчете использованы cледующие значения параметров: a=0,25386 см2/мин.; l=0,008856 Вт/см×град.; N=6; ℓ=0,6 см; εs1=εs2=0,85; α1=α2=0,001506 Вт/см2×град.; t0(r)= 444,0; 457,7; 468,0; 474,0; 475,1; 471,1; 461,7 °C.

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

clear

close all

clc

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

TTR = xlsread('ТВР период отжига');

Nzones=size(TTR,2)

 yst2=TTR(3,:); yst1=yst2; yst4=TTR(7,:);

 yst3=yst4;% среда и нагреватели слева, справа

 tl_exp=TTR(1,:);

coord=TTR(9,:); L=coord(Nzones); % координаты и длина (м)

 tau_finish=TTR(10,Nzones); V=L/tau_finish % время (мин); скорость ленты (м/мин)

tau_end=coord/V; % текущее время конца этапа ТВР (мин)

tau_zone(1)=tau_end(1); Vyst1(1)=0; Vyst2(1)=0; Vyst3(1)=0; Vyst4(1)=0;

 for i=2:Nzones

 tau_zone(i)=tau_end(i)-tau_end(i-1); % длительности этапов ТВР (мин)

 Vyst1(i)=(yst1(i)-yst1(i-1))/tau_zone(i); Vyst2(i)=(yst2(i)-yst2(i-1))/tau_zone(i);

 Vyst3(i)=(yst3(i)-yst3(i-1))/tau_zone(i); Vyst4(i)=(yst4(i)-yst4(i-1))/tau_zone(i);

 % скорости изменения температур среды и нагревателей на этапах (град/мин)

 end;

m = 0; sig=5.67e-12; Kel=273.15; es2=0.85; es1=es2; rs1=0.8; rs2=rs1;

alfa2=0.001506; alfa1=alfa2; a=0.6*0.4231; lambd=0.008856;

l=0.6; N=7; x = linspace(0,l,N); Ntau=20; N1=round(N/2);

tbeg=TTR(14,1:N);

par_prev=[alfa1; alfa2; rs1; rs2]; TolX=0.001;

A=[]; b=[]; nonlcon=[]; optset=[];

lb=[0.0015; 0.0015; 0.05; 0.05]; ub=[0.001506; 0.001506; 0.95; 0.95];

Aeq=[0 0 1 -1]; beq=[0]; par_all=[];

 for i=1:Nzones

 tau_fin=tau_zone(i);

 tau = linspace(0,tau_fin,Ntau);

 if i==1

 a1=yst1(i); a2=yst2(i); a3=yst3(i); a4=yst4(i);

 else

 a1=yst1(i-1); a2=yst2(i-1); a3=yst3(i-1); a4=yst4(i-1);

 end;

 b1=Vyst1(i);b2=Vyst2(i); b3=Vyst3(i);b4=Vyst4(i);

 

 par0=par_prev; tl_iden=tl_exp(i);

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

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

 par_all=[par_all; par'];

%         options=odeset('RelTol',1e-4);

%         sol = pdepe(m,@pdedif,@pdebeg,@pdebound, x,tau,options);

%         t = sol(:,:,1)

 t1=[t(:,1) t(:,N1) t(:,N)]

 if i==1

 surf(x,tau,t)

 title('Температурное поле в период отжига')

 xlabel('Координата x, см')

 ylabel('Время \tau, мин')

 grid on

 hold on

 else

 figure(1)

 surf(x,tau+tau_end(i-1),t)

 end;

 tbeg=t(end,:);

 par_prev=par;

 figure(2)

 plot(x,t(end,:),'r')

 title(strcat('Решение при \tau =', num2str(tau_end(i)),'(мин)'))

 xlabel('Координата x, см')

 ylabel('t(x,\tau)')

 grid on

 hold on

 end;

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

Из таблицы 2 видно, что отклонения рассчитанных температур верхней поверхности пластины от ее экспериментальных значений составляют десятые доли градуса, за исключением первой зоны, что можно объяснить неравномерностью условий теплообмена на входе в печь отжига. Полученные степени черноты rs2 представляют часть энергии радиационного теплообмена, попадающей с поверхности печи на изделие в каждой из зон. Их значения в диапазоне от 0,5 до 0,95 хорошо согласуются с физическими представлениями о радиационном теплообмене в печах отжига стеклоизделий.

Таблица 2

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

Table 2

Heat-exchange condition identification results

Зона отжига	Экспериментальная температура tl_iden, °C	Рассчитанная температура tfin(end), °C	Степень черноты поверхности печи rs2
1	477	479,2	0,7394
2	496	496,4	0,7722
3	513	512,8	0,8083
4	531	530,9	0,9207
5	532	531,8	0,9106
6	533	533,0	0,9153
7	526	526,0	0,9075
8	500	500,2	0,8101
9	473	472,5	0,7080
10	446	446,0	0,6423
11	430	430,0	0,7528
12	413	413,0	0,8820
13	371	371,0	0,5973
14	352	351,9	0,8225
15	329	329,0	0,7119
16	306	306,0	0,7414
17	284	284,1	0,7389
18	260	259,8	0,5361

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

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

Литература

1.     Рубанов В.Г. Автоматизация и управление объектами промышленности строительных материалов // Строительные материалы. 1996. № 2. С. 180–19.

2.     Михеев М.А., Михеева М.И. Основы теплопередачи. М.: Энергия, 1977. 344 с.

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

4.     Лыков А.В. Теория теплопроводности: учеб. пособие. М.: Высш. школа, 1967. 599 с.

5.     Gardon R. Calculation of temperature distributions in glass plates undergoing heat-treatment. J. Amer. Ceram. Soc., 1958, vol. 41, no. 6, pp. 200–209.

6.     Марголис Б.И. Программы моделирования температурных полей в изделиях плоской формы // Программные продукты и системы. 2016. № 2. С. 124–127.

7.     Мазурин О.В., Белоусов Ю.Л. Отжиг и закалка стекла: учеб. пособие. М.: Изд-во МИСИ и БТИСМ, 1984. 114 с.

8.     Лазарев Ю. Моделирование процессов и систем в MATLAB: учеб. курс. СПб: Питер, 2005. 512 с.

9.     Дьяконов В.П. MATLAB 7.*/ R2006/ R2007: самоучитель. М.: ДМК Пресс, 2008. 768 с.

10.   Марголис Б.И. Нахождение оптимального режима отжига стеклоизделий, обеспечивающего минимальные энергозатраты // Стекло и керамика. 2003. № 5. С. 12–13.