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

Download Full-Text PDF Cite this Publication

Text Only Version

 

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)20040060080010001200
    Courant i (A)54,7116184270365473

    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
    VariableCoefficientStd. Errort-StatisticProb.
    C-48.6666717.38394-2.7995180.0488
    0(°)0.4177860.02231918.718900.0000
    R-squared0.988713Mean dependent var243.7833
    Adjusted R-squared0.985892S.D. dependent var157.2108
    S.E. of regression18.67335Akaike info criterion8.953273
    Sum squared resid1394.776Schwarz criterion8.883860
    Log likelihood-24.85982F-statistic
    Durbin-Watson stat1.055198Prob(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)0100300500800
    Thermal conductivity (W/m.K)45,35745,35743,03137,21630,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
    VariableCoefficientStd. Errort-StatisticProb.
    C52.953651.92632427.489480.0001
    -0.0211150.002845-7.4208690.0051
    R-squared0.948338Mean dependent var40.00720
    Adjusted R-squared0.931117S.D. dependent var6.958590
    S.E. of regression1.826325Akaike info criterion4.331662
    Sum squared resid10.00638Schwarz criterion4.175438
    Log likelihood</>-8.829156F-statistic
    Durbin-Watson stat1.384629Prob(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)3004005006008001000
    Specific heat

    (J/Kg.K)

    4504915245556921034

    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
    VariableCoefficientStd. Error t-StatisticProb.
    C297.297328.1853110.547950.0018
    u0.4713510.0514599.1597260.0028
    R-squared0.965478Mean dependent var542.4000
    Adjusted R-squared0.953970S.D. dependent var92.27296
    S.E. of regression19.79671Akaike info criterion
    Sum squared resid1175.730Schwarz criterion8.941859
    Log likelihood-20.74521F-statistic
    Durbin-Watson stat1.976317Prob(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)300400500600800
    Density (Kg/m3)78607830780077607690

    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
    VariableCoefficientStd. Errort- StatisticProb.
    C7966.4865.4054051473.8000.0000
    U-0.3432430.009869– 34.780380.0001
    R-squared0.997526Mean dependent var7788.000
    Adjusted R- squared0.996702S.D. dependent var66.10598
    S.E. of regression3.796632Akaike info criterion5.795280
    Sum squared resid43.24324Schwarz criterion5.639055
    Log likelihood-12.48820F-statistic1209.675
    Durbin-Watson stat1.652027Prob(F- statistic)0.000052
    Dependent Variable:
    Method: Least Squares
    Date: 10/29/18 Time: 21:12
    Sample: 1 5
    Included observations: 5
    VariableCoefficientStd. Errort- StatisticProb.
    C7966.4865.4054051473.8000.0000
    U-0.3432430.009869– 34.780380.0001
    R-squared0.997526Mean dependent var7788.000
    Adjusted R- squared0.996702S.D. dependent var66.10598
    S.E. of regression3.796632Akaike info criterion5.795280
    Sum squared resid43.24324Schwarz criterion5.639055
    Log likelihood-12.48820F-statistic1209.675
    Durbin-Watson stat1.652027Prob(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
    3160.74
    10930.87

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

    Temperature (°C)Emissivity
    1000.87
    6000.87
    13000.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. 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.

Leave a Reply

Your email address will not be published. Required fields are marked *