Inverse boundary design of a radiative smelting furnace with ablative phase change phenomena

Applied Thermal Engineering 98 (2016) 1140–1149

Contents lists available at ScienceDirect

Applied Thermal Engineering j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / a p t h e r m e n g

Research Paper

Inverse boundary design of a radiative smelting furnace with ablative phase change phenomena H. Farzan, S.M. Hosseini Sarvari *, S.H. Mansouri Mechanical Engineering Department, Shahid Bahonar University of Kerman, 76175-133, Kerman, Iran

H I G H L I G H T S

• • • •

The ablation phenomenon in a reverberatory smelting furnace is simulated numerically. The results are verified by comparing with exact analytic solution. Inverse design problem is solved to construct the desired melting rate. The conjugate gradient method is used to solve the inverse phase change problem.

A R T I C L E

I N F O

Article history: Received 10 September 2015 Accepted 6 January 2016 Available online 15 January 2016 Keywords: Inverse design Ablative phase change Smelting furnace Conjugate gradient method

A B S T R A C T

An inverse analysis is employed to control the time rate of heaters in a 2-D smelting furnace to provide the specified radiative heat flux across the design surface to establish a desired melting rate. The design surface in the smelting furnace is the melting surface of the metal concentrate bank, and the melting process is considered to occur as an ablation phenomenon. The net radiation method is used to determine the radiation exchange between the elements of the furnace surfaces and the melting surface. The conjugate gradient method is employed to minimize the objective function, which is the sum of square residuals between the estimated and the desired heat fluxes over the design surface. It is shown that the proposed inverse technique is reliable and accurate for predicting the heater power distribution. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Phase change phenomena are a crucial part of metal industries such as steel and copper industries. Therefore, achieving the desired parameters inside the smelting furnaces such as temperature or heat flux distribution, suitable solid thickness on the furnace wall, and melting rate are appealing subjects to the researchers. Appropriate melting rate across the melting surface in radiant smelting furnace has significant effect on the efficiency and the performance of furnace. The thickness of protective banks which are formed on the furnace wall and the progressive solid/liquid interface can be determined by controlling the heater powers. Measuring the great part of these parameters is technically difficult, time consuming and expensive. However, numerical techniques have proposed significant methods to handle these problems. Among the numerical methods, inverse heat transfer techniques released suitable procedures to obtain the desired parameters. In the previous studies, inverse heat transfer methods were used to obtain the heat flux load with known location of solid/liquid in-

* Corresponding author. Tel.: +98 341 2111763; fax +98 341 2111763. E-mail address: [email protected] (S.M. Hosseini Sarvari). http://dx.doi.org/10.1016/j.applthermaleng.2016.01.029 1359-4311/© 2016 Elsevier Ltd. All rights reserved.

terface [1,2]. Moreover, heat flux load inside the furnace was determined with the temperature measurement in the phase change material by many researchers [3,4]. Among the methods that were developed to estimate the unknown parameters in heat transfer problems, the conjugate gradient method (CGM) or the Levenberg– Marquardt method (LMM) were widely used [5–7]. These methods propose appropriate and straightforward procedure; in addition, implementation of such methods is computationally efficient and inexpensive. Howell et al. [8] applied the inverse heat transfer analysis to determine the approximate initial design of complicated tools. Federov et al. [9] and Duan et al. [10] applied inverse heat transfer techniques to obtain desired heat flux or temperature distribution on the design surface. Hosseini Sarvari et al. [11–15] utilized inverse methods to determine desired heat flux and temperature distributions on the design surface in two and three-dimensional enclosures with absorbing–emitting radiant media. The optimal design of a radiant enclosure with genetic algorithm was investigated by Hosseini Sarvari [16]. Payan et al. [17] used the CGM to determine the strength of heaters over the heater surface of a radiant enclosure with natural convection to create the desired condition on the design surface. The desired temperature and heat flux distributions over the design surface in a channel formed by two parallel plates with forced convection and surface radiation exchange were

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

obtained by estimation of the unknown heat flux distribution over the heater surface [18]. Payan et al. [19] used the CGM to obtain the temperature distribution near the flame by estimating the heat flux over the opposite wall in the direction of flame diffusion. Wang et al. [20] estimated the geometry of inner wall surface of a furnace by the CGM. One of the industrial applications of a radiant furnace is the melting of metals. Inverse heat transfer techniques were applied to phase change problems to obtain the unknown boundary condition. Even though, due to the nonlinearity of ablation problems and the ill-posed nature of inverse problems, few studies have been done in this field. Petrushevsky and Cohen [21] used in-depth temperature measurement to obtain the heat flux and recession rate history. They calculated the heat flux at the ablating surface, and then the recession rate was approximated by the estimated heat flux. Oliveira and Orlande [22] estimated the heat flux at the surface of ablating material by using temperature and surface position measurements. The CGM was used to estimate the heat flux distribution over an insulated ablator surface by Hakkaki-Fard and Kowsary [23]. Molavi et al. [24] solved one dimensional ablation problem to achieve the heat flux at the receding surface of ablative material with temperature-varying thermal properties. LeBreux et al. [25] used Kalman filter to predict the time-varying bank thickness on the wall surface of a high temperature furnace. An inverse analysis was employed to determine the unknown effective specific heat through the temperature data collected in a muffle furnace [26]. The previous studies are classified in two main categories. The first category is concerned with the inverse phase change problem without considering the power source which causes phase change. Another category is focused on the inverse design problems to produce desired thermal conditions (desired temperature and heat flux distributions) over the design surface. In the present study, we try to couple these two categories to solve a practical inverse design problem. The inverse heat transfer analysis is applied to solve a practical heat transfer problem to obtain an appropriate set of heaters over the heater surface of a radiant furnace in order to achieve a desired melting rate. The major heat transfer mechanism inside the furnace is assumed to be radiation and the phase change along the melting surface is considered to occur as an ablation phenomenon. Therefore, the melting process is assumed to be isothermal over the ablative surface. The radiation problem inside the enclosure is solved by utilizing the surface radiation estimation and the furnace media is supposed to be transparent. 2. Description of problem A two dimensional radiant furnace, demonstrates at Fig. 1. The furnace medium is transparent, and all the walls are gray-diffuse. This geometry can be considered as a simplified smelting furnace whose task is melting the concentrate bulk of metals. The heaters are located at the top and the left walls of the furnace. The walls are adiabatic and the slag surface is isothermal at the molten slag temperature. The goal of the inverse problem is to find a set of heaters over the heater surface to create the desired uniform melting rate along the melting surface (design surface). 3. Direct problem The direct problem consists of three parts: the radiation problem, the conduction problem, and the ablative phase change problem. 3.1. Radiant heat transfer problem in the furnace Radiation exchange between the furnace’s surface elements is simulated by the net radiation method (NRM). The equation of ra-

1141

Fig. 1. The schematic shape of smelting furnace.

diation exchange for the surface elements with specified temperature is described by the following equation: NS

∑ [δ

kj

j =1

− (1 − ε k ) Fk− j ] J j (t ) = ε kσ Tk 4 (t )

(1)

whereas the radiative transfer equation for other surface elements with specified heat flux is given by: NS

∑ [δ j =1

kj

− Fk− j (t )] J j (t ) = ϕ k (t )

(2)

Where Fk− j is the view factor, which is obtained by the Hottel’s crossed string rule [27]. By solving the linear system of equations achieved by Eqs. (1) and (2), the radiosities that leave the surfaces are calculated. Then the unknown boundary condition (temperature or heat flux) is determined by the following equation:

1− ε k ϕ k (t ) + J k (t ) = σ Tk4 (t ) εk

(3)

In Eqs. (1) and (2), the δ kj is the Kronecker delta defined by

⎧1 k = j δ kj = ⎨ ⎩0 k ≠ j

(4)

3.2. The conduction problem Regarding the geometry of the smelting furnace and the entrance location of metal concentrate, the produced liquid layer may remove from the solid bank as soon as the phase change occurs. Therefore, the convection mechanism may be ignored, as compared to radiation and conduction mechanisms. Hence the ablation is a suitable approximation for phase change phenomenon. The solid bank in the smelting furnace is restricted by the walls from the right wall and slag bath from the bottom. The walls are assumed to be adiabatic, and the slag bath is considered to be isothermal at the melting temperature. The dominant mode of heat transfer through the solid zone is conduction and the radiation load is applied on the ablative interface. As the ablated material is removed from the bank surface, the ablative surface is considered at constant melting temperature. The solid bank properties are assumed to be constant and homogenous during ablation process. The governing heat transfer equation to determine the temperature distribution and heat fluxes in the solid bank is as follows:

1142

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

∂T = α∇ 2T ∂t

(5)

where T is temperature, and α is the thermal diffusivity. Eq. (5) with its subjected boundary conditions is solved by the finite element method (FEM). In the FEM, the domain of interest is discretized by a triangular mesh. Then the unknown temperature distribution T ( x, y, t ) is approximated over the entire domain by: 3

T ( x, y, t ) = ∑ T jγ j ( x, y )

(6)

j =1

where T j ’s are the temperatures at nodes with coordinates ( x, y ) . The γ j ’s are the shape functions which interpolate the solution between the discrete values obtained at the nodes with the following property:

γ i ( x j, y j ) = δ ij

(7) Fig. 2. Receding of ablating surface in two consecutive time steps.

The FEM discretization of Eq. (5) for each element is as follows:

⎡ Kˆ ⎤ ⎣ ⎦

(e )

{T }(e) = { fˆ }

(e )

(8) termined. The governing heat transfer equation along the ablative surface is given by:

where

⎡ Kˆ ⎤ ⎣ ⎦

{ fˆ }

(e )

(e )

(e )

ˆ ⎤ + ωΔt ⎡ Kˆ ⎤ = ⎡⎣ M ⎦ ⎣ ⎦

(e )

(

(9a)

){ }

ˆ ⎤(e) − (1 − ω ) Δt ⎡ Kˆ ⎤(e) Tˆ = ⎡⎣ M ⎣ ⎦ ⎦

(e )

old

(

{}

+ Δt ω fˆ

(e )

{}

+ (1 − ω ) fˆ

(e ) old

)

(9b) (e ) Here, ⎡ Kˆ ⎤ is the stiffness matrix which contains the contribu⎣ ⎦ (e ) is the tion of boundary conditions and the shape functions, Tˆ

vector of unknown temperatures, and

{} fˆ

(e )

{}

is the thermal load

vector. Δt is the time step and the subscript old indicates the temperatures evaluation at previous time step. The time parameter (0 ≤ ω ≤ 1) determines the solving procedure. By selecting ω = 1 the solving method is implicit and by choosing ω = 0 the method is completely explicit. The element of defined matrix and vectors are:

K ij(e) = ∫

Λ (e )

Mij(e) = ∫

⎛ ∂γ i ∂γ j ∂γ i ∂γ j ⎞ ⎜⎝ ∂x ∂x + ∂y ∂y ⎟⎠ dxdy

(10a)

⎛ γ iγ j ⎞ ⎜⎝ ⎟ dxdy α ⎠

(10b)

Λ (e )

fi(e) = ∫

Γ (e )

ψ (e)γ i ds

(10c)

3.3. The ablation problem The solid bank begins to smelt as the result of radiation heat transfer exchange from the heaters located on furnace walls. Through the concentrate bank, the phase change occurred on those surfaces which are exposed to radiation. The molten materials then removes away from the exposed surface and drain to the melt bath at the bottom furnace. At each time step, the exposed ablative surface is considered at isothermal melting temperature, Tmelt . A fully implicit discretization of Eq. (8) is used to determine the temperature distribution and temperature gradient in the solid zone with moving grid approach. Obtaining the temperature gradient at melting surface, the amount of conductive heat fluxes toward the solid zone is de-

h

ds = ϕ −ψ dt

(11)

where h is the latent heat of melting, and

ψ = −k

∂T ∂n

(12)

is the conductive heat flux. Computing the conductive and radiative heat fluxes over exposed ablative surface, the recession rate of m-th surface element over the ablative interface is calculated by the finite difference formula as:

Δsm =

Δt (ϕ m − ψ m ), m = 1, … , M h

(13)

where ϕ m and ψ m are the radiative and conductive heat fluxes over the m-th surface element over the ablative surface. A least square approach based on the RCURVE routine from the IMSL library of FORTRAN is used to approximate the ablative surface at each time step. The schematic of melting surface recession toward the solid zone is shown in Fig. 2 in two consecutive time steps due to radiative heat load. 4. Direct problem solution procedure The solution domain of the solid zone has an irregular boundary, therefore solving the conduction heat transfer using FEM requires a triangular mesh in the solid zone. Moreover, the solid zone must be remeshed at each time step to track the interface. A fast and efficient method based on the algorithm developed by Lo [28] is used for triangulation of the solid zone. Due to remeshing of the solid (e ) zone, the initial temperature distribution, {Told }, from the previous time step must be found for the new grid points. Fig. 3 shows the triangular mesh for two consecutive time steps. The initial temperature for grid P at each time step is updated by the following equation: 3

1 TPν,+old = ∑ T jν,oldγ j ( xP , yP ) j =1

(14)

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

1143

should be controlled, are the time rate of heater powers along the heater surface. The desired parameter, nodal melting rate, is related to the radiative heat flux and conductive heat flux over the design surface by Eq. (13). As the ablative interface is assumed as an isothermal surface at melting temperature, and all the walls are adiabatic, the conductive heat flux through the moving surface may be obtained by solving the conduction problem, defined by Eq. (5). Then knowing the desired ablative rate, Δs Δt , the radiative heat flux may be obtained by Eq. (13). Hence, the inverse problem is simply reduced to solve the radiation inverse problem to estimate the radiative heat flux over the design surface (ablative interface), by an appropriate set of heaters over the heater surface. The desired and estimated radiative heat fluxes over the design surface, respectively, can be expressed by two vectors as below:

Fig. 3. Mesh regeneration of solid zone in two consecutive time steps. Red solid lines and black dashed lines show the new and old meshes, respectively.

where subscript j shows the node numbers of the triangular element for the previous time step that surrounds the node P . The computational algorithm for the direct solution problem can be summarized as follows: Step 1. Specify the final time, t f , and the time increment, Δt , Set t = 0. Step 2. Triangulate the solid zone, divide the furnace walls to the N S surfaces and apply the boundary conditions for the solid zone and the furnace walls. Note that the number of surfaces over the melting interface must be equal from the both side. Step 3. Solve the radiation heat transfer equations for the furnace walls and find the radiative heat flux distribution along the melting interface, using Eqs. (1)–(4). Step 4. Solve the heat transfer problem in the solid zone and find the conductive heat flux distribution along the melting interface, using Eqs. (5)–(10). Step 5. If the time is less than t f continue, otherwise stop. Step 6. Given the radiative and conductive heat flux distributions along the melting interface, update the moving boundary by using Eq. (13). Step 7. Set t new = t old + Δt . Step 8. Given the melting interface from step 6, remesh the solid zone and update the furnace wall location along the melting interface. Step 9. Update the initial temperature distribution throughout the solid zone by using Eq. (15). Go to the step 3.

5. Inverse heat transfer problem The goal of the inverse problem is to construct the melting rate in a desired pattern by controlling the time rate of heaters over the heater surface of smelting furnace. The unknown variables, which

q dT (q h ) = {ϕ d ,m , m = 1, … , M }

(15a)

q eT (q h ) = {ϕ e,m , m = 1, … , M }

(15b)

where ϕ d ,m and ϕ e,m are the desired and estimated heat fluxes over the design surface, respectively. Here q h is the vector of heat powers on the heater surface defined by:

q hT = {ϕ h,n , n = 1, … , N }

(16)

where ϕ h,n’s are the heater powers. The solution of inverse problem is based on the minimization of a least square objective function prescribed by:

G = [q d − q e ] [q d − q e ] T

(17)

The CGM is chosen to minimize the objective function, Eq. (17). The CGM is an iterative procedure in which at each iteration a suitable step size, β , is taken along a direction of descent, d , in order to minimize the objective function, so that:

qνh+1 = qνh − β dν

(18)

where the superscript ν is the iteration number. The direction of descent vector can be determined as a conjugation of gradient direction, ∇G , and the direction of descent from the previous iteration as follows:

dν = ∇G ν + λ dν −1

(19)

here λ is the conjugation coefficient given by:

λν =

∇G (qνh ) ∇G T (qνh ) with λ 0 = 0 ∇G (qνh−1 ) ∇G T (qνh−1 )

(20)

The gradient direction vector is determined as:

∇G (qνd ) = −2 S T [qνd − qνe ]

(21)

here S is the sensitivity matrix. The elements of the sensitivity matrix are:

Smn =

∂ϕ e,m ∂ϕ h,n

(22)

The estimated heat fluxes can be linearized with Taylor series expansion and then the minimization, with respect to step size β ν , is performed to yield the following expression for the step size:

[Sν dν ] [qνd − qνe ] T [ Sν dν ] [ Sν dν ] T

βν =

(23)

1144

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

5.1. Sensitivity problem

Table 1 Thermal properties of ablating material.

To minimize the objective function given by Eq. (17), it is required to calculate the components of sensitivity matrix defined by Eq. (22). The sensitivity problem is obtained by differentiating the direct problem, which included the radiation problem and the ablative phase change problem given by set of Eqs. (1) and (2), with respect to the nodal heat fluxes over the heater surfaces NS

∑ [δ

kj

− (1 − ε k ) Fk− j ]

kj

− Fk− j ]

j =1

NS

∑ [δ j =1

∂J j (t ) =0 ∂ϕ h,n (t )

∂J j (t ) = δ kn ∂ϕ h,n (t )

(24a)

(24b)

Parameter

Typical value

Unit

k α

0.22 9.11 × 10−8 44.7 × 108 833 273

W m.K

h Tmelt T0

m2 s J m3 K K

10. Replace ν by ν + 1 and go to Step 3. 7. Results and discussion 7.1. Verification of direct problem

Solving Eqs. (24), the components of the sensitivity matrix can be obtained by differentiating Eq. (4) with respect to heat powers over the heater surface as follows:

∂ϕ e,m εm ∂J m (t ) =− , m = 1, … , M , n = 1, … , N ∂ϕ h,n (1 − ε m ) ∂ϕ h,n (t )

(25)

The sensitivity matrix is a matrix with M rows and N columns. Obviously, Eqs. (24) and (25) are similar to Eqs. (1)–(3), and they solve in a similar fashion. Hence, by solving an equivalent radiative transfer problem, defined by Eqs. (24) and (25), the equivalent heat fluxes over the design surface, ∂ϕ e,m ∂ϕ h,n , are in fact the components of the n-th column of sensitivity matrix. Therefore, in order to complete the sensitivity matrix, the Eqs. (24) and (25) must be solved N times. 6. Overall computational procedure The computational procedure for the inverse problem at each time step is given by following expressions: Step 1. Set ν = 0 and guess a set of heater powers over the heater surface, qνh . Step 2. Solve the direct problem, and compute the conductive heat fluxes, {ψ m , m = 1, … , M } , and estimated radiative heat fluxes, qνe . Step 3. Knowing the conductive heat fluxes over the elements of design surface, and the desired melting rate, calculate the desired radiative heat fluxes over the elements of design surface, q d , using Eq. (13). Step 4. Knowing q d and q e , calculate the objective function G given by Eq. (17). Terminate the iteration procedure if the objective function is less than a small prescribed value. Otherwise go to step 5. Step 5. Determine the sensitivity matrix by solving Eqs. (24) and (25). Step 6. Compute the gradient direction, ∇G , from Eq. (21), then compute the conjugate coefficient, λ , from Eq. (20). Step 7. Compute the direction of descent, dν , from Eq. (19). Step 8. Compute the search step size, β ν , from Eq. (23). Step 9. Compute a new set of heaters, q hv+1, over the heater surface, using Eq. (18). Step

To evaluate the performance of moving grid scheme in this paper, the numerical results are compared with the exact analytic solution of one dimensional ablating problem. The material properties, initial and boundary conditions are shown in Table 1 [29]. The comparison of ablation front location in numerical solution with exact approximate solution is shown in Fig. 4. In addition, the numerical temperature history during the ablation process is compared with the approximate exact solution of benchmark problem, as depicted in Figs. 5 and 6. As it is shown, the ablated interface location and the temperature history along the time have good agreement with the approximate exact solution. 7.2. Inverse problem The accuracy and robustness of inverse solution are examined by a test case. The geometry of the furnace (with 3.5m height and 5m width) and the concentrate bank are shown in Fig. 1. The prop-

Fig. 4. The schematic shape of one-dimensional ablation problem.

Fig. 5. Comparison of front location by exact analytic solution and present numerical solution [29].

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

1145

Table 4 Thermal properties of solid zone.

Fig. 6. Solid phase temperature at t = 0.2 s, 1 s, 2 s, 3 s, and 4 s. The solid and dashed lines indicate the exact analytical solution and the numerical solution, respectively [29].

erties of metal concentrate bank are approximated by the average properties of copper. The boundary conditions over all surfaces of furnace, including heater surface, slag surface, and ablative interface, are shown in Table 2. To study the mesh independency of numerical solution, three cases with 1098, 1420, and 2404 elements to triangulate the solid zone at the first time step are studied. Table 3 shows the maximum relative error in temperature with respect to the case with 2404 elements. It is seen that the relative errors are less than 2% for all cases. Hence, in order to decrease the numerical expenses of numerical calculations the case with 1098 elements is selected. Since the solid zone is shrinked during the time by ablation over the moving surface, for a given expansion factor [28], the mesh becomes finer over the time. Although the molten concentrate ablates and drains into the slag bath, but at each time, a thin layer of molten metal covers the moving interface. Hence, the emissivity of ablative interface is considered to be equal to the emissivity of the slag surface. The right wall and the step wall, where the copper concentrate bank is located, are adiabatic. The heaters are placed alternatively on the upper and the left walls, and the walls between two consecutive heaters are assumed to be adiabatic. The heaters on the top wall and the left wall have 0.25m and 0.175m widths, respectively. The ablated molten material is flown to the furnace bed at the bottom of the furnace, which creates a molten slag bath. Therefore, the bottom

Table 2 Furnace boundary conditions. Quantity

Value

Heat flux over the furnace walls and roof Temperature over the melting interface Temperature over the surface of slag bath Heat powers of the heaters Emissivity of furnace walls Emissivity of the slag and moving surfaces

ϕ =0 Tmelt = 1500°C Ts = 1450°C ϕ h ≡ unknown ε w = 0.9 ε s = 0.05

Table 3 Maximum relative error in temperature with respect to the case with 2404 elements for triangulation of solid zone. Case number

Relative error (%)

1098 Elements 1420 Elements

1.9% 1.4%

Parameter

Typical value

Unit

k α

0.021 0.0175

W m.K

h Tmelt T0

5 × 105 1500 450

m2 s J m3 °C °C

surface is assumed to be isothermal at melting temperature. The melting interface is the design surface which is assumed to be at constant melting temperature. The thermal properties of solid zone and the initial condition are shown in Table 4. The aim of the inverse problem is to find the set of heater powers over the heater surface of the furnace to obtain a desired uniform melting rate over the design surface. During the inverse procedure by the CGM, negative values for heater powers may be obtained. In these cases, the negative heater powers may be replaced by zero, and the iteration procedure is continued until convergence. To achieve the specified melting rate over design surface, a distinct heat flux distribution should be imposed over the design surface by the heaters. In the present test case, a uniform melting rate, 1mm s , is considered for the entire design surface, so the concentrate mass recedes monotonously. Fig. 7 shows the temperature distribution in the concentrate mass during 6 minutes ablation process in the furnace. Figs. 8 and 9 show the rate of change for powers of heaters over the upper and left walls of the furnace to achieve the uniform melting rate. Inasmuch as the view factors between the elements over the heater and design surfaces are decreased by recession of design surface, to have a desired uniform ablation rate, the power of heaters must be increased. Since the left wall has a greater view factor with respect to the design surface, the heaters on the left wall provide major part of required heat flux over the design surface, in comparison with those on the upper wall. Moreover, as seen, the power of heaters over some locations of the upper wall is decreased during the ablation process, which causes to decrease the effect of the heaters over the upper wall. In comparison with those heaters on the upper wall, the heaters on the left wall behave in the same trend, in such a way that the power of heaters do not vary significantly from one heater to other. Since the conductive heat flux over the design surface decreases, due to decrease of temperature gradient, the greater portion of received radiative heat flux consumes to ablate the bank surface. Fig. 10 shows the time variation of required radiative heat fluxes over the design surface. As seen, the required radiative heat flux over the design surface reaches an asymptote. In the current investigation, the root mean square error (RMSE) is used to show the efficacy and the performance of the inverse method. RMSE is defined as:

RMSE =

s −s ⎞ 1 M⎛ ∑ 100 × d,msd,m e,m ⎟⎠ M m=1⎜⎝

2

(26)

Fig. 11 depicts the RMSE between the desired melting rate and the achieved melting rate along the design surface. As seen, RMSE is less than 1% for major part of the melting process, which shows an appropriate agreement between the desired values and the estimated values of melting rate. Fig. 12 shows the volume change of concentrate during the melting process. Because of constant ablating rate and reduction of design surface during the melting process, the ablated volume decreases by the recess of design surface, but the volume rate of ablation (the slop of the curve) is constant during the process.

1146

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

Fig. 7. Temperature contours and recession at (a) 60 s , (b) 120 s , (c) 240 s , (d) 360 s .

Fig. 8. The rate of change for the heaters over the upper wall.

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

1147

Fig. 9. The rate of change for the heaters over the left wall.

8. Conclusion

design surfaces decreases, in such a way that the radiative heater powers are increased.

By controlling the melting process in metallurgical smelting furnace, not only an optimized procedure may be achieved, but also the energy consumption by the heaters can be controlled. In the present paper, we show that the desired melting rate may be achieved by an inverse approach. To create the desired melting rate on the moving design surface, the radiative power of heaters should be controlled to create an appropriate heat flux distribution over the design surface. It was shown that the temperature gradient through the melting surface is decreased and reaches an asymptote. By decreasing the temperature gradient, the conductive heat flux toward the concentrate bank diminishes. However, as the ablative surface recesses, the view factors between the heater and

Nomenclature

c d F G h J k M N

Specific heat [ J kg .K ] Direction of decent vector View factor Objective function Latent heat of melting [ J m3 ] Radiosity [ W m2 ] Conductivity [ W m.K ] Number of surface elements over the design surface Number of surface elements over the heater surface

Fig. 10. Elemental radiative heat fluxes over the design surface.

1148

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

Superscripts e Element T Transpose of matrix ν Iteration number Subscripts 0 Initial at t = 0 d Desired, design e Estimated h Heater m Design surface element n Heater surface element s Slag w Wall

References Fig. 11. Root mean square errors between desired and estimated melting rates.

n NS s S t T x, y, z

Normal direction Total number of surface elements over the furnace walls Ablation front location [ m ] Sensitivity coefficient Time [ s] Temperature [ K ] Cartesian coordinates [ m ]

Greek symbols ε Emissivity α Thermal diffusivity [ m2 s ] γ Shape function ϕ Radiative heat flux [ W m2 ] ψ Conductive heat flux [ W m2 ]

δ ρ ω λ σ β

Kronecker delta function Density [ kg m3 ] Time parameter Conjugation coefficient Stefan–Boltzmann constant [ W m2 .K 4 ] Search step size

Fig. 12. The rate of change for ablated volume of concentrate.

[1] R. Crzymkowski, D. Slota, One-phase inverse Stefan problem solved by Adomian decomposition method, Comput. Math. Appl. 51 (2006) 33–40. [2] D. Slota, Direct and inverse one-phase Stefan problem solved by the variational iteration method, Comput. Math. Appl. 54 (2007) 1139–1146. [3] O. Tadrari, M. Lacroix, Prediction of protective banks in high temperature smelting furnaces by inverse heat transfer, Int. J. Heat Mass Transf. 49 (2006) 2180–2189. [4] M. LeBreux, M. Désilets, M. Lacroix, Fast inverse prediction of phase change banks in high temperature furnaces with a Kalman filter coupled with a recursive least-square estimator, Int. J. Heat Mass Transf. 53 (2010) 5250– 5260. [5] C.H. Huang, B.H. Chao, An inverse geometry problem in identifying irregular boundary configurations, Int. J. Heat Mass Transf. 40 (1997) 2045–2053. [6] J.V. Beck, B. Blackwell, A. Haji-Sheikh, Comparison of some inverse heat conduction methods using experimental data, Int. J. Heat Mass Transf. 39 (1996) 3649–3657. [7] P.C. Tuan, C.C. Ji, L.W. Fong, An input estimation approach to on-line twodimensional inverse heat conduction problems, Numer. Heat Tr. B-Fund 29 (1996) 345–363. [8] J.R. Howell, O.A. Ezekoye, J.C. Morales, Inverse design model for radiative heat transfer, J. Heat Transfer 122 (2000) 492–502. [9] A.G. Fedorov, V.R. Voller, R. Viskanta, Inverse optimal design of the radiant heating in materials processing and manufacturing, J. Mater. Eng. Perform. 7 (1998) 719–726. [10] K.J. Daun, J.R. Howell, D.P. Morton, Optimization of Radiant Enclosures through Non-linear Programming, in: Proc. 12th Int. Heat Transfer Conf, Grenoble, France, 2002. [11] S.M. Hosseini Sarvari, S.H. Mansouri, J.R. Howell, Inverse design of threedimensional enclosures with transparent and absorbing-emitting media using an optimization technique, Int. J. Heat Mass Transf. 3 (2003) 149–162. [12] S.M. Hosseini Sarvari, S.H. Mansouri, Inverse design for radiative heat source in two-dimensional participating media, Numer. Heat Transf. 46 (2004) 283–300. [13] S.M. Hosseini Sarvari, S.H. Mansouri, J.R. Howell, Inverse boundary design radiation problem in absorbing-emitting media with irregular geometry, Numer. Heat Transf. 43 (2003) 565–584. [14] S.M. Hosseini Sarvari, J.R. Howell, Inverse boundary design conduction-radiation problem in irregular two-dimensional domains, Numer. Heat Transf. 44 (2003) 209–224. [15] S.M. Hosseini Sarvari, Inverse determination of heat source distribution in conductive–radiative media with irregular geometry, J. Quant. Spectrosc. Radiat. Transf. 93 (2005) 383–395. [16] S.M. Hosseini Sarvari, Optimal geometry design of radiative enclosures using the genetic algorithm, Numer. Heat Transf. 52 (2007) 127–143. [17] S. Payan, S.M. Hosseini Sarvari, H. Ajam, Inverse boundary design of square enclosures with natural convection, Int. J. Therm. Sci. 48 (2009) 682–690. [18] A.A. Shokouhi, S. Payan, A. Shokouhi, S.M. Hosseini Sarvari, Inverse boundary design problem of turbulent forced convection between parallel plates with surface radiation exchange, Heat Transfer Eng. 36 (2015) 488–497. [19] S. Payan, S.M. Hosseini Sarvari, A. Behzadmehr, Reconstruction of temperature distribution in the combustion region of a non-gray medium, Numer. Heat Tr. A-Appl. 68 (2015) 908–924. [20] S. Wang, S.C. Lin, Y.C. Yang, Geometry estimation for the inner surface in a furnace wall made of functionally graded materials, Int. J. Heat Mass Transf. 67 (2015) 1–7. [21] V. Petrushevsky, S. Cohen, Nonlinear inverse heat conduction with a moving boundary: heat flux and surface recession estimation, J. Heat Transfer 121 (1999) 708–711. [22] A.P.D. Oliveira, H.R.B. Orlande, Estimation of the heat flux at the surface of ablating materials by using temperature and surface position measurements, Inverse Probl. Sci. Eng. 12 (2004) 563–577. [23] A. Hakkaki-Fard, F. Kowsary, Heat flux estimation in a charring ablator, Numer. Heat Tr. A-Appl. 53 (2007) 543–560.

H. Farzan et al./Applied Thermal Engineering 98 (2016) 1140–1149

[24] H. Molavi, A. Hakkaki-Fard, M. Molavi, R.K. Rahmani, A. Ayasoufi, S. Noori, Estimation of boundary conditions in the presence of unknown moving boundary caused by ablation, Int. J. Heat Mass Transf. 54 (2011) 1030–1038. [25] M. LeBreux, M. Désilets, M. Lacroix, An unscented Kalman filter inverse heat transfer method for the prediction of the ledge thickness inside hightemperature metallurgical reactors, Int. J. Heat Mass Transf. 57 (2013) 265–273. [26] K.S. Jhajj, S.R. Slezak, K.J. Daun, Inferring the specific heat of an ultra high strength steel during the heating stage of hot forming die quenching, through inverse analysis, Appl. Therm. Eng. 83 (2015) 98–107.

1149

[27] J.R. Howell, R. Siegel, M.P. Menguc, Thermal Radiation Heat Transfer, Taylor & Francis Group, Florida, 2011. [28] S. Lo, A new mesh generation scheme for arbitrary planar domains, Int. J. Numer. Methods Eng. 21 (1985) 1403–1426. [29] S.L. Mitchell, T.G. Myers, Heat balance integral method for one-dimensional finite ablation, J. Thermophys. Heat Transfer 22 (2008) 508–514.