ISSN 0236-235X (P)
ISSN 2311-2735 (E)

Journal influence

Higher Attestation Commission (VAK) - К1 quartile
Russian Science Citation Index (RSCI)

Bookmark

Next issue

4
Publication date:
09 September 2024

Simulator programs for temperature fields in flat form products

Date of submission article: 25.11.2015
UDC: 658.52
The article was published in issue no. № 2, 2016 [ pp. 124-127 ]
Abstract:The article considers a model of a temperature field in the flat shape product at asymmetric convecti ve-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 convectiveradiant 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.
Аннотация:Рассмотрена модель температурного поля в изделии плоской формы при несимметричном конвективно-радиационном теплообмене поверхностей изделия с окружающей средой и ограждающими поверхностями (нагревательными элементами) технологического оборудования. Получены соотношения для расчета температурного распределения в изделии с использованием численных конечно-разностных методов. Рассмотрен пример расчета конвективно-радиационного охлаждения стеклянной пластины. Показано, что результаты расчета хорошо согласуются с известными литературными источниками. Сформулирована возможность решения поставленной задачи в среде программирования Matlab. На основе стандартной функции pdepe в среде Matlab разработана программа, позволяющая по заданным теплофизическим характеристикам материала (коэффициентам теплопроводности, температуропроводности) и параметрам конвективно-радиационного теплообмена (коэффициентам конвективной теплоотдачи и приведенным степеням черноты) моделировать температурное поле в изделии. Рассмотрены особенности разработки программы, связанные с заданием функций типа дифференциального уравнения, начальных и граничных условий. Приведены программные коды функций, основной программы и результаты расчета температурного распределения. Произведен сравнительный анализ решения задачи в среде Matlab с результатами, полученными с помощью конечно-разностных соотношений, и продемонстрировано их хорошее совпадение. Показана перспективность использования стандартных функций Matlab для решения задач идентификации условий теплообмена и оптимизации режимов термообработки изделий в различных технологических процессах.
Authors: Margolis B.I. (borismargolis@yandex.ru) - Tver State Technical University, Tver, Russia, Ph.D
Keywords: heat treatment optimization, identification of heat transfer conditions, differential equation, initial and boundary conditions, modeling, finite difference methods, convective-radiant heat transfer, flat shape product, a tempering
Page views: 12018
Print version
Full issue in PDF (7.11Mb)
Download the cover in PDF (0.37Мб)

Font size:       Font:

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

                                    (1)

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

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

Здесь  – номер временной точки расчета;  – номер точки по толщине  D= = /N – толщина слоя; Dt – шаг по времени [3].

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

                 (2)

При расчете конвективного теплообмена необходимо приравнять тепловые потоки на поверхности, например r=0, в соответствии с законами Фурье: q(0, t)=–l׶t(0, t)/¶x и Ньютона: q(0, t) = = a1[t(0, t) – tc1(t)] [6]. Тогда градиент ¶t(0, t)/¶x= = a1[tc1(t) – t(0, t)]/l, а вторая производная на поверхности r=0:

 

и по аналогии для поверхности r=N:

                             (3)

где a1, a2 – коэффициенты конвективной теплоотдачи для нижней и верхней поверхностей изделия (Вт/(м2град)); tc1(t), tc2(t) – температуры среды под и над пластиной (°C); l – коэффициент теплопроводности материала (Вт/(м·град)).

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

      (4)

где es1, es2 – приведенные степени черноты нижней и верхней ограждающих поверхностей технологического оборудования; s = 5,67·10-12 (Вт/(м2K4)) – постоянная Стефана–Больцмана; TH1(t), TH2(t) – температуры ограждающих поверхностей под и над изделием (K).

C учетом передачи тепла точкам поверхностей изделия от внутренних соседних точек за счет теплопроводности получим окончательные выражения для расчета температур поверхностей плоского тела в виде:

(5)

Для проверки модели температурного поля (2)–(5) был рассмотрен cлучай конвективно-радиационного охлаждения стеклянной пластины, приведенный в [3, 7]. При расчете использованы cледующие значения параметров: a=0,2 (см2/мин.); l=0,00838 (Вт/(см∙град.)); = 0,6 (см); εs1=εs2=0,91; N=6; =0,1 (см); a1=a2=0,001256 (Вт/см2град.); t0(r)=707 (°C); tс1=tс2=27(°C); tн1=tн2=27 (°C); Dt=0,5 (с). Результаты расчета, приведенные в таблице 1, хорошо согласуются с источниками.

Таблица 1

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

Table 1

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

Время расчета t, с

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

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

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

                      (6)

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

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

p(x, t, t)+ q(x, t)× f(x, t, t, ¶t/¶x)=0,                      (7)

где q(0, t)= q(, t)=1 – коэффициенты при градиентах ¶t(0, t)/¶x  и ¶t( t)/¶x;

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

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

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 alfa1 alfa2 es1 es2 sig Kel lambd

tc1=27; tc2=27;

tn1=27; tn2=27;

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

ql = 1;

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

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;

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

 tau_fin=2; alfa1=0.001256; alfa2=0.001256;

 x = linspace(0,l,N);

 tau = linspace(0,tau_fin,Ntau);

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

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

 t = sol(:,:,1)

 surf(x,tau,t)

 title(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

Время расчета t, с

0

6

10

20

30

60

120

¥

Поверхность t(0)

707,0

624,7

605,2

565,7

531,8

451,1

345,4

27,0

Центр t(3)

707,0

693,8

673,9

625,0

582,9

486,3

365,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.


Permanent link:
http://swsys.ru/index.php?page=article&id=4158&lang=&lang=en&like=1
Print version
Full issue in PDF (7.11Mb)
Download the cover in PDF (0.37Мб)
The article was published in issue no. № 2, 2016 [ pp. 124-127 ]

Perhaps, you might be interested in the following articles of similar topics: