CC BY
25
УДК 658.52 Дата подачи статьи: 02.12.16
DOI: 10.15827/0236-235X^30.1.148-151 2017. Т. 30. № 1. С. 148-151
ПРОГРАММА ИДЕНТИФИКАЦИИ УСЛОВИЙ ТЕПЛООБМЕНА ДЛЯ ИЗДЕЛИЙ ПЛОСКОЙ ФОРМЫ
Б.И. Марголис, д.т.н.., зав. кафедрой, [email protected] (Тверской государственный технический университет, наб. Аф. Никитина, 22, г. Тверь, 1 70026, Россия)
Рассмотрена постановка задачи идентификации условий теплообмена для изделия плоской формы при несимметричном конвективно-радиационном теплообмене поверхностей изделия с окружающей средой и ограждающими поверхностями (нагревательными элементами) технологического оборудования. Сформулирована возможность решения поставленной задачи в среде программирования Matlab.
На основе стандартной функции/ттеоп в среде Matlab разработана программа, позволяющая по заданным тепло-физическим характеристикам материала (коэффициентам теплопроводности, температуропроводности), параметрам конвективно-радиационного теплообмена (коэффициентам конвективной теплоотдачи и приведенным степеням черноты) и параметрам температурно-временного режима в печи отжига идентифицировать часть энергии радиационного теплообмена, попадающей с поверхности печи на изделие в каждой из зон.
Приведен пример идентификации параметров радиационного теплообмена в печи отжига листового прокатного стекла на основе программы моделирования температурного поля в среде Matlab. Рассмотрены особенности разработки программы, связанные с необходимостью учета изменяющихся начальных и граничных условий на каждом из этапов температурно-временного режима отжига изделия с помощью функций pdebeg и р<!еЪои^ стандартной функции pdepe Matlab. Приведены программные коды функций, основной программы и результаты расчета температур поверхности ленты и степеней черноты поверхности печи.
Произведен анализ результатов работы программы, и продемонстрировано хорошее совпадение полученных степеней черноты с их физически обоснованными значениями при радиационном теплообмене в печах отжига стеклоиз-делий.
Показана перспективность использования стандартных функций Matlab для решения задач оптимизации режимов термообработки изделий в различных технологических процессах.
Ключевые слова: идентификация условий теплообмена, изделие плоской формы, конвективно-радиационный теплообмен, степень черноты, начальные и граничные условия, моделирование, температурное поле, температурно-временной режим, оптимизация режимов термообработки.
Определение параметров конвективно-радиационного теплообмена играет важнейшую роль в большинстве технологических процессов, связанных с термообработкой изделий [1]. Например, при описании процесса отжига листового стекла после создания математической модели температурного поля [2-6] возникает необходимость проверки совпадения результатов расчетов с экспериментально измеренными температурами характерных точек изделия (для ленты это обычно температура верхней поверхности).
Для решения этой задачи необходимо идентифицировать параметры конвективного (коэффициенты теплоотдачи для нижней и верхней поверхностей изделия а.1, а,2 (Вт/(м2град))) и (или) радиационного (степени черноты нижней и верхней ограждающих поверхностей технологического оборудования е^, е^) теплообмена. Выбор определяемых параметров зависит от типа оборудования в соответствии с преобладающими механизмами внешнего теплопереноса. Различают конвективные, радиационные и конвективно-радиационные печи [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 alfal 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)A4-(tl+Kel)A4); ql = 1;
pr = -alfa2/lambd*(tc2-tr)-es2*sig/lambd*(rs2*(tn2+ +Kel)A4-(tr+Kel)A4); qr = 1;
Это связано с необходимостью учета изменяющихся начальных и граничных условий на каждом из этапов температурно-временного режима отжига изделия при его перемещении в следующую зону печи отжига. В приведенных выше функциях a1, a2, a3, a4(°C) и b1, b2, b3, b4(°С/мин.) - начальные температуры и скорости изменения температур среды и печи соответственно под и над лентой для каждой из зон отжига.
Вспомогательной при идентификации условий теплообмена является приведенная далее функция Iden, определяющая отклонение рассчитанной по модели температуры верхней поверхности ленты от экспериментально полученной tliden:
function y =Iden (par,m,tau)
global l al b1 a2 b2 a es1 es2 sig Kel lambd alfal alfa2 a3 b3 a4 b4 tbeg x tl iden t rsl 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); tfn=t(end,:); y=abs(fn (end)-tlid en);
Исходные данные для решения задачи идентификации, описывающие условия в цехе отжига листового прокатного стекла, приведены в таблице 1.
Таблица 1
Исходные данные по зонам печи отжига
Table 1
Initial data by annealing furnace zones
Скорость движения изделия в печи составляет К=1,2 (м/мин.), поэтому по данным таблицы можно рассчитать время нахождения в зоне отжига и скорости изменения температур среды и печи на каждом этапе температурно-временного режима. Для упрощения расчетов рассмотрена симметричная задача конвективно-радиационного теплообмена,
то есть ai=a,2 и 6*1=6*2. Данные таблицы 1 в среде Matlab удобно читать из файла Excel. При расчете использованы следующие значения параметров: а=0,25386 см2/мин.; Х=0,008856 Вт/смхград.; N=6; 1=0,6 см; е^1=е^2=0,85; ai=a2=0,001506 Вт/см2хград.; to(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 alfal alfa2 a3 b3 a4 b4 tbeg x tl iden t rsl rs2
TTR = xlsread('ТВР период отжига'); Nzones=size(TTR,2)
yst2=TTR(3,:); ystl =yst2; yst4=TTR(7,:); yst3=yst4;% среда и нагреватели слева, справа tl_exp=TTR(1,:);
coord=TTR(9,:); L=coord(Nzones); % координаты и длина (м)
tau_finish=TTR(10,Nzones); V=L/taufinish % время (мин); скорость ленты (м/мин)
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); ifi==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);
Зона Коор- Темпера- Темпера- Температура
от- дината тура среды тура печи поверхности
жига печи, м ai= ai, °C «з= «4, °C tl iden, °C
i 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
ii 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
l49
par0=par_prev; tl_iden=tl_exp(i); optset=optimset('Display', 'iter', 'TolX',TolX) [parfval]=fmin-con (@Jden,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)
шиСТемпературное поле в период отжига') xlabel('Координата x, см') ylabel('Время \tau, мин') grid on hold on else figure(1)
surf(x, tau+tauend (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 49б 49б,4 0,7722
3 513 512,8 0,8083
4 531 530,9 0,9207
5 532 531,8 0,910б
б 533 533,0 0,9153
7 52б 52б,0 0,9075
8 500 500,2 0,8101
9 473 472,5 0,7080
10 44б 44б,0 0,б423
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
1б 30б 30б,0 0,7414
17 284 284,1 0,7389
18 2б0 259,8 0,53б1
Представленные результаты идентификации условий конвективно-радиационного теплообмена с окружающей средой и ограждающими поверхностями можно использовать для оптимизации режимов термообработки изделий плоской формы с целью сокращения длительности технологического процесса или энергетических затрат на него [10].
Наличие в МайаЬ наряду с pdepe функции ттсоп, решающей задачи оптимизации функции нескольких переменных с ограничениями, позволяет эффективно использовать ее для оптимизации режима отжига стеклоизде-лий плоской формы.
Время t, мин
Координата х, см
Температурное поле в ленте в период отжига Temperature field in a band during annealing
Литература
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. N° 5. С. 12-13.
Software & Systems Received 02.12.16
DOI: 10.15827/0236-235X.030.1.148-151 2017, vol. 30, no. 1, pp. 148-151
PROGRAM FOR HEAT CONDITIONS IDENTIFICATION IN FLAT PRODUCTS B.I. Margolis1, Dr.Sc. (Engineering), Head of Chair, [email protected] 1 Tver State Technical University, Nikitin Quay 22, Tver, 170026, Russian Federation
Abstract. The article considers the problem statement on identification of heat transfer conditions for flat products with asymmetrical shapes convective-radiative heat transfer surfaces to the environment and enclosing surfaces (heating elements) of the process equipment. The paper formulates a possible solution of the problem in Matlab.
The program, which has been developed based on the standard fmincon function in MatLab, allows identifying the part of radiation heat transfer energy, which gets from a furnace surface to the product in each zone. For this purpose it uses predetermined thermo-physical characteristics of the material (thermal conductivity, thermal diffusivity), the parameters of convective-radiative heat transfer (coefficients of convective heat transfer and an emissivity factor) and the temperature and time parameters on an annealing furnace.
The paper presents an example of radiative heat transfer parameter identification in an annealing furnace for rolled glass sheet on the basis of the temperature field simulation program in Matlab. The authors consider the features of program development related to the need to take into account changing initial and boundary conditions at each stage of the temperature-time mode of an annealing product using pdebeg andpdebound of the standard functionspdepe in MatLab. There are software codes of functions and main program, as well as the results of calculating band surface temperatures and emissivity factors of a furnace surface.
There is the analysis of the results of the program. The paper demonstrates good agreement between the obtained emissivity factors and their physically reasonable values for radiative heat transfer in glass annealing furnaces. The paper shows the prospects of using standard Matlab functions to solve product heat treatment mode optimization problems in various technological processes.
Keywords: heat transfer conditions identification, flat shape product, convective-radiant heat transfer, emissivity factor, initial and boundary conditions, simulation, temperature field, temperature-time mode, heat treatment optimization.
References
1. Rubanov V.G. Automation and management of industry facilities of construction materials. Stroitelnye materialy [Construction Materials]. 1996, no. 2, pp. 18-19 (in Russ.).
2. Mikheev M.A., Mikheeva M.I. Osnovy teploperedachi [Heat-Transfer Principles]. Moscow, Energiya Publ., 1977, 344 p.
3. Mazurin O.V., Lalykin N.V. A mathematical model of a flat glass annealing process. Steklo i keramika [Glass and Ceramics]. Moscow, Stroyizdat Publ., 1984, vol. 41, iss. 1, pp. 9-13 (in Russ.).
4. Lykov A.V. Teoriya teploprovodnosti [Thermal Conductivity Theory]. Study guide, Moscow, Vysshaya shkola, 1967, 599 p.
5. Gardon R. Calculation of temperature distributions in glass plates undergoing heat-treatment. Jour. Amer. Ceram. Soc. 1958, vol. 41, no. 6, pp. 200-209.
6. Margolis B.I. Simulator programs for temperature fields in flat form products. Programmnyeprodukty i sistemy [Software & Systems]. 2016, no. 2, pp. 124-127 (in Russ.).
7. Mazurin O.V., Belousov Yu.L. Otzhig i zakalka stekla [Glass Flashing and Toughening]. Study guide, Moscow, MISI i BTISM Publ., 1984, 114 p.
8. Lazarev Yu. Modelirovanieprotsessov i sistem v MATLAB [Modeling processes and systems in MATLAB]. Training course, St. Petersburg, Piter Publ., 2005, 512 p.
9. Dyakonov V.P. MATLAB 7. */R2006/R2007: samouchitel [MATLAB 7.*/ R2006/ R2007: Teach Yourself]. Moscow, DMK Press, 2008, 768 p.
10. Margolis B.I. An Optimum Regime for Glass Annealing with Minimum Energy Consumption. Steklo i keramika [Glass and Ceramics]. Moscow, Stroyizdat Publ., 2003, vol. 60, iss. 5, pp. 12-13 (in Russ.).
