Physical Process Modeling


    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.


Demo Movies - Heat Transfer Process:
Process of establishing of the temperature 1
Process of establishing of the temperature 2
Process of establishing of the temperature 3
Demo - heat process modeling
Heat transfer process - Isosurface 1
Heat transfer process - Isosurface 2
Heat transfer process - Streamline 1

See also:
Measurement Devices
Three-Dimensional Model Of System Furnace-Heated Body
Energy Analysis Of ERF
Process Of Melting
Cooling Of The Bodies
Regulation Of The Process Of Heating
Regulation Of The Temperature
Model Of Regulation By Assigning Hysteresis
A Model Of PID Regulation
A Model Of PID Program Regulator
A Distribution Of The Installed Capacity In ERF
Electro Resistance Furnace Part I
Electro Resistance Furnace Part II
Electro Resistance Furnace Part III
Created by Physical Process Modeling

oven modelling FEM

MONO-DIMENSIONAL AND TWO-DIMENSIONAL MODEL
/System Differential Equation/





A Summarized Model For Studying The Operation Of Electric-Resistance Furnaces

    The classical methods for analysis of electric-resistance furnaces operation are connected with usage of thermo-technical parameters on the basis of which algebraic or differential equations are made depending on the modes under investigation – whether set or transition ones with or without a load. Their solving presents the development of the processes and shows the properties of the constructed furnace for which the relevant assignment has been received. Depending on the working conditions and the requirements of the consumer, assigned are and therefore are liable to assessment:

- the characteristics of the heated body;
- the temperature, time or speed of heating;
- the construction particularities of ERF (el. resistance furnaces) influencing the process;

    The design which is done by using various methods ensures the construction and mechanic sizing of all elements but it does not assess of the development of the temperature fields in heating and cooling of the complex furnace–detail.
    A common-place method of design is the graphic-analytical which uses criteria and graphic dependencies:


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)
Nusselt, Grasshoff formula



    Some of the criteria given above are accounted by graphic dependencies which leads to a lesser preciseness, has mostly tentative character and does not contain concrete data about the running of the transition process.
    In modern conditions, this question can be solved by complex software which requires the relevant hardware and intellectual potential.
    The purpose of this project is to propose a summarized model and an appropriate approach for observation of the ERF operation. The transition process of heating taking place in the walls and the heated detail is solved by a model of parameters concentrated in the relevant elements of the system furnace–detail. For this purpose a system of differential equations (DE) and Finite Element Method (FEM) is formed. Software has been developed so that the system of DE be automated; thus need of expensive software has been eliminated. One of the primary advantages of the proposed software is its specialized character intended for the discussed questions and hence its easy usage and setting.
    The analysis is carried out under the following assumptions:

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


One-layer wall with one heater

    This case examines a one-layer wall with a heater which touches it; this way we obtain the basis substitution scheme from two resistances and a capacity (concerning the wall) and the resistance of the environment. The graph is presented when discussing thermal transfer through a solid body.
    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.

Rkl1 = 0.01;

    It should be noted that the possibility of accounting the functional dependency of the coefficient of thermal conductivity and the specific heating of the materials by the temperature is of significant importance. Further on in the project, when a three-dimensional model is discussed comparisons are made between the operation with variable and constant values where the difference is expressed more clearly.     The coefficient of thermal conductivity () 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);
R12 = (1/Lambda)*(0.2);

    Similarly to 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].

ode45, ode32 Transient process
Figure 1


One-layer wall with one heater

accounting of the resistance of the radiation between them(RklS )


    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.

Heat Transfer - Transient process
Figure 2


One-layer wall with one heater

investigation into the influence of RklS


    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.

temperature/iteration - Matlab modeling
Figure 3

furnace - Matlab modeling
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.


One-layer wall with one heater
the data for Rkl are dependent of the temperature


    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
AIR [oC] WALL [oC]
30 4050 6080100150 200300400
1010.010.6 11.512.213.4 14.517.220.0 26.734.6
209.810.4 11.212.113.3 14.517.420.2 27.035.0


    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:

Heat Transfer - Matlab modeling
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.

Multi-layer wall with one heater

    Substitution scheme in this case is a linear electric chain built on the basis of electro-thermal analogy.
    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:

a11 = -((1/((R11)*C1))+ (1/((R12+R13 )*C1)) ) ; - basic equation of the first layer of the wall
a12 = ( (1/((R12+R13)*C1)) ) ; - link of the first layer with the second layer
a21 = ( (1/((R12+R13)*C2)) ) ; - link of the second layer with the first layer
a22 = -( (1/((R12+R13)*C2)) + (1/((R14+Rkl1)*C2)) ) ; - basic equation of the second layer

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.

Heat Transfer - Matlab modeling
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:

Heat Transfer - Matlab modeling
Figure 7


    A system of DE for a three-layer wall is described as follows:

S11 = -( (1/((RklS1+R11)*C11)) + (1/((R12+R21)*C11)) ) ; - basic equation of the first layer
S12 = (1/((R12+R21)*C11)) ; - link of the first layer with the second layer
S13 = 0 ; - link of the first layer with the third layer
S21 = (1/((R12+R21)*C22)) ; - link of the second layer with the first layer
S22 = -( (1/((R22+R31)*C22)) + (1/((R12+R21)*C22)) ) ; - basic equation of the second layer
S23 = (1/((R22+R31)*C22)) ; - link of the second layer with the third layer
S31 = 0 ; - link of the first layer with the third layer
S32 = (1/((R31+R22)*C33)) ; - link of the second layer with the third layer
S33 = -( (1/((R31+R22)*C33)) + (1/((R32+Rkl1)*C33)) ) ; - basic equation of the third layer
NG11 = (1/((RklS1 + R11)*CN1)) ; - link of the heater with the first layer
NG12 = -( (1/((RklS1 + R11)*CN1)) ) ; - basic equation of the heater
S_1 = (1/(( RklS1 + R11)*C11)) ; - link of the first layer with the heater



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].

Multi-layer wall, heated body and heater

    The linear thermal chain describing the process of heating of a body is built in analogy to the way given above where the equations of the body are added to the system. When building a system of equations the following assumption is made - the heater is located evenly along the whole area of the wall; this way the link wall - heated body is not described for the moment. The substitution scheme in this case is shown here. The equations of the layers of the wall remain the same; only the equations of the body are added:

T11 = -( (1/((RklT1+Rt11)*Ct11)) + (1/((Rt12+Rt21)*Ct11)) ); - basic equation of the first layer
T12 = (1/((Rt12+Rt21)*Ct11)) ; - link of the first layer with the second layer
T13 = 0 ; - link of the first layer with the centre
T21 = (1/((Rt12+Rt21)*Ct22)) ; - link of the second layer with the first layer
T22 = -( (1/((Rt22+Rt31)*Ct22)) + (1/((Rt12+Rt21)*Ct22)) ); - basic equation of the second layer
T23 = (1/((Rt22+Rt31)*Ct22)) ; - link of the second layer with the centre
T31 = 0 ; - link of the centre with the first layer
T32 = (1/((Rt31+Rt22)*Ct33)) ; - link of the centre with the second layer
T33 = -( (1/((Rt31+Rt22)*Ct33)) ) ; - basic equation of the centre


    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).

AlfaT1= 0.0000000567* 0.7*(( y(4)+273)+(y(7)+273))*((y(4)+273)^2+(y(7)+273)^2);
RklT1 = 1/(AlfaT1);


    The equations of the heater in the two directions are as follows:

NG11 = (1/((RklS1+R11)*CN1)); - link of the heater with the first layer of the wall
NG12 = (1/((RklT1+Rt11)*CN1)); - link of the heater with the first layer of the body
NG13 = -(NG11+NG12); - basic equation of the heater
S_1 = 1/(C11*(RklS1+R11)); - link of first layer of the wall with the heater
S_2 = 1/(Ct11*(RklT1+Rt11)); - link of the first layer of the body with the heater


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).

Heat Transfer - Matlab modeling
Figure 8

    The results from the remaining ode methods are: ode45[3312981], ode23[1093556], ode23s[72], ode113[1714010]. From the realized computing process are also obtained the data of the thermal difference surface - center of the heated body. The maximal value reaches 329.84oC. The thermal difference as function of the number of iterations of the computing process is shown graphically on figure 9.

thermal shock
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).


chamber resistance furnace





TWO-DIMENSIONAL MODEL OF THE SYSTEM FURNACE-HEATED DETAIL
/System Differential Equation/



    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.

ode15s ode23s FEM
Figure 1

    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.

Finite Element Method

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:
Y1(1,k+1)= X(1,1) ; first wall first layer
Y2(1,k+1)= X(2,1) ; first wall second layer
......
Y12(1,k+1)= X(12,1) ; sixth wall second layer

Arrays containing the data for the heated body:
Yt1(1,k+1)= X(13,1) ; first wall
.......
Yt6(1,k+1)= X(18,1) ; sixth wall
YCC (1,k+1)= X(19,1) ; center of the body

Arrays containing the data for the heaters:
Tn1 (1,k+1)= X(20,1) ; heater to first wall
.......
Tn6 (1,k+1)= X(25,1) ; heater to sixth wall

numerical solving - own numbers numerical solving - own numbers



    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