|
|
|
Making of a mathematical model of the operation of electro-resistance furnaces is based on the
laws of thermal transfer and on the processes taking place in the mode of heating as discussed
in this project. In all models described below are used substitution schemes made on the
principle of electro-thermal analogy. The analysis is made "step by step" starting from
the most trivial example and reaching up to the dimensional modeling of the device.
The making of a model and its realization through computer simulation enables a detailed exploration of the processes. By examining the process through mono-dimensional (linear) and two-dimensional substitution scheme we can observe all the particularities in the process of making the mathematical model of the furnace.
|
/System Differential Equation/ |
|
Biot's criterion (Bi)
Fourier's criterion (Fo) Nusselt's criterion (Nu) Grasshoff's criterion (Gr) Pecle's criterion (Pe) Reynolds's criterion (Re) Prandtl's criterion (Pr) |
|
- the detail is a homogenous body the differences in which are accounted by their thickness;
- the thermal transfer through radiation is assumed as the basic one when computing the process in the chamber of the furnace;
- the thermal transfer to the environment is made through natural convection;
Both the heated body and the furnace are seen as dimensional data arrays with their relevant characteristics:
i-coefficient of thermal transfer,
,S-thickness and area of the layers of the walls,
qi–power of the individual heaters, n1-number of layers of the wall of the furnace, n2-number
of layers of the detail. Through computer simulation of the process of thermal exchange the relevant temperature picture is
obtained. The software enables examination of two-dimensional and three-dimensional model.
The obtained results enable a more precise evaluation and correction during the stage of design and are used
for correcting some of the parameters – geometric and power characteristics.
| The resistance of convection and radiation toward the environment is accepted as a constant value. In the next examples this part of the thermal chain is discussed with the real values after which analyzed is the difference in the results obtained in a constant mode and in a mode dependent on the momentary temperature. )
is obtained as function of the temperature where to each iteration a new value is given depending on the momentary temperature.
In catalog data is usually presented as a value in the
form const1+const2*Tsr, where sr is the
average temperature. In this case the average temperature is substituted by the momentary temperature y(1). Lambda = 0.520+0.349*0.001*y(1); The values of resistances are obtained depending on the momentary value of
and the geometric parameters of the wall. The general syntax is of the following form: Rnm = (1/Lambda)*(dnm/Snm) ; Rnm - consecutive thermal resistance in the system Lambda - coefficient of thermal conductivity ( )dnm - thickness of the n-th layer in the m-th side Snm - area of the n-th layer in the m-the side (in the case of a one-layer wall as it is discussed here, the indexes (n,m) for d and S make no sense but this syntax is used further on, too, when three-dimensional models of multi-layer furnaces are discussed) ; - according to the required syntax of Matlab, a line ending in a semicolon does not appear on the screen which is convenient in computing processes with multiple iterations. The quotient dnm/Snm now and in the following equations is presented by its result. For the discussed version the resistances of the wall are obtained in the form: R11 = (1/Lambda)*(0.2);
Similarly to R12 = (1/Lambda)*(0.2); the specific thermal capacity is accounted
by catalog data in the form const1 + const2 *Tsr. The latter will be discussed as function of the momentary
temperature y(1). The thermal capacity is expressed by the product of the mass of the material and the specific thermal
capacity. The mass of the material of the relevant node is obtained from the product of the wall volume and specific thickness.
Due to its triviality these computations will not be presented.C11 = 200 * (837+0.41*y(1));
The system discussed here is described by one differential equation which is an only member of the matrix yp:S1 = -((1/(( R11)*C11))+ (1/((R12+Rkl1)*C11)) ) ;
In this case the matrix has only one value and disturbing influence (q1/C11) of a heater of certain power q1.yp = [ S1*y(1) + (q1/C11) ]; In the case of a heater in close contact with the wall the description of the resistance of radiation between the two elements is avoided. This is done further on when the heater is described as an independent element in the thermal chain with its own capacity. The reason for this assumption is to render opportunity for observation the description of the system in detail by describing all the steps of the modeling process. The computing process of the model thus described is realized by the MatLab program. In Figure 1 are shown the graphs obtained from the different methods for numerical computing of the differential equations. The results present the process of heating of one-layer wall with one heater the temperature being function of the time. As it can be seen, for the conducted simulation ode the methods show similar results. Differences are observed in exploring the modes of regulation when assigning different conditions. This comparison is made in order to determine an optimal computing method depending on the description of the system. The numbers in the square brackets show the number of the iterations realized in the process (designations used). The chosen power of the heater is 2[kW].  Figure 1 The case examines a one-layer wall with a heater located in the vicinity to it with accounting of the resistance of radiation between them (RklS1). The thermal substitution scheme and the way of its obtaining is discussed in the analysis thus made of the laws of thermal transfer. The initialization of RklS1 is consistent with the deduction of the equations thus made in radiation thermal exchange in order to realize them through the applied software. AlfaS1= 0.0000000567* 0.7*(( y(2)+273)+(y(1)+273))*((y(2)+273)^2+(y(1)+273)^2); 0.0000000567 - Stephan - Boltzman's constant. 0.7 - specific coefficient of blackness. y(2), y(1) - respectively momentary temperature of the heater and monetary temperature of the wall. RklS1 is expressed by: Sst - the area of the wall (in this case 1[m2]) As in the example above, the resistance of convection and radiation toward the environment is a constant (Rkl1=0.01). The values of resistances of the wall R11 and R12, the coefficient and the thermal capacity C11
remain unchanged, i.e. the function of the momentary temperature. The basic equation of the wall obtains
the following form: S1 = -((1/((RklS1+R11)*C11))+ (1/((R12+Rkl1)*C11)) ) ; The equations of the heater are as follows: NG11=(1/((RklS1+R11)*CN1)); - The equation of the link between the heater and the wall. As a component of the matrix yp is a21 i.e. row second - basic equation of the heater, column first - basic equation of the wall: NG12=-((1/((RklS1+R11)*CN1))); - Basic equation of the heater (a22). This, as well as all the next basic equations of the elements of the thermal chain are placed on the main diagonal of the matrix: S_1= 1/ (C11*(RklS1+R11)); - Equation of the link between the wall and the heater. As a component of the matrix yp is a12, i.e. row first-basic equation of the wall, column second-basic equation of the heater. The thermal capacity CN1 is insignificant due to which it is accepted as a constant value.
yp = [
S1*y(1) + S_1*y(2) + 0
NG11*y(1) + NG12*y(2) + (q1/CN1)
] ;
In this case the disturbing influence (q1/CN1) presents the ratio of the power to the capacity of the heater.
The small value of CN1 compared with the value of the thermal capacity of the walls of the furnace determines the large number
of iterations with ode15s method.  Figure 2 The presentation of RklS1 as a constant value gives some idea about the influence of this part of thermal chain in the entire process of heating of the system wall - heater as discussed here. The resistance as function of the temperature change its value the data being presented graphically on figure 3. The dependency is the function of RklS1 of the number of computing iterations. The constant value is initialized with the obtained average value of the same data RklS1 = 0.09555485.  Figure 3  Figure 4 The remaining equations of the heater and the wall remain without any change. The obtained temperature picture is shown on of figure 4.
The resistance of the convection and radiation of the wall toward the environment is expressed by the sum of the respective coefficients. The convective thermal exchange is a complex process of exchange of heat and mass for the description of which criteria are used obtained on the basis of experimental results. When making a model with a concrete assignment - a resistance furnace or, as it is in this case - part of it, it is convenient to apply step-by-step approximation. The value of Rkl concerning the convection is updated for each iteration depending on y(1) according to the data given in table 1: Table 1
The thermal exchange through radiation is expressed according to the syntax required by Matlab where: AlfaL = 0.0000000567* 0.7 *(( y(1)+273)+(273))*((y(1)+273)^2+(273)^2 ); For Rkl it should be written: Rkl = (1/(AlfaK+AlfaL)*S ) ;
AlfaK - the coefficient of convection AlfaL - the coefficient of radiation S - the area of the wall. In this case there is possibility for studying the influence of the two coefficients in the transition process of establishing of the temperature. The graphs:  Figure 5 are obtained with accounting only of AlfaK in the equation for Rkl. In this case the temperature of the heater is established at 1150.2oC, and the temperature of the wall at 591.7oC. Performing the same simulation with accounting only of AlfaL leads correspondently to temperatures 1151.3oC and 593.27oC. Due to the small difference in the obtained values, graphs are not to be given. The computation of Rkl, on accounting the sum of the coefficients of radiation and convection, leads to establishing the temperature of the heater at 1136.7oC and of the wall at 572.88oC. In the so made description the data are incorrect since y(1) is the temperature in the center of the wall and not on its surface. The process thus described can be used for preliminary (tentative) computations. As it has been noted when accounting of the dependency Rkl=f(T) it is necessary to find the momentary temperature of the surface of the medium located at the side of the environment. For obtaining correct data the following approach is applied: the wall is divided into three layers the thickness of the end ones being determined by additional accounts (for example one millimeter). Thus the temperature of its two surfaces is obtained with an insignificant error. The forming of system of differential equations for a multi-layer wall with N layers is done layer by layer with the relevant basic equations and links. The accepted approach for finding the temperature of the two surfaces of the wall requires making a system which describes a three-layer wall. The examining of a two-layer wall aims at observing the processes of heat passage and clarifying the mechanism of forming the differential equations. The process described below excludes the heater as a separate element of the thermal chain adding the disturbing influence in the first row of the matrix, i.e. to the first layer of the wall. The description is analogous to this which was done for one wall at the beginning of the chapter. In this case the emphasis is placed on the number of iterations obtained by the different ode methods. A similar approach is used later too, when discussing the processes of regulation for the purposes of a substantial shortening of the computing time. The equations are made in the following way:
yp = [ a11*y(1) + a12*y(2) + (q1/C11) a21*y(1) + a22*y(2) + 0 ]; The graphs of the results from the different methods are presented in Figure 6. The methods for numerical solving of the differential equations ode45, ode23, ode113 in case of accounting of the parameters of the heater carry out the computing process with thousands of iterations. As it can be seen, when the heater is switched off for computing a system of DE, the different methods work with a similar number of iterations. The result, which is confirmed when examining a three-dimensional model with multi-layer walls, is that the number and sizes of the walls influence slightly the computing time of the process, the primary influence being the data for the different nodes in the chain.  Figure 6 In the beginning of the process of heating the temperature of the outward layer of the wall decreases which is shown on the graph:  Figure 7 A system of DE for a three-layer wall is described as follows:
yp = [
S11*y(1) + S12*y(2) + S13*y(3) + S_1*y(4) + 0
S21*y(1) + S22*y(2) + S23*y(3) + 0 + 0
S31*y(1) + S32*y(2) + S33*y(3) + 0 + 0
NG11*y(1) + 0 + 0 + NG12*y(4) + (q1/CN1)
] ; The basic equations in the main diagonal are marked in red. It should be noted that it is of great importance that the diagonal should not be disconnected, i.e. each of its elements should be initialized with the consecutive basic equation: ode15s, ode23, ode23s, ode45, ode113. The performed simulation gives the picture of heating of the system as well as of the thermal loading of the wall. The beginning of the process in the individual elements of the thermal chain: outward layer, heater 400[s], heater 40[s], heater 4[s].
The resistance of convection and radiation between the heated body and the heater RklT1 is determined by analogy to RklS1 with accounting of the momentary temperatures of the heater (in this case y(7)) and the first layer of the body y(4).
yp = [
S11*y(1) + S12*y(2) + S13*y(3) + 0 + 0 + 0 + S_1*y(7) + 0
S21*y(1) + S22*y(2) + S23*y(3) + 0 + 0 + 0 + 0 + 0
S31*y(1) + S32*y(2) + S33*y(3) + 0 + 0 + 0 + 0 + 0
0 + 0 + 0 + T11*y(4) + T12*y(5) + T13*y(6) + S_2*y(7) + 0
0 + 0 + 0 + T21*y(4) + T22*y(5) + T23*y(6) + 0 + 0
0 + 0 + 0 + T31*y(4) + T32*y(5) + T33*y(6) + 0 + 0
NG11*y(1) + 0 + 0 + NG12*y(4) + 0 + 0 + NG13*y(7) + (q1/CN1)
] ;
The zeroes in the matrix determine the lack of a link between the relevant components. The disturbing influence is
recorded in the last column of the line of the basic equation of the heater. The remaining elements of the same
column are zeroes, i.e. the influence refers to the heater only. The obtained results from the computing process
realized with the method ode15s are shown in figure 8
(walls of the furnace with heater,
layers of the body with heater). Figure 8  Figure 9 Furnace Gallery
The gallery shows furnaces made by various manufacturers. The analysis suggested by us refers to
similar constructions of an electric-resistance chamber furnace (ECF).
The realization of a two-dimensional model of the system furnace-detail is described in a way analogous to the linear version, the links between the separate walls being added for each layer. Without accounting of the links between the four walls, the process is to develop in such a way as it would in four detached linear systems. In case of heaters with equal power, the result will be the same as that of a mono-dimensional model. The matrix which does not account the links between the walls is shown here. This case presents interest only for clarifying the building of the matrix; for this reason the results will not be given. In correct description with accounting of the links between the layers of the walls, the system of differential equations and the matrix assume the form shown here. Figure 1 is built from the results obtained by the explained two-dimensional model of a furnace without a heated body. The description of the examined system with a heated body leads to a system of equations and a matrix presented here. The simulation performed with the method ode15s yields results whose graphs are given in figures: two-dimensional model of a furnace without a heated body, walls of the furnace with heater, layers of the body with heater.
The outward layers of the walls of the furnace preserve the character obtained by multi-layer linear scheme of wall with concern to the beginning of the process of heating: beginning of the process of heating , the outward layer of the wall , layers of the body. The descriptions thus made of a linear and two-dimensional substitution scheme through a system of differential equations give a general idea of the operation of the system furnace – heated body. In fact, the obtained results serve for studying the laws of thermal transfer through a one-layer and multi-layer flat wall and are a preparation for the description of a three-dimensional model of electric resistance furnace. The primary purpose of the project – investigation of the operation of electro-technological devices through forming a mathematical model and its realization as a computer simulation could be achieved by means of a dimensional presentation of the thermal unit.
 Figure 2
|
| Realization Of The Computing Process Through Finding The Own Numbers And Vectors Of Matrices
The mathematical realization of the topic discussed here is implemented in chapter Numerical Methods.
In the three-dimensional description of the system furnace – detail as examined here, the matrices are made according to the syntax necessary for solving the differential equations through the methods in MatLab program (files ode45.m, ode15s.m, ode23.m, ode23s.m, ode113.m). The numerical solving of systems DE for the problematics under observation here can be done by using the functions for finding of own vectors of matrices, matrix exponents, the method of the variables of the condition, realized in the same program. The following changes are necessary: the plus should not be written between the elements in the rows of the matrix, the elements should not be multiplied by y(*). After such an editing the file containing the matrices obtains the form shown here. It is necessary to add the array of the initial conditions for the process of heating. As a private case, if the temperature of all elements of the thermal chain is 20oC, the linear array X0 which assigns the starting values of the simulation for a two-layer furnace and a body with two layers is made in the following way:
X0=[ 20;20; 20;20; 20;20; 20;20; 20;20; 20;20; starting conditions for the walls of the furnace 20;20;20;20;20;20;20; starting conditions for the body – 6 sides + center 20;20;20;20;20;20; starting conditions for the heaters ]; The disturbing influences are given in a separate matrix: MATRIX_B = [ 0;0;0;0;0;0;0;0;0;0;0;0; to the walls of the furnace 0;0;0;0;0;0;0; to the walls of the body (q1/CN1); (q2/CN2); (q3/CN3); (q4/CN4); (q5/CN5); (q6/CN7); to the heaters ]; The computing process is realized through a cycle in which the data are computed from each iteration through several functions. The time of simulation is obtained by the product of the number of iterations of the cycle and the step dt. Detailed information about the operation of the functions used in Matlab is given here. The computing cycle is written in the following way:
eig(MATRIX_A)% checks whether the basic matrix is square
dt = 2; % step of the process
for k = 1: 2000, % for cycle according to the Matlab syntax
Time(k+1)=k*dt; % array of the running time. The values in Time() serve for drawing up of the graphs
t=Time(k+1); % temporary variable
E1=expm(MATRIX_A.*t); % computes the matrix exponent using Pade’s approximation
X = E1*X0+(E1-eye(size(MATRIX_A)))*inv(MATRIX_A)*MATRIX_B ;
% the obtained values are written in the basic array X. The functions in use are: eye–it forms a single matrix
% (along the main diagonal are located ones, and the remaining members are zeroes)
% whose size is obtained from the function size, inv – finds the inverse matrix of the assigned one
end % end of the cycle
It is convenient that the data for each of the elements of thermal chain be written in separate arrays so as the operation with them could be easier when they are examined and the graphs drawn. For the example discussed here, in the body of the cycle the initialization of the following arrays should be added:
Arrays containing the data for the walls of the furnace:
Arrays containing the data for the heated body:
Arrays containing the data for the heaters:
Evidently, the choice of step plays a crucial role in obtaining correct data at the beginning of the process. A question remains: how far goes the influence of the step upon the end results. Finding the answer to this question will allow determining of the possible limits of the assigned value and thence the necessary processor time for carrying out of the simulation. According to the results given here the choice of step hardly influences the end result. When discussing the three-dimensional model a comparison is made between the process in constant and variable l and thermal capacity. In order to obtain the data depending on the momentary temperature it is necessary that to each iteration of the cycle new values be given consistent with the obtained momentary temperature. The whole file, containing the initialization of the thermal resistances and capacity, system of differential equations and the matrices, is added into the body of for cycle. The shortening of the processor time is achieved through removing the heaters from the system of DE. For this purpose their temperature is computed by the dependency:
Tn1(1,k+1)=(q1+(X(1,1)/(RklS1+R11))+(X(13,1)/(RklT1+Rt11)))/((1/(RklS1+R11))+(1/(RklT1+Rt11))); Tn1(1,k+1) – array of obtained data X(1,1) – momentary temperature of the first layer of the wall X(13,1) – momentary temperature of the first layer of the heated body. |
| Physical Process Modeling Resources: Mathematics Physics Electronics Programming Heat transfer |
