Software & Systems

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

Скорость движения изделия в печи составляет 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

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