Программы моделирования температурных полей в изделиях плоской формы Текст научной статьи по специальности «Физика»

УДК 658.52 Дата подачи статьи: 25.11.15

DOI: 10.15827/0236-235X.114.124-127

ПРОГРАММЫ МОДЕЛИРОВАНИЯ ТЕМПЕРАТУРНЫХ ПОЛЕЙ В ИЗДЕЛИЯХ ПЛОСКОЙ ФОРМЫ

Б.И. Марголис, д.т.н.., зав. кафедрой, [email protected] (Тверской государственный технический университет, наб. Аф. Никитина, 22, г. Тверь, 1 70026, Россия)

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

Ключевые слова: температурное поле, изделие плоской формы, конвективно-радиационный теплообмен, конечно-разностные методы, моделирование, начальные и граничные условия, дифференциальное уравнение, идентификация условий теплообмена, оптимизация режимов термообработки.

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

дг (X, т) д2г (X, т)

— = а ~дьХХ~ , (1)

где ^, т) - значение температуры в координате по толщине x в момент времени т (°С); a - коэффициент температуропроводности материала изделия (м2/с) [3].

Большую роль для точного моделирования температурных распределений в изделиях играет корректный учет начальных и граничных условий, которые в большинстве стадий технологических процессов термообработки имеют сложный характер.

Так, начальные условия обычно неравномерны, граничные условия соответствуют несимметричному конвективно-радиационному теплообмену поверхностей изделия с окружающей средой и ограждающими поверхностями (нагревательными элементами) технологического оборудования [4, 5]. В этом случае наиболее удобны численные конечно-разностные методы, при использовании которых частные производные в уравнении (1) по времени и координате заменяются на их численные аналоги соответственно [и+1^) - ^(г)]/Дт и [^(г-1) -- 2^г) + ^(Г+1)]/(Д £ )2.

Здесь I = 0, т - номер временной точки расчета; г = О, N - номер точки по толщине 1\ Д £ = = £ /Ы £ - толщина слоя; Дт - шаг по времени [3].

Тогда формулы для расчета температурного распределения во внутренних точках плоского тела будут выглядеть следующим образом:

г Лг)=г (г)+Дт- а х

х (г -1) - 2г,. (г)+г,, (г+1) ]/ (м )2, (2)

г = 1, N -1.

При расчете конвективного теплообмена необходимо приравнять тепловые потоки на поверхности, например г=0, в соответствии с законами Фурье: q(0, т)=-Х-д^0, т)/дx и Ньютона: q(0, т) = = а,1^(0, т) - 41(т)] [6]. Тогда градиент д^0, т)/дx= = а,1[41(т) - ^0, т)]/Х, а вторая производная на поверхности г=0:

(3)

й (0) = Л (0, т)/йг2 = 2а, [/с, (т) - (0)]/(ХД£) и по аналогии для поверхности г=Ы: й (Ж) = т)/сТ = = 2а2 Ь2 (т) -I, (Ж)]/(ХМ),

где а1, а2 - коэффициенты конвективной теплоотдачи для нижней и верхней поверхностей изделия (Вт/(м2град)); /С1(т), 42(т) - температуры среды под и над пластиной (°С); X - коэффициент теплопроводности материала (Вт/(мград)).

При расчете радиационного теплообмена с учетом закона Стефана-Больцмана по аналогии с (3) и (4) вторые производные на поверхностях г=0 и г=Ы находятся по выражениям:

Й (0) = 28л1а[7» - Т, 4(0)]/(ХМ), (Ж) = 28,2а[Гв42(т) -Т4 (Ж)]/(Х^),

(4)

(5)

где s5i, zs2 - приведенные степени черноты нижнеи и верхнеИ ограждающих поверхностей технологического оборудования; ст = 5,6710-12 (Вт/(м2K4)) -постоянная Стефана-Больцмана; Thi(x), T#2(x) -температуры ограждающих поверхностей под и над изделием (K).

C учетом передачи тепла точкам поверхностей изделия от внутренних соседних точек за счет теплопроводности получим окончательные выражения для расчета температур поверхностей плоского тела в виде: 't.+l (0) = t. (0) + Дх- a х

х [2 • (t, (1) - tt(0)V(M)2 + Q^ (0) + Qr (0)],

t.+l (N) = t (N) + Дх^ a х

х[2 • (t, (N -1) - tt (N) )/(Де)2 + Qc, (N) + Qr

Для проверки модели температурного поля (2)-(5) был рассмотрен cлучай конвективно-радиационного охлаждения стеклянной пластины, приведенный в [3, 7]. При расчете использованы следующие значения параметров: a=0,2 (см2/мин.); Х=0,00838 (Вт/(смтрад.)); £ = 0,6 (см); е,1=е,2=0,91; N=6; М =0,1 (см); ai=a2=0,001256 (Вт/см2град.); fo(r)=707 (°С); 4i=42=27(°C); tHl=tH2=ll (°С); Дх=0,5 (с). Результаты расчета, приведенные в таблице 1, хорошо согласуются с источниками.

Таблица 1

Распределение температур в неограниченной пластине при симметричном конвективно-радиационном охлаждении

Table 1

Temperature distribution in a unconfined plate at symmetrical convective-radiant cooling

Время расчета х, с 0 6 10 20 30 60 120 œ

Поверхность t(0) 707,0 624,7 605,1 565,5 531,4 450,6 344,8 27,0

Центр t(3) 707,0 694,5 674,4 625,1 582,7 485,8 365,1 27,0

Решение поставленной задачи также может быть получено с помощью стандартной функции рёере среды программирования Matlab [8, 9]. Рассмотрим характерные особенности разработки такой программы.

Для использования рёере необходимо корректное задание трех функций: типа дифференциального уравнениярёеЛ/, начальных условийpdebeg и граничных условий pdebound. Решаемое дифференциальное уравнение вpdepe описывается следующим образом:

с ( х, х, t, dt/ dx )

dt_

дх

d

-s (

дх L х, х,t, dt/ дх )

xmf ( х, х, t, dt/ дх )]-

(6)

где m=0 - константа формы, равная 0, 1, 2 соответственно для пластины, цилиндра, шара; c(x, х, t, dt/dx)=1/a - коэффициент при dt/dx в формуле (1); fx, х, t, dt/dx)= dt/dx; s(x, x, t, dt/dx)=0 - отсутствующий член с источником излучения.

Граничные условия в pdepe описываются выражением для х=0 на нижней и х= t на верхней поверхностях плоского изделия:

p(x, х, t)+ q(x, х)- fx, х, t, dt/dx)=0, (7)

где q(0, х)= q( I, х)=1 - коэффициенты при градиентах 3/(0, х)/йг и dt{ I, т)/дх;

p(0,х, t) =«! [С) - t(0,х^/Х+е^Т,» -T4(0,х)]Д ;

p((, х, t) = -а2 [tc 2(х) -1((, х)]/Х-е, 2аТ Т» - T х)]/Х

- формулы для градиентов на поверхностях пластины при несимметричном конвективно-радиационном теплообмене.

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

function [c,f,s]=pdedif(x, tau, t, dtdx) global a c = 1/a; f = dtdx; s = 0;

function t0=pdebeg(x) t0= 707;

function [pl, ql,pr, qr]=pdebound(xl, tl,xr,tr, tau) global alfal alfa2 es1 es2 sig Kel lambd tc1=27; tc2=27; tn1=27; tn2=27;

pl=alfa1/lambd*(tc1-tl)+es1*sig/lambd*((to1+Kel)4-(tl+Kel)4); ql = 1;

pr=-alfa2/lambd*(tc2-tr)-es2*sig/lambd*((tn2+Kel)A4-(tr+Kel)A4); qr = 1;

Приведем текст основной программы с использованием функции pdepe: clear

close all clc

global a alfa1 alfa2 es1 es2 sig Kel lambd m = 0; sig=5.67e-12; Kel=273.15; es1=0.91; es2=es1;

= x

l=0.6; N=15; a=0.2; Ntau=241; lambd=0.00838;

taujin=2; alfa1=0.001256; alfa2=0.001256;

x = linspace(0,l,N);

tau = linspace(0,taufin,Ntau);

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

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

t = sol(:,:,1)

surf(x,tau,t)

Шle(strcat('Температурное поле для числа точек по координате N=',num2str(N))) xlabel(Координата x, см') ylabel(Время \tau, мин )

Полученная поверхность рассчитанного температурного распределения приведена на рисунке.

Результаты расчета температурного поля в пластине в среде Matlab приведены в таблице 2. Разница с результатами, полученными с помощью конечно -разностных соотношений (2)-(5), не превышает 1°C.

Таблица 2

Распределение температур в неограниченной пластине при симметричном конвективно-радиационном охлаждении, полученное в среде Matlab

Table 2

Temperature distribution in a unconfined plate at symmetrical convective-radiant cooling in Matlab

Время расчета z, с 0 6 10 20 30 60 120 œ

Поверхность t(0) 707,0 б24,7 б05,2 5б5,7 531,8 451,1 345,4 27,0

Центр t(3) 707,0 б93,8 б73,9 б25,0 582,9 48б,3 3б5,7 27,0

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

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

Литература

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

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

3. Марголис Б.И. Моделирование полей температур и напряжений в стеклоизделиях: монография. Тверь: Изд-во ТвГТУ, 2001. 100 с.

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

5. Марголис Б.И. Моделирование температурных полей в стеклоизделиях при их отжиге // Стекло и керамика. 2003. № 2. С. 3-5.

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

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

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

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

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

Температурное поле для числа точек по координате N=15

Температурное поле в стеклянной пластине A temperature field in a glass plate

DOI: 10.15827/0236-235X.114.124-127 Received 25.11.15

SIMULATOR PROGRAMS FOR TEMPERATURE FIELDS IN FLAT FORM PRODUCTS Margolis B.I., Dr.Sc. (Engineering), Head of Chair, [email protected] (Tver State Technical University, Nikitin Quay 22, Tver, 170026, Russian Federation) Abstract. The article considers a model of a temperature field in the flat shape product at asymmetric convective-radiant heat transfer of product surfaces to the environment and protecting surfaces (heating elements) of production equipment. The paper presents equations for calculating product temperature distribution using numerical finite-difference methods. There is a calculation example of convective-radiant cooling of a glass plate. The article shows that the calculation results are in a good

agreement with published sources. It also states possible solutions for the problem in Matlab. The standard function pdepe in MatLab has become a basis for development of a program that allows simulating a product temperature field according to specified thermal and physical characteristics of a material (thermal conductivity, temperature conductivity) and convective-radiant heat transfer parameters (a convective heat transfer coefficient and emissivity factors). The paper considers program development features related to setting functions such as differential equations, initial and boundary conditions. There are program codes of functions, of the main program and calculation results of the temperature distribution. A comparative analysis of problem solution in Matlab and the results obtained using the finite-difference relations demonstrates their good agreement. The paper shows the prospects of using standard Matlab functions to solve the identification problems of heat transfer conditions and optimization of product thermal treatment in various technological processes.

Keywords: temperature field, flat shape product, convective-radiant heat transfer, finite difference methods, modeling, initial and boundary conditions, differential equation, identification of heat transfer conditions, 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. Mazurin O.V., Belousov Yu.L. Otzhig i zakalka stekla [Glass Flashing and Toughening]. Study guide, Moscow, MISI i BTISM Publ., 1984, 114 p.

3. Margolis B.I. Modelirovaniepoley temperatur i napryazheny v stekloizdeliyakh [Modeling temperature fields and voltage in glass products]. Monograph. Tver, TSTU Publ., 2001, 100 p.

4. Mazurin O.V., Lalykin N.V. A mathematical model for the process of annealing flat glass. Glass and ceramics. Moscow, Stroyizdat Publ., 1984, vol. 41, iss. 1, pp. 9-13 (in Russ.).

5. Margolis B.I. Modeling Temperature Fields in Glass Articles in Annealing. Glass and ceramics. Moscow, Strojizdat Publ., 2003, vol. 60, iss. 1-2, pp. 33-35 (in Russ.).

6. Lykov A.V. Teoriya teploprovodnosti [The Theory of Heat Conduction]. Study guide, Moscow, Vysshaya shkola Publ., 1967, 599 p.

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

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. Glass and ceramics. Moscow, Strojizdat Publ., 2003, vol. 60, iss. 5, pp. 12-13 (in Russ.).

НИИ «Центрпрограммсистем» на Международном военно-техническом форуме

«АРМИЯ-2016»

С 6 по 11 сентября в г. Кубинка Московской области пройдет Международный военно-технический форум «АРМИЯ-2016».

НИИ «Центрпрограммсистем» представит на нем свои новейшие наиболее перспективные разработки в области информационных и автоматизированных систем, тренажеростроения и защиты информации.

Партнеры и заказчики могут ознакомиться с ними и обсудить перспективы сотрудничества в располагающей обстановке.