Numerical Simulation of a Combined Radiation-Conduction Heat Transfer in An Electric Furnace

Numerical Simulation of a Combined Radiation-Conduction Heat Transfer in An Electric Furnace

Mamadou Alouma Diallo Department of Physics University Cheikh Anta Diop, Dakar, Senegal

Aboubacar Chedikh Beye Department of Physics University Cheikh Anta Diop, Dakar, Senegal

Cheikh Mbow

Department of Physics University Cheikh Anta Diop, Dakar, Senegal

Abstract – The aim of this article is to present a numerical procedure for the computation of the temperature distribution inside the volume of a heated material and the energy flows at its limits in an electric furnace. For the heat transfer modelling, we use a temperature dependant thermal properties model to get a nonlinear diffusion equation with a uniform radiation from the temperature heat source. The numerical solution of the combined heat transfer problem is obtained through a finite difference discretization in two dimensions with the help SCILAB scripts.

Keywords – Combined radiation and conduction heat transfer, finite difference discretization, Hottles net radiation method

1. INTRODUCTION:

   Generally, combined radiation and conduction heat transfer takes place in high temperature furnaces to achieve the desired quality of material. Radiation transfer takes place from the source to the boundary of the material while transfer by conduction takes place in the bulk of the material. The goal of solving the problems of combined radiation and conduction heat transfer is the calculation of the temperature distribution in the volume of the medium and the energy flows at its limits.

   Three important challenges are faced by the modeler who wishes to study these problems of combined transfers: the complexity and the nonlinearities intrinsic to the radiative transfers, the nonlinear conductive model induced by the properties of the material to be heated, the ignorance of the heat flow and the external surface temperature of the material to be heated.

   First, for the radiative transfer, a model that is compatible with resolution techniques and other equations governing the process has to be choosen. The model should also be reliable and able to accurately predict the radiative flux and divergence of radiative flux distributions in the medium. In addition, the model should be efficient from a computational point of view.

   The radiative transfer equation is a conservation equation of radiative energy. It is a complex integro-differential equation. There is no analytical solution available at the moment in its general form. In order to solve it, physical and mathematical approximations must be introduced.

   We can consider possible approximations under three different types of categories:

   1. Simplification of the spectral nature of the properties by means of spectral average radiative properties;
   2. Use of a radiative heat transfer model based on similar integrated flows or quantities.

   These simplifications will then greatly reduce the mathematical complexity required of the model and we will use this approach in this study.

   There are various numerical techniques for solving the radiative transfer equation, for example zonal, spherical harmonics, Monte Carlo, flow, discrete ordinate, finite volume and discrete transfer (DTM) methods, finite element method etc. It should be noted, however, that almost all the methods listed above have certain disadvantages. But it is sometimes possible to combine the characteristics of two or more methods to develop a more efficient technique for modeling radiation heat transfer in furnaces.

   In this study, we will use the Hottel zonal method. In its most simple version, it consists of decomposing a medium and its boundaries into a large number of isothermal surface element exchange areas and isothermal volume elements in order to calculate the net radiative flux exchanged between all these exchange areas [1] (Goheneche et al. Sacadura, 2002).

   For the conduction problem, many practical engineering situations require solving nonlinear transient heat transfer problems. Due to the limitations of analytical solutions for nonlinear heat transfer problems, a number of numerical methods have been developed to solve such problems. The nonlinear resulting heat equation, has applications in various branches of science and engineering, including thermal processing of materials, liquid movement in porous media.

   These non linearities in the governing equations and the boundary conditions describing the temperature distribution are due to the fact that most metallic materials have thermal properties (thermal conductivity, specific heat, and density) that are usually temperature-dependent. However, because of the difficulties associated with the solution of these nonlinear heat transfer problems, simplifying assumptions are usually made to linearize such problems. For example, in the case of materials that have thermal conductivity which varies slightly with

   temperature, constant thermal conductivity is generally assumed. However, if temperature change is substantial or the thermal conductivity varies greatly with temperature, the assumption of constant thermal conductivity may lead to significant error in the solution. Therefore, when modeling and simulating temperature distribution for such problems, nonlinearities caused by temperature dependent thermal proprieties have to be accounted for by the numerical computation.

   Concerning the solution of the combined set of equations, precise and detailed temperature distributions must be obtained so that derivatives can be accurately evaluated.

   To calculate the distribution of the temperature over the volume of the material, it is necessary to know the distribution of the intensity of the radiation if it cannot be considered uniform. In addition, according to the radiation transfer equation and its boundary conditions, it follows that in order to calculate the radiation intensity distribution; the temperature distribution must be known for all the points on the beam inside the material and at its limits. This is how the combination of radiative and conductive heat transfer is mathematically manifested.

   In combined radiation transfer and conduction calculations, the radiation transfer equation must be iteratively resolved. With the advancement of computers, it is theoretically possible to solve any partial differential equation using an appropriate numerical discretization scheme. The challenge here is to design a scheme that preserves the physics of the problem.

   Among the few earlier works on combined conduction- radiation problems, Razzaque et al. [2] analyzed the two- dimensional coupled conduction-radiation problems with finite element method limited to only non-scattering medium with isothermal black walls. The product integration method found its application in the work of Tan [3] for solving combined conduction-radiation problem in square enclosure with isothermal walls. Kim et al. [4] employed the discrete ordinates method (DOM) for coupled radiative and conductive heat transfer in rectangular enclosures. Rousse [5] and Rousse and al. [6] used the Control-Volume Finite Element Method (CVFEM) for the solution of combined mode of heat transfer in two dimensional cavities. Collapse dimension method has been used to analyse combined conduction- radiation problem by Talukdar et al. [7]. Mahapatra et al.

   [8] investigated a new hybrid method where the concepts of modified differential approximation were employed by blending discrete ordinate method and spherical harmonics method on combined heat transfer in two-dimensional planar geometry. Recently,Amiri et al. [9] analyzed the problem of combined conduction and radiation heat transfer in 2D irregular geometries by using DOM and blocked-off method with both temperature and heat flux boundary conditions.
2. MODELLING OF THE CONDUCTION HEAT TRANSFER:

   Electric furnaces consist of a heating chamber with electricity as source of heat to achieve high temperatures

   for molten metals, alloys, ceramics, refractories and so on. For metalworking, electricity generally does not have an electrochemical effect on the metal but simply heats it.

   Although modern electric furnaces are usually arc furnaces or induction furnaces, in this application we will treat a furnace in which heat is produced by resistance elements lining the inside of the furnace.

   To model the heat transfer problem, first define a number of parameters.

   0()- temperature of the radiating source (within the resistors) at time t

   [0()]- Resistivity in (ohm.m) of the nichrome 80-20 corresponding to the temperature (t)

   ()- the electrical current flowing through the electrical resistances at time t.

   We have the following formula:

   [0()]= (1 + [0() ]) (1)

   With:

   – reference temperature in ° C.

   – Resistivity in (ohm.m) of the nichrome 80-20 at the reference temperature

   temperature coefficient of nichrome 80-20

   We will choose in this study =20°C; = 1.08 106. ; = 58 106

   We will express 0() as a function of the current ()

   using a linear regression:

   We have the following data from a vendor catalog:

   Table 1: tabulation of the temperature as function of electric current in a Nichrome 80-20 resistance of 7mm diameter.

   Température (°C) (0)	200	400	600	800	1000	1200
   Courant i (A)	54,7	116	184	270	365	473

   Table 2: Results of the statistical estimation of I as function of 0

   350.3972

   Dependent Variable: I(A)
   Method: Least Squares
   Date: 12/09/18 Time: 00:50
   Sample: 1 6
   Included observations: 6
   Variable	Coefficient	Std. Error	t-Statistic	Prob.
   C	-48.66667	17.38394	-2.799518	0.0488
   0(°)	0.417786	0.022319	18.71890	0.0000
   R-squared	0.988713	Mean dependent var		243.7833
   Adjusted R-squared	0.985892	S.D. dependent var		157.2108
   S.E. of regression	18.67335	Akaike info criterion		8.953273
   Sum squared resid	1394.776	Schwarz criterion		8.883860
   Log likelihood	-24.85982	F-statistic
   Durbin-Watson stat	1.055198	Prob(F-statistic)		0.000048

   = 0.41779 0 48.667 (69) (2)

   And the correlation coefficient R = 0.99.

   The power per furnace resistance unit length () is expressed as follows:

   () = (1+[0()]) 2() (3)

   The instantaneous power generated by the furnace is then written: () = ()

   With = for an oven like the one shown in Figure 4 below.

   With

   – length of a groove

   – Number of grooves inside the oven

   The electrical energy consumed during the time T of the heating process is:

   k being the thermal conductivity coefficient of Steel steel, 99.2% Fe, 0.2% C.

   The conduction or diffusion of heat in the metal is modeled by the following equation:

   ()() . (()) = 0 (6)

   u being the temperature in the metal to be heated,

   k is the thermal conductivity coefficient of steel,

   the density of steel

   the specific heat of steel

   So the governing equations of the model are therefore the following set of equations :

   ()() . (()) = 0, (, ) × ]0, [

   (1+[0()]) 2

   (, )

   () = , (, ) × [0, ] (7)

   = 0

   () (4)

   1

   { (, 0) = 0

   The heating process occurs inside an Electric furnace which consist of a heating chamber with electricity as the heat source for achieving the required source temperatures. Because of the lack of accurate radiative properties we limit consideration to gray diffuse surfaces with uniform radiosity

   Figure 1: Schematization of the furnace with its load

   Heat is transmitted to the object by radiation. The flux transmitted by radiation becomes predominant at temperatures above 400°C. This heat transfer is performed according to the following relation:

   (,)

   In order to obtain a continuous function of the temperature u, a correlation is made between the temperature and the Thermal conductivity. We obtain the following table:

   Table 3: Thermal conductivity of steel as function of temperature

   Température (°C)	0	100	300	500	800
   Thermal conductivity (W/m.K)	45,357	45,357	43,031	37,216	30,238

   The following figure shows the variations of the thermal conductivity k as a function of the temperature u.

   50

   45

   40

   ()

   ()

    

   35

   30

   25

   400 600 800 1000

   Figure 2 : Evolution of the thermal conductivity of steel as a function of temperature

   ()

   = 1 (5)

   This relationship reflects the fact that the radiated heat is transmitted to the metal by conduction at the boundary.

   being the outer unitary normal of the surface boundary;

   The below table gives the parameters of the regression made with the Eviews software.

   Table 4: Results of the statistical estimation of k as function of temperature

   55.06929

   Dependent Variable:
   Method: Least Squares
   Date: 10/29/18 Time: 21:03
   Sample: 1 5
   Included observations: 5
   Variable	Coefficient	Std. Error	t-Statistic	Prob.
   C	52.95365	1.926324	27.48948	0.0001
   	-0.021115	0.002845	-7.420869	0.0051
   R-squared	0.948338	Mean dependent var		40.00720
   Adjusted R-squared	0.931117	S.D. dependent var		6.958590
   S.E. of regression	1.826325	Akaike info criterion		4.331662
   Sum squared resid	10.00638	Schwarz criterion		4.175438
   Log likelihood</>	-8.829156	F-statistic
   Durbin-Watson stat	1.384629	Prob(F-statistic)		0.005063

   () = 52.9536539 0.02111466019 (8)

   As above, a correlation is made between the temperature and the specific heat. We obtain the following table:

   Table 5: specific heat of steel function of temperature

   Temp (°K)	300	400	500	600	800	1000
   Specific heat (J/Kg.K)	450	491	524	555	692	1034

   The following figure shows the variations of the specific heat as a function of the temperature u.

   1100

   1000

   900

   Table 6: Results of the statistical estimation of as function of temperature

   9.098083

   83.90058

   Dependent Variable:
   Method: Least Squares
   Date: 10/29/18 Time: 21:18
   Sample: 1 5
   Included observations: 5
   Variable	Coefficient	Std. Error t	-Statistic	Prob.
   C	297.2973	28.18531	10.54795	0.0018
   u	0.471351	0.051459	9.159726	0.0028
   R-squared	0.965478	Mean dependent var		542.4000
   Adjusted R-squared	0.953970	S.D. dependent var		92.27296
   S.E. of regression	19.79671	Akaike info criterion
   Sum squared resid	1175.730	Schwarz criterion		8.941859
   Log likelihood	-20.74521	F-statistic
   Durbin-Watson stat	1.976317	Prob(F-statistic)		0.002751

   = 297.2972973 + 0.4713513514 (9)

   Applying the same procedure as above, we obtain the above table (table 7) and figure (figure 4).

   Table 7: Density of steel as a function of temperature

   Temp (°K)	300	400	500	600	800
   Density (Kg/m3)	7860	7830	7800	7760	7690

   7900

   7850

   800

   ()

   ()

    

   700

   600

   7800

    

   7750

   500

   400

   300 400 500 600 700 800 900 1000

   Figure 3: Evolution of the specific heat of steel as a function of temperature

   The below table gives the parameters of the regression made with the Eviews software.

   7700

   7650

   300 400 500 600 700 800

   Figure 4: Evolution of the density of steel as function of temperature

   The statistics of the regression are summurized below :

   Table 8: Results of the statistical estimation of as L

   L

   L

    

   Surface 2

   Surface 1

   Surface 2

   Surface 1

    

   Dependent Variable:
   Method: Least Squares
   Date: 10/29/18 Time: 21:12
   Sample: 1 5
   Included observations: 5
   Variable	Coefficient	Std. Error	t- Statistic	Prob.
   C	7966.486	5.405405	1473.800	0.0000
   U	-0.343243	0.009869	– 34.78038	0.0001
   R-squared	0.997526	Mean dependent var		7788.000
   Adjusted R- squared	0.996702	S.D. dependent var		66.10598
   S.E. of regression	3.796632	Akaike info criterion		5.795280
   Sum squared resid	43.24324	Schwarz criterion		5.639055
   Log likelihood	-12.48820	F-statistic		1209.675
   Durbin-Watson stat		1.652027	Prob(F- statistic)	0.000052

   Dependent Variable:
   Method: Least Squares
   Date: 10/29/18 Time: 21:12
   Sample: 1 5
   Included observations: 5
   Variable	Coefficient	Std. Error	t- Statistic	Prob.
   C	7966.486	5.405405	1473.800	0.0000
   U	-0.343243	0.009869	– 34.78038	0.0001
   R-squared	0.997526	Mean dependent var		7788.000
   Adjusted R- squared	0.996702	S.D. dependent var		66.10598
   S.E. of regression	3.796632	Akaike info criterion		5.795280
   Sum squared resid	43.24324	Schwarz criterion		5.639055
   Log likelihood	-12.48820	F-statistic		1209.675
   Durbin-Watson stat		1.652027	Prob(F- statistic)	0.000052

    

   function of temperature.

   L

   = 7966.486486 0.3432432432 (10)

3. MODELLING OF THE RADIATION HEAT TRANSFER :

   For the calculation of the flux q1 in equation (5), the enclosure boundary is subdivided into areas so that over each such area the following restrictions are met:

   1. All surfaces are opaque.
   2. The temperature is uniform.
   3. The surface properties are uniform.
   4. The , , and are independent of wavelength and direction so that at any surface temperature, the hemispherical total absorptivity and emissivity are equal and depend only the temperature : = = 1 where is the reflectivity.
   5. All energy is emitted and reflected diffusely.
   6. The incident and hence reflected energy flux is uniform over each individual area.

   Figure 5: Schematization of the furnace and surfaces 1 and 2

   Although no real surface is truly gray, it often happens that is relatively constant over that part of the spectrum where the blackbody emissive power is substantial, making the simplifying assumption of a gray surface warranted.

   The below tables (9) and (10) give the values of total emissivity as function of temperature. But thanks to hypothesis 4, we will consider these values as spectral emissivity as well.Furthermore, for simplification purposes and since the furnace will operate at 1000°C in this present case, we will choose 1 = 0.85 and 2 = 0.87

   Table 9: Emissivity of Oxidized Stainless Steel 303 as function of temperature

   Temperature (°C)	Emissivity
   316	0.74
   1093	0.87

   Table 10: Emissivity of Nichrome 80Ni-20Cr, oxidized as function of temperature

   Temperature (°C)	Emissivity
   100	0.87
   600	0.87
   1300	0.89

   A complex radiative exchange occurs inside the enclosure as radiation leaves a surface, travels to other surfaces, is partially reflected, and is then rereflected many times within the enclosure with partial absorption at each contact with a surface. It is complicated to follow the radiation as it undergoes this process; fortunately, this is not always necessary. A convenient analysis can be formulated by using the net-radiation method. In this method, radiative energy balances are constructed for each surface, and the resulting set of equations is then solved. This method was first devised by Hottel and later developed in a different manner by Poljak [10]. An alternative approach was given by Gebhart in [11] and [12].

   In its most simple version, the walls of the enclosure are divided into many uniform-property, uniform-temperature surface elements. Consider one element of this subdivision which we call the kth inside surface area Ak of the enclosure. The qi,k and qo,k are the rates of incoming and outgoing radiant energy per unit area of Ak.

   The outgoing radiation energy flux from a given location

   on surface k, qo,k is made up of the emitted and reflected flux from that surface :

    

   , = 4 + ,

   (11)

   where is the temperature of surface k. Note that all quantities are evaluated at a particular location on surface

   1. The quantity qi,k is the radiation flux incident at the given location from all other surfaces in the enclosure, including surface k itself, if it is concave. The quantity qo,k is often called the radiosity of the surface. It is derived from the portions of radiant energy leaving the surfaces inside the enclosure that reach the kth surface. The quantity

      enclosure are incident to any of the surfaces Ak, the incident energy is then equal to:

      , = 1,11 + 2,22 + + , + +

      , + + , (14) From configuration-factor reciprocity, we have:

      11 = 1; 22 = 2; ; =

      ; ; = (15)

      Then taking into account equation (28), Equation (27) can be written so the only area appearing is :

      , = ,11 + ,22 + + ,+ +

      , + + , (16)

      so that, by dividing both parts of the above equality by ,

      the incident flux is:

      qi,k, is often referred to as the irradiance. If the kth surface can view itself (is concave), a portion of its outgoing flux

      , = 1 ,

      (17)

      =

      =

       

      will contribute directly to its incident flux.

      Note that, contrary to the practice in most of heat transfer, these energy fluxes carry a directionality–the radiosity is the portion of the radiant energy flux with the component away from the surface, while the irradiance is the portion

      Substituting equation (13) into (12) and (17) into (12) to eliminate , will provide two energy balance equations for in terms of and ,. Therefore, we obtain:

      directed toward the surface.The net radiative heat flux

      =

      =

      (4

      ) (18)

      leaving surface k qk is the difference between the radiosity and the irradiance :

      1

      ,

      =

      =

      =

      (

      )

      (19)

      =

      = (,

      ,

      )

      (12)

      ,

      =1 ,

      =1

      ,

      ,

      This flux corresponds to the usual concept used in heat

      Application : N=2 ;

      transfer, as the net energy flux is taken as positive if in the

      = 1

      (4

      ) (20)

      direction parallel to the surface normal of the position on k.

      1 11 1

      ,1

      The quantity qk is the energy flux supplied to Ak by some

      = 2

      (4

      ) (21)

      means other than the radiation inside the enclosure, to make up for the net radiative gain or loss and thereby

      2 12 2

      ,2

      maintain the specified inside surface temperature. For example, if Ak is the inside surface of a wall of finite thickness, Qk could be the heat conducted through the wall from the outside to Ak.

      The final equation for energy transfer quantifies the

      irradiance as the sum of the radiant energies reaching a location on surface k from all other areas on the enclosure surface. This relation can have various forms depending on the degree of approximation used in the analysis.

      Equation (11) can be re-written as (noting k = k for a gray surface):

      , = 4 + , = 4 + (1 ), =

      1 = ,1 ,111 ,212 (22)

      2 = ,2 ,121 ,222 (23)

      In radiative heat transfer, a view factor is the proportion of the radiation which leaves surface that strikes surface .

      The view factor are as follows: 11 = 0, 12 = 1

      Introducing the values of 11 and 12 in equation (22), we obtain:

      1 = ,1 ,2 (24)

       

      4 + (1 ), (13)

      On the other hand, 22 = 1 21 112 = 221

      where k = 1 k = 1 k has been used for opaque gray surfaces (refer to assumption 4). Since all the outgoing fluxes from the surfaces A1, A2, , Ak, , AN in the

      Which gives the value of 21

      : 21

      = 1 ;

      2

      So the view factor 22 can be expressed as follows:

      1

      (44)

      2 1

      2 1

       

      =

      1 +( 1 1)1

      (28)

      1

      22 = 1

      Since 1 = =

      1 2 2

      2

      By replacing geometrical form factors with their

      2

      2 1

      2 1

       

      4 4

      expressions, equation (23) can be written as:

      1

      ( )

      = 1 1

      +( 1)

      (29)

      = 1 (1 1) = 1 (

      ) (25)

      1 2

      2 ,2

      ,1 2

      ,2

      2

      2

      ,2

      ,1

      And finally,

      Eliminating ,1 and ,2 between equations (20), (21) and

      =

      1 +( 1 1)

      (30)

      (24), and eliminating ,1 and ,2 between equations (24)

      1 2

      and (25), we get :

      The sign minus in formula (28) means that the flux is in the

      = 4 11

      4 + 12

      direction of ( being the outer unitary normal of the

      { 1 1 1 1

      2 2 2

      (26)

      surface boundary).

      = 1

      2 2 1

4. NUMERICAL RESOLUTION OF THE COMBINED HEAT TRANSFER :

   And this will lead to the following system of equations :

   12

   = 4 4 + 1

   This step consists of solving numerically the set of

   { 2

   2 2 1

   = 1

   1 1

   (27)

   equations (7).

   2 2 1

   Eliminating 2 between the two above equations leads to the following equation :

   12 ( 1 ) = 4 4 + 1

   1. 1. Discretization of the PDE:

         ,

         ,

          

         Let be the solution given by the model over the grid i,j including the boundary at time = and source temperature 0.

         FD-grids in our present case are square. The values of the

         2

         2 1

         2 1 1 1

         dependent variable are calculated at the nodes, while

         And this relation will lead to the following one :

         <>parameters are specified for the spacing between the nodes (node centered grid).

         Figure 6 : spatial grid for the finite difference method

         We obtain the following scheme :

         +1 = + [ ( ) ( ) + ( ) ( )]

         ,

         ,

         2 +1/2,

         +1,

         ,

         1/2,

         ,

         1,

         ,+1/2

         ,+1

         ,

         ,1/2

         ,

         ,1

         ([]4 [ ]4)

         ,

         +1, =

         0

         ,

         ,

          

         ( )

         ,

         , = 0;

         ([]4 [ ]4)

         ,

         (31)

         1, =

         0

         ,

         ,

          

         ( )

         ,

         , = + 1;

         ([]4 [ ]4)

         ,

         ,+1 =

         0

         ,

         ,

          

         ( )

         ,

         , = 0;

         ([]4 [ ]4)

          

         , ,1

         0 ,

         = , = + 1;

         ,

         ,

          

         ( )

         = u0

         = u0

          

         0

         0

          

         { ,

         With :

         +1/2, =

         +1/2, =

         (, +1, ) + (, , ) 2

         (, , +1) + (, , ) 2

         1/2, =

         1/2, =

         (, 1, ) + (, , ) 2

         (, , 1) + (, , )

         2

         ,

         ,

          

         The problem is to compute the approximate solution of the PDE (7) at location (, ) and time =

         corresponding to the control 0 with the system of equations (31).

         Now let us address the problem of computation of boundary temperatures.

         Lets consider the second equation of the above system

         (31) which we re-write as below :

         Just note that the value = in (33) is known when computing thanks to the first equation of (31).

         +1,

         +1,

          

         Now, (33) can be solved with the NewtonRaphson

         0

         0

          

         method to provide the unknown temperatures at the boundary. The initial guess for the unknown temperatures at the boundary at time can be chosen equal to the source temperature .

         For example at t = 0, just after power is put on in the furn,ace, we have = 0 = 0. We can then solve the

         4 4

         +1,

         , +

         , +

          

         1,

         =

         ([0 ] [,] )

         ,

         ,

          

         ( )

         , = 0 (32)

         equation (33) and find = 0 for the points of the

         ,

         ,

          

         boundary.

         ,

         ,

          

         At = (n = 1) the values 1 that are inside the domain

         Note that for this equation, the points on the boundary have index = 0 and the next closest points (that are not in the boundary) have index = 1. Furthermore, the symetry of the problem allows us to deal only with this boundary equation since the others are similar.

         In (32), let us set = , = and =

         (out of the boundary) will be calculated with the numerical scheme (first equation of the system (31)) but from the values 0 of the points of the boundary for step 1.

         ,

         ,

          

         ,

         ,

          

         Then, since all the points inside the domain are known, we calculate = 1 for the points of the boundary for step 1

         by solving equation (33).

         +1, 0

         ,

         And so on at step n corresponding to the time = , the

         Since ( ) = + , (where the coefficients and

         values which are inside the domain (outside the

         ,

         ,

         ,

         are found by identification with (8)), equation (32) will

         boundary) will be calculated with the numerical scheme

         yield

         =

         (44) which will lead to the following

         +

         (the first equation of the system (31)) but from the values

         ,

         ,

          

         1 of the points of the boundary for step n-1. Then after

         equation in x (which is the unknown boundary temperature to be computed) for every time t: 4 + ( + ) ( + ) 4 = 0 (33)

         ,

         ,

          

         Assuming that at time n, all the temperatures inside the heated body are known thanks to the first equation of the system (31), we are looking for = which is on the

         boundary.

         knowing all the points inside the domain, one calculates

         ,

         ,

          

         = for the points of the boundary by solving the equation (33).

      2. Numerical analysis of the scheme:

   Note that we use the same discrétisation in x as in y. Lets suppose that the solution u is twice differentiable with respect to the space variables and once with respect to the

   Vol. 8 Issue 06, June-2019

   time variable. It is relatively easy to verify that the consistency error associated to the scheme (31) is such that [() + 2] except for the boundary where one has [() + ] and consequently, the scheme is consistent.

   For the stability in , the necessary and sufficient condition is:

   1

   (34)

   2 4

   Where = 0 () = (0) if is a decreasing function of .

   C. Simulation:

   C1. Input data:

   0 = 20° (initial temperature of the load);

   0 = 1000°; = 7200

   1 = 0.85; 2 = 0.87

   Stephane Bolzmann coefficient which value is 5.67*10-8

   For the functions (), () and () ; refer to formulas (8), (9) and (10).

   = 3; (length of the furnace squared section side)

   = 1; (length of the load squared section side) n number of time steps

   = /( + 1) (space step);

   = 30 (time step);

   C2. Results:

   921.	921.	885.	861.	846.	842.	846.	861.	885.	921.	921.
   921.	817.	741.	692.	663.	654.	663.	692.	741.	817.	921.
   885.	741.	640.	573.	536.	524.	536.	573.	640.	741.	885.
   861.	692.	573.	496.	453.	440.	453.	496.	573.	692.	861.
   846.	663.	536.	453.	407.	393.	407.	453.	536.	663.	846.
   842.	654.	524.	440.	393.	377.	393.	440.	524.	654.	842.
   846.	663.	536.	453.	407.	393.	407.	453.	536.	663.	846.
   861.	692.	573.	496.	453.	440.	453.	496.	573.	692.	861.
   885.	741.	640.	573.	536.	524.	536.	573.	640.	741.	885.
   921.	817.	741.	692.	663.	654.	663.	692.	741.	817.	921.
   921.	921.	885.	861.	846.	842.	846.	861.	885.	p>921.	921.

   Figure 7: distribution of temperatures in the heated metal

   Figure 8: temperatures level curves

   Level curves are circles centered at the center of the metal. squarred section. The temperatures decrease as one gets far from the boundary. In the boudary, the temperatures are getting bigger as one gets closer to the corners. So the physics of the problem is respected.

5. CONCLUSION :

The small time step approximation combined with the finite difference discretization allowed us to compute all the required temperatures in the heated material. And furthermore, the resulting temperatures respect the physics of the problem.

REFERENCES :

1. Goheneche, J.-M., & Sacadura, J.-F. (2002). The Zone Method: A New Explicit Relation to Calculate the Total Exchange Areas in Anisotropically Scattering Medium Bounded by Anisotropically Reflecting Walls. ASME Journal of Heat Transfer, Vol. 124, pp. 696-703.1960
2. Razzaque, M. M., Howell, J. R., and Klein, D. E.: Coupled radiative and conductive heat transfer in a twodimensional rectangular enclosure with gray participating media using finite elements, JHT, 106(3), 613619, 1984.
3. Tan, Z.: Combined radiative and conductive heat transfer in two-dimensional emitting, absorbing, and anisotropically scattering square media, Int. Commun. Heat Mass Transfer, 16, 391401, 1989a.
4. Kim, T. Y. and Baek, S. W.: Analysis of combined conductive and radiative heat transfer in a two-dimensional rectangular enclosure using the discrete ordinates method, IJHMT, 34(9), 22652273, 1991.
5. Rousse, D. R. (2000). Numerical predictions of two- dimensional conduction, convection, and radiation heat transfer I. Formulation, International Journal of Thermal Sciences, Vol. 39, pp. 315-331.

   Talukdar, P. & Mishra, S. C. (2001). Transient conduction and radiation heat transfer with heat generation in a

   participating medium using the collapsed dimension method. Numerical Heat TransferPart A, Vol. 39, No. 1,

   pp. 79-100

6. Rousse, D. R., Gautier, G. & Sacadura, J. F. (2000). Numerical predictions of two-dimensional conduction, convection and radiation heat transfer II. Validation, International Journal of Thermal Sciences, Vol. 39, pp. 332-353.
7. Talukdar, P. & Mishra, S. C. (2001). Transient conduction and radiation heat transfer with heat generation in a participating

   Vol. 8 Issue 06, June-2019

   medium using the collapsed dimension method. Numerical Heat TransferPart A, Vol. 39, No. 1, pp. 79-100.

8. Mahapatra, S. K., Nanda, P. & Sarkar, A. (2005). Analysis of coupled conduction and radiation heat transfer in presence of participating medium using a hybrid method. Heat Mass Transfer, Vol. 41, pp. 890-898.
9. Amiri, H., Mansouri, S. H. & Safavinejad, A. (2010). Combined conductive and radiative heat transfer in an anisotropic scattering participating medium with irregular geometries. International Journal of Thermal Sciences, Vol. 49, pp. 492-503.
10. Poljak, G.: Analysis of heat interchange by radiation between diffuse surfaces, Tech. Phys. USSR 1(56), 555590, 1935.
11. Gebhart, B.: Surface temperature calculations in radiant surroundings of arbitrary complexityFor gray, Diffuse radiation, IJHMT, 3(4), 341346, 1961.
12. Gebhart, B.: Heat Transfer, 2nd edn., McGraw-Hill, New York, pp. 150163, 1971.