A general method for predicting the bank thickness of a smelting furnace with phase change

Applied Thermal Engineering 162 (2019) 114219

Contents lists available at ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

Research Paper

A general method for predicting the bank thickness of a smelting furnace with phase change ⁎

Bowen Zhanga, Jie Meia, Chunyun Zhanga, Miao Cuia, , Xiao-wei Gaoa, Yuwen Zhangb,

T

⁎

a

State Key Laboratory of Structural Analysis for Industrial Equipment, Key Laboratory of Advanced Technology for Aerospace Vehicles, Dalian University of Technology, Dalian 116024, China b Department of Mechanical and Aerospace Engineering, University of Missouri, Columbia, MO 65211, USA

H I GH L IG H T S

thickness of a smelting furnace involving phase change is predicted. • Bank transient nonlinear inverse model with phase change is developed for 1st time. • 2D and ABAQUS with UEL are coupled to ensure the good inversion performances. • CVDM • The prediction is accurate, efficient, stable and robust.

A R T I C LE I N FO

A B S T R A C T

Keywords: Inverse heat conduction Solid-liquid phase change Melting furnace

A new methodology is proposed for predicting the bank thickness covering and protecting the refractory brick walls of the smelting furnaces. The inverse method predicts the bank thickness changing with both time and coordinates, by using a two-dimensional physical model of the furnace that is more general and challenging than the previous studies. Moreover, the thermal conductivities of the phase change material are identified simultaneously because they are poorly known, which is more challenging and practical. The inverse problem is solved by the Levenberg-Marquardt Method (LMM). The direct problem of the smelting furnace and the sensitivity matrix coefficients of the LMM are calculated by ABAQUS, using a complex user-defined element that is set up based on the complex-variable-differentiation method (CVDM) and the user element subroutine (UEL). Finally, numerical examples are given to examine the performances of the approach for predicting the bank thickness covering the refractory brick walls of the smelting furnaces.

1. Introduction Smelting furnaces are commonly employed to melt materials such as steel, cooper, and nickel calcine, which require high powers and temperatures. Fig. 1 shows the schematic of a typical smelting furnace [1], in which the materials are melted by the large heat provided by the high voltage electrodes (only one electrode is shown here) and solidified by the cold brick wall [2]. The solid zone near the brick wall is called a “bank” that could protect the inner surface of the brick wall from the highly corrosive molten materials [1,3]. The thickness of the bank affects the life of the brick wall, which is difficult to be probed due to the hostile conditions [1]. Moreover, the thickness would change with both time and coordinates, which makes the problem more complex and difficult. An ingenious method was proposed to solve this problem to

⁎

determine the bank thickness, in which the inverse heat transfer problem was involved [1–3]. The main idea was installing sensors in the brick wall to obtain temperature or flux information, and then recovered the bank thickness using an inverse heat transfer algorithm based on the temperature or flux information. Solving an inverse heat conduction problem is composed of two main parts: solutions to both the direct heat conduction and the inverse problems, and the connection/couple of the solutions. For heat conduction problems, many numerical methods could be used, such as the finite difference method (FDM) [4,5], the finite element method (FEM) [6,7], the boundary element method (BEM) [8,9], the meshless method (MLM) [10,11], and the finite volume method (FVM) [1–3,12]. Every method has its advantages and disadvantages, and the FEM is a powerful method in application and there is available commercial software, such as ANSYS, NASTRAN, ABAQUS, and so on. Among these software

Corresponding authors. E-mail addresses: [email protected] (M. Cui), [email protected] (Y. Zhang).

https://doi.org/10.1016/j.applthermaleng.2019.114219 Received 20 May 2019; Received in revised form 17 July 2019; Accepted 6 August 2019 Available online 07 August 2019 1359-4311/ © 2019 Elsevier Ltd. All rights reserved.

Applied Thermal Engineering 162 (2019) 114219

B. Zhang, et al.

Nomenclature

Greek

a b c E(x,y,t) Erel

δ ζ γ λ ξ ρ μ

e H h J M N n p Q q s T t x y

real part imaginary part mass specific heat, J kg−1 K−1 bank thickness, m relative error between the identified and the exact/real values, % a small imaginary part heat of fusion, J kg−1 heat convective coefficient, W m−2 K−1 sensitivity matrix total number of measurements total number of identified parameters number of degrees of freedom iteration number heat source, W kg−1 heat flux, W m−2 inverted parameter vector temperature, K time, s x-coordinate, m y-coordinate, m

updated vector of recovered/inverted parameters random measurement error random number thermal conductivity, W m K−1 small positive number density, kg m−3 damping factor

Subscripts 0 i ∞ S L Im Re

initial time the ith component of a vector environmental temperature, K solid (PCM) liquid (PCM) imaginary part value real part value

Superscripts # * 0

measured value complex value initial guess

take a considerable amount of time. The gradient methods are more suitable because of their high accuracy and efficiency [26–33]. For the gradient method, the accurate evaluation of sensitivity matrix coefficients is a key issue. In the recent years, the authors’ group have successfully employed the complex-variable-differentiation method (CVDM) to calculate the sensitivity matrix coefficients [25–27]. In this study, a complex element is set up through user element subroutine (UEL) [35] to combine the complex-variable-differentiation method with ABAQUS. Then, the Levenberg–Marquardt method [32] is applied to solve the inverse problem. The direct heat conduction problem with phase change must be solved/called many times for solving an inverse problem, and this will take a considerable amount of time. Thus, the physical model of the smelting furnace is commonly reducible to one-dimensional [1–3,5,36], which is reasonable in some special cases. In this study, the model is two-dimensional, which is more general and challenging than the 1-D problem. In addition, the thermal conductivities of the phase change material are identified simultaneously because they are poorly known [36], which is more challenging and practical. The remainder of this paper is as follows. In Section 2, the twodimensional physical model of the smelting furnace and basic assumptions are described. The direct problem is solved by using FEM in the Section 3. Section 4 briefly reviews the Levenberg–Marquardt method, and describes a method to combine the complex-variable-differentiation method with ABAQUS using the user element subroutine (UEL). In Section 5, the efficiency, accuracy coupled with effects of initial guesses and measurement errors for predicting the bank thickness of the smelting furnace are investigated in detail. Conclusions are drawn in Section 6.

packages, the ABAQUS is best known for its nonlinear solution capability and it is very suitable for heat conduction problems with phase change occurs in a smelting furnace. Moreover, the ABAQUS provides users abundant subroutines to carry out the secondary development. Therefore, the direct problem is solved by using the FEM in the ABAQUS in the present work. To solve an inverse heat transfer problem, many methods were proposed, such as the conjugate gradient method [13–17], the multiple model adaptive inverse method [18], the genetic algorithm [19], the Levenberg–Marquardt method [20–22], the Kalman filter method [23,24] and so on. These methods could be classified into two categories overall [25]: the stochastic and the gradient ones. Since the problem under consideration in this paper is highly nonlinear, it would

2. Two-dimensional model of the smelting furnace The two-dimensional sketch of a smelting furnace is shown as Fig. 2. Here, the model is delineated as an annulus, and the physical model consists of brick and phase change material parts. The heat flux imposed on the inner face q(x,y,t) that represents the furnace heat load varies with both time and coordinates; the outer surfaces of the brick are subjected to convective boundary conditions. The air temperature

Fig. 1. Cross view of a typical smelting furnace [1]. 2

Applied Thermal Engineering 162 (2019) 114219

B. Zhang, et al.

∂T ∂T nx + λ y n y = h (Ta − T ) (outer boundary) ∂y ∂x

λx

(4)

The effective heat capacity of mushy zone can be obtained based on the basic energy balance [39]:

∫T

Tl

s

ca (T ) dT =

∫T

s

Tl + Ts 2

cs (T ) dT +

Tl

∫T +T

s

l

2

cl (T ) dT + H

(5)

where cs and cl are the heat capacity of solid and liquid, Ts and Tl are the solidus and liquid temperature of PCM, and H is the latent heat of fusion, respectively. Suppose that the cs and cl do not depend on temperature, then the effective heat capacity ca can be written as follows based on Eq. (5).

ca = Fig. 2. Schematic diagram of a smelting furnace.

(1) The temperature gradients across the x and y directions are much greater than that in the z direction. Thus, the two-dimensional model is adopted in this study. (2) The phase change problem is non-isothermal. (3) The thermal properties of these two materials are temperature independent. But for the phase change material (PCM), the thermal properties are different for the solid and liquid phase. (4) The heat transfer inside the liquid phase of the phase change material (PCM) is conduction dominated [1–3]. (5) The thermal contact resistance between these two materials is neglected.

λx

⎧ [C ] + [K ] ⎫ {T }t + Δt = {P }t + Δt + [C ] {T }t ⎨ ⎬ Δt ⎩ Δt ⎭

(6)

(7)

(8)

For each time step, the heat capacity matrix [C], conduction matrix [K], temperature vector {T}t and load vector {P}t+Δt are known, while the temperature vector {T}t+Δt at time t + Δt is unknown. Representing [C ] + [K], {T}t+Δt, and {P}t+Δt + [C ] {T}t with K¯ , T¯ and P¯ , respectively, Δt Δt it follows that

¯ = P¯ KT

(9)

In this study, the FEM code is set up using ABAQUS user element subroutine (UEL) and the system equation Eq. (9) is solved by ABAQUS/Standard solver. The following examples will validate the accuracy of this numerical model. Consider a typical one-dimensional solidification problem which is reported in Voller and Swaminathan [42]. The geometry and boundary conditions are shown in Fig. 3, and the thermal properties of the alloy are shown in Table 1. The initial temperature is 969 K, and the left side is specified with temperature of 573 K. The mesh size is 0.001 m and the time step is set to 5 s. The movement of both the solidus and the liquidus fronts compared with Voller and Swaminathan [42] are shown in Fig. 4. It can be seen from Fig. 4 that the results calculated by this numerical model are in excellent agreement with those of Voller and Swaminathan [42], which proves the effectiveness and the accuracy of

(1) (2)

∂T ∂T nx + λ y n y = q (inner boundary) ∂x ∂y

T > Tl

where [C] is the heat capacity matrix, {∂ T/∂ t} is a vector that contains the derivatives of temperature with respect to time, [K] is the heat conduction matrix, {T} is the temperature vector, and {P} is the vector of thermal load. Using the Euler backward time difference scheme [37], one could obtain:

Generally, there are two kinds of methods to solve the heat conduction problem with phase change, front-tracking methods and fixeddomain methods [37]. But for this problem, the front-tracking methods is not suitable, because the methods cannot deal with the phase change problems with a finite freezing range [38]. In this study, the problem is solved by effective heat capacity method, which is a widely used fixeddomain method and suitable for the problems with phase change [38]. Based on the assumptions illustrated in Section 2, the governing equation for the effective heat capacity method can be written as follows [39]:

T = T¯ (initial temperature)

Ts ⩽ T ⩽ Tl

[C ]{∂T / ∂t } + [K ]{T } = {P }

3. The direct problem

∂T ∂ ∂T ∂ ∂T = (λ x )+ (λ y ) ∂t ∂x ∂x ∂y ∂y

T < Ts H Tl − Ts

In severe latent heat cases, this method may result in numerical instabilities, as there is a discontinuous step jump of ca, which will cause the difficulties in evaluating the capacitance matrix [40]. Bonacina et al. [39] suggested that some numerical smooth approximations could be used to avoid such instabilities and to improve the calculated accuracy. In this study, the ca is modified to a secant term during the early iterations of the solution to a time step. Then, the FEM form of this problem can be written as follows by using the variational principle [41]:

and the convective heat transfer coefficient are denoted by T∞ and h, respectively. In the working process of the smelting furnace, the phase change material part can be divided to three zones: solid layer, mushy layer, and liquid layer. The thickness of the solid varies with both time and coordinates. The actual phase change process is complicated; Thus, the mathematical model is simplified based on the following assumptions:

ρ (T ) ca

cs ⎧ ⎪ cs + cl + ⎨ 2 ⎪ cl ⎩

(3)

Fig. 3. The geometry and boundary conditions of a typical 1D solidification problem. 3

Applied Thermal Engineering 162 (2019) 114219

B. Zhang, et al.

Table 1 The thermal properties of the aluminum-copper alloy [42]. Parameter

Value

Unit

cp, solid λsolid ρ cp, liquid λliquid H TS TL TF

900 200 2800 1100 90 3.9 × 105 821 919 933

J kg−1 K−1 W m K−1 kg m−3 J kg−1 K−1 W/m K J kg−1 K K K

Fig. 5. A quarter FEM model of the smelting furnace.

ensure that this grid and time step are reasonable, and the results are shown in Fig. 6. Fig. 7 shows the computed temperatures at points E, F and G. It can be seen that the results calculated with UEL are in excellent agreement with those calculated by traditional ABAQUS, which further validates the accuracy of the proposed numerical model. Fig. 4. The movement of both the solidus and the liquidus fronts.

4. The inverse problem Different from the direct problem, some parameters are unknown and need to be inverted in the inverse problem, which are depicted in Fig. 8. For this problem, the parameters which need to be inverted are the value of a, b and c in Eq. (10), and the thermal conductivity of the PCM, λsolid and λliquid. Once the value of these parameters are inverted, the thickness E(x,y,t) can be obtained by solving the direct problem. The inverse problem can be constructed as a problem of minimization of the following objective function:

Table 2 The thermal properties of the smelting furnace. Parameter

Value

Unit

ρbrick cbrick λbrick ρ cp, solid cp, liquid λsolid λliquid H TS TL

2600 875 16.8 2100 1600 1600 2 13 4.2 × 105 1115 1155

kg m−3 J kg−1 K−1 W m K−1 kg m−3 J kg−1 K−1 J kg−1 K−1 W m K−1 W m K−1 J kg−1 K K

F (s1, s2, ...,sN ) =

4

3

2πt ) tmax 4

M

∑ (Ti# − Ti (s1, s2, ...,sN ))2 i=1

(11)

where M is the number of measured temperatures, N is the total number of the inverted parameters, The vector s = (s1, s2,…, sN) is made up of the inverted parameters. si is the ith inverted parameter. Ti# and Ti are the measured and calculated temperatures, respectively, i = 1, 2,…,M. The inverted parameters vector can be updated by:

this method. Next, taking a smelting furnace problem as the direct problem for this study, the geometry and boundary conditions are shown in Fig. 2. For the direct problem, all the parameters and conditions are known. The initial temperature is 1100 K, the air temperature and the convective heat transfer coefficients are 300 K and 15 W m−2 K−1, respectively. And the time is 5 × 104 s, the thickness of the PCM and brick wall are 0.12 m and 0.1 m, respectively. The thermal properties of the direct problem are shown in Table 2. Finally, the flux imposed on the inner face is given by

q (x , y, t ) = a·cos2 (10x − 2) + (b·sin(10y + 2) + c )·sin2 (

1 M

skp + 1 = skp + δ p

(12)

where p is the iteration number, k = 1, 2, …, N, and δ is determined by the following equation [32]:

[J T J + μ diag (J T J )] δ = J T [Ti# − Ti (s )]

(13)

where J is sensitivity coefficients matrix as shown in Eq. (14), μ is the damping factor which is adjusted at each iteration, and diag represents diagonal elements.

(10)

∂T (s )

⎡ ∂1s ⎢ 1 ⎢ ∂T2 (s) J = ⎢ ∂s1 ⎢ ⋮ ⎢ ∂TM (s) ⎢ ⎣ ∂s1

−2

where a, b and c are 1.2 × 10 , 5 × 10 and 1.4 × 10 W·m , respectively. Fig. 5 shows a quarter FEM model of the smelting furnace and the positions of points A, B and C. The total number of the nodes and elements of the whole FEM model are 5524 and 5540. The time step is fixed to 200 s. Several numerical simulations were performed to 4

∂T1 (s ) ∂s2

⋯

∂T1 (s ) ∂sN

∂T2 (s ) ∂s2

⋯

∂T2 (s ) ∂sN

⋮ ∂TM (s ) ∂s2

⎤ ⎥ ⎥ ⎥ ⋮ ⋮ ⎥ ⎥ ⋯ ∂TM (s) ⎥ ∂sN ⎦

(14)

Applied Thermal Engineering 162 (2019) 114219

B. Zhang, et al.

Fig. 6. The bank thickness of the 270° cross section in different mesh sizes and time steps: (a) mesh size; (b) time step.

Fig. 7. Computed temperatures at points E, F and G. Fig. 9. Sketch of the complex variable element.

Fig. 8. Schematic of the 2D phase-change inverse problem.

The determination of matrix J is very important in the computational procure of the inverse problem. There are many methods to compute this matrix. In this study, the matrix is computed by the complex-variable-differentiation method (CVDM) [34]. In CVDM, the variable X of a real function f (X) is replaced by the complex variable X + ie, with the imaginary part e being very small. The function f (X + ie) can be expanded in a Taylor series as:

f (X + ie ) = f (X ) + ief ′ (X ) −

e2 f ″ (X ) + o (e3) 2

(15) Fig. 10. Convergence curves of the inverted parameters.

Since e is very small, the derivative of f (X) can be written as:

f ′ (X ) =

Im(f (X + ie )) e

coefficients, either for linear or nonlinear inverse problems [25–27]. In the present study, the CVDM is combined with ABAQUS by UEL, which can compute the sensitivity coefficients matrix and Ti with high

(16)

The CVDM is with high accuracy for calculating sensitivity 5

Applied Thermal Engineering 162 (2019) 114219

B. Zhang, et al.

Fig. 11. The bank thickness of the smelting furnace at different time: (a) 12600 s; (b) 25000 s; (c) 37600 s; (d) 50000 s. Table 3 Effects of initial guessed values on the results. Initial guessed values a Test Test Test Test

1 2 3 4

1500 3000 18,000 36,000

b 625 1250 7500 15,000

c 1750 3500 21,000 42,000

λsolid 0.25 0.5 3 6

λliquid 1.625 3.25 19.5 39

Table 4 The inverted results and relative errors with different measurement errors. Total iteration number 10 11 9 10

Objective function

Inverted parameters

−7

a b c λsolid λliquid

2.7837 × 10 2.4152 × 10−7 2.3404 × 10−7 2.8706 × 10−7

ζ = 1%

ζ = 3%

ζ = 5%

results

Erel

results

Erel

results

Erel

11,951 4951.4 14,058 2.011 12.7

0.41% 0.97% 0.41% 0.55% 2.31%

11,861 4864.9 14,212 2.039 12.12

1.16% 2.70% 1.51% 1.95% 6.77%

11,767 4782.1 14,372 2.067 11.54

1.94% 4.36% 2.66% 3.35% 11.23%

Fig. 12. Convergence of the parameters based on different initial guessed values: (a) 0.125 s_exact; (b) 0.25 s_exact; (c) 1.5 s_exact; (d) 3.0 s_exact. 6

Applied Thermal Engineering 162 (2019) 114219

B. Zhang, et al.

temperature perturbation for the imaginary node. A complex number a − b⎤ a + bi can be represented as ⎡ , then, Eq. (17) could be expressed ⎣b a ⎦ as follows.

⎡ K¯ Re − K¯ Im ⎤ ⎛ TRe ⎞ = ⎛⎜ P¯Re ⎞⎟ ⎢ K¯ Im K¯ Re ⎥ ⎝TIm ⎠ ¯ ⎦ ⎝ PIm ⎠ ⎣ ⎜

⎟

(20)

Finally, the sensitivity coefficients and Ti can be obtained by computing Eq. (20) using ABAQUS/Standard solver. The iteration is stopped until the objective function or the difference between Fp+1 and Fp is within a specified tolerance [22].

F (s1, s2, …, sN ) ⩽ ξ or |F p + 1 − F p| ⩽ ξ

(21)

The inversion procedure is summarized as follows. Step 1: Calculate the sensitivity coefficients and Ti by ABAQUS/ Standard using the complex element based on the guessed values of s0. Step 2: Determine the objective function by Eq. (11). Step 3: Check the convergence criterion Eq. (21). Terminate the iteration if the convergence criterion is achieved. Otherwise, continue the following procedure. Step 4: Obtain δ by solving Eq. (13). Step 5: Update the initial guess s0, and update p by p + 1. Then, return to Step 1.

Fig. 13. The inverted bank thickness of the 270° cross section with different random measurement errors.

5. Inversion of the bank thickness of the smelting furnace For this problem, the thickness E (x, y, t) can be obtained by solving the direct problem based on the inverted parameters which are described in Section 4, and everything else is kept the same as in the direct problem. As shown in Fig. 2, four sensors (A, B, C and D) are set on the outer surfaces of the brick, and the measurement is implemented every 1000 s. Thus, there are 200 measurements. Fig. 10 shows the convergence curves of the inverted parameters, for which the initial guessed values are set to a half of their actual values, and it can be seen that the convergence could be achieved. The inverted parameters become stable after 4 iterations, and reach to their actual values after 11 iterations, which shows the high efficiency of the identification. Then, the bank thickness E (x, y, t) can be predicted by solving the direct problem. Fig. 11 shows the bank thickness of the smelting furnace at different times. To investigate the effects of initial guessed values, five sets of initial guessed values are tested. Table 3 lists the inversion results and Fig. 12 shows the convergence with different initial guessed parameters. From Table 3 and Fig. 12, it can be seen that convergence can be achieved for each set of initial guessed values, and every parameter reached to its real value, which shows initial guessed values have weak effects on inversion results. In the above analysis, the temperature measurements are exact and without any errors. However, there must be measurement errors in the actual applications. In order to investigate the effect of measurement errors on inversion results quantitatively, the relative errors are defined as

Fig. 14. The inverted bank thickness at 30000 s with different measurement errors: (a) 0%; (b) 1%; (c) 3%; (d) 5%.

accuracy and efficiency. To use the CVDM in ABAQUS, the traditional FEM in ABAQUS should be transformed into complex. Then, Eq. (9) could be written as:

K¯ *T * = P¯ *

(17)

where the complex temperature vector can be written as follows

T * = TRe + TIm·i

(18)

Erel =

In Eq. (17), K¯ * and P¯ * have the same form with T*. Then the sensitivity coefficients can be computed by the following equation:

∂TRe T = Im ∂s e

sestimasted − sexact × 100% sexact

(22)

where sexact and sestimated represent the exact and the recovered values, respectively. A random error term is added to the exact temperature to account for the random measurement error [32].

(19)

ABAQUS is a real-valued FEM, which does not have a provision for complex operations. Therefore, a complex user element is set up to implement the complex variable FEM within ABAQUS by UEL [35]. The complex element is set up by the following method described in [43]. As shown in Fig. 9, two sets of nodes (one real set and one imaginary set) are defined for each element. For this element, every node has one degree of freedom. TRe is temperature for the real node and TIm means

Tmeasured = Texact (1 +

γ ζ) 2.576

(23)

where ζ is the random measurement error, γ is a random number between −1 and 1. Three random measurement errors are considered, i.e., ζ = 1%, ζ = 3% and ζ = 5%, respectively and the random numbers applied to the exact measurement are the same for each value of ζ. The inverted 7

Applied Thermal Engineering 162 (2019) 114219

B. Zhang, et al.

results are listed in Table 4. Fig. 13 shows the bank thickness of 270° cross section, and Fig. 14 shows the bank thickness of whole model at 30000 s. From these results showed in Table 4 and Fig. 14, it can be seen that identification error increases with the increase of the measurement error, as expected. Moreover, it can be found that the inverted bank thickness extremely close to the direct, which means that the bank thickness can be predicted with high accuracy, even with certain measurement errors. The numerical example verifies the robustness of this method, and this method can be applied to the actual application.

disc in a disc brake system, Appl. Therm. Eng. 31 (2011) 2439–2448. [14] C.H. Huang, C.T. Lee, An inverse problem to estimate simultaneously six internal heat fluxes for a square combustion chamber, Int. J. Therm. Sci. 88 (2015) 59–76. [15] C.S. Liu, C.L. Kuo, W.S. Jhao, The multiple-scale polynomial Trefftz method for solving inverse heat conduction problems, Int. J. Heat Mass Transf. 95 (2016) 936–943. [16] H. Rahideh, P. Malekzadeh, M.R. Golbahar Haghighi, M. Vaghefi, Two-dimensional inverse transient heat conduction analysis of laminated functionally graded circular plates, Numer. Heat Transf. Part A: Appl. 62 (2012) 992–1014. [17] M.R. Golbahar Haghighi, P. Malekzadeh, M. Afshari, Inverse internal pressure estimation of functionally graded cylindrical shells under thermal environment, Acta Mech. 225 (2014) 3377–3393. [18] C. Lv, G.J. Wang, H. Chen, S.B. Wan, Estimation of the moving heat source intensity using the multiple model adaptive inverse method, Int. J. Therm. Sci. 138 (2019) 576–585. [19] S. Mohasseb, M. Moradi, T. Sokhansefat, F. Kowsari, A. Kasaeian, O. Mahian, A novel approach to solve inverse heat conduction problems: coupling scaled boundary finite element method to a hybrid optimization algorithm, Eng. Anal. Bound. Elem. 84 (2017) 206–212. [20] F.B. Liu, A hybrid method for the inverse heat transfer of estimating fluid thermal conductivity and heat capacity, Int. J. Therm. Sci. 50 (2011) 718–724. [21] E. Garcia, D. Meresse, I. Pombo, S. Harmand, J.A. Sanchez, Identification of heat partition in grinding related to process parameters, using the inverse heat flux conduction model, Appl. Therm. Eng. 66 (2014) 122–130. [22] T. Xie, Y.L. He, Z.X. Tong, W.X. Yan, X.Q. Xie, An inverse analysis to estimate the endothermic reaction parameters and physical properties of aerogel insulating material, Appl. Therm. Eng. 87 (2015) 214–224. [23] H.M. Park, An efficient Karhunen-Loeve Galerkin method for sequential solution of inverse problems, Int. J. Therm. Sci. 73 (2013) 8–17. [24] S. Wen, H. Qi, Y.T. Ren, J.P. Sun, L.M. Ruan, Solution of inverses radiation-conduction problems using a Kalman filter coupled with the recursive least-square estimator, Int. J. Heat Mass Transf. 111 (2017) 582–592. [25] M. Cui, W.W. Duan, X.W. Gao, A new inverse analysis method based on a relaxation factor optimization technique for solving transient nonlinear inverse heat conduction problems, Int. J. Heat Mass Transf. 90 (2015) 491–498. [26] M. Cui, X.W. Gao, J.B. Zhang, A new approach for the estimation of temperature dependent thermal properties by solving transient inverse heat conduction problems, Int. J. Therm. Sci. 58 (2012) 113–119. [27] X.C. Yu, Y.G. Bai, M. Cui, X.W. Gao, Inverse analysis of thermal conductivities in transient non-homogeneous and non-linear heat conductions using BEM based on complex variable differentiation method, Sci. China Phys. Mech. 56 (2013) 966–973. [28] M. Cui, Q.H. Zhu, X.W. Gao, A modified conjugate gradient method for transient nonlinear inverse heat conduction problems: a case study for identifying temperature-dependent thermal conductivities, ASME J. Heat Transf. 136 (2014) 091301. [29] W.W. Han, H.B. Chen, T. Lu, Estimation of the time-dependent convective boundary condition in a horizontal pipe with thermal stratification based on inverse heat conduction problem, Int. J. Heat Mass Transf. 132 (2019) 723–730. [30] M. Cui, N. Li, Y.F. Liu, X.W. Gao, Robust inverse approach for two-dimensional transient nonlinear heat conduction problems, J. Thermophys. Heat Transf. 29 (2015) 253–262. [31] M. Cui, K. Yang, X.L. Xu, S.D. Wang, X.W. Gao, A modified Levenberg-Marquardt algorithm for simultaneous estimation of multi-parameters of boundary heat flux by solving transient nonlinear inverse heat conduction problems, Int. J. Heat Mass Transf. 97 (2016) 908–916. [32] M. Cui, Y. Zhao, B.B. Xu, X.W. Gao, A new approach for determining damping factors in Levenberg-Marquardt algorithm for solving an inverse heat conduction problem, Int. J. Heat Mass Transf. 107 (2017) 747–754. [33] M. Cui, L. Zhou, J. Mei, B.W. Zhang, Estimation of slab surface radiative emissivities by solving an inverse coupled conduction, convection and radiation problem, Numer. Heat Transf. Part A: Appl. 72 (2017) 765–779. [34] J.N. Lyness, C.B. Moler, Numerical differentiation of analytic functions, SIAM J. Numer. Anal. 4 (1967) 202–210. [35] R. Fielder, A. Montoya, H. Millwater, P. Golden, Residual stress sensitivity analysis using a complex variable finite element method, Int. J. Mech. Sci. 133 (2017) 112–120. [36] M. Hafid, M. Lacroix, An inverse heat transfer method for predicting the thermal characteristics of a molten material reactor, Appl. Therm. Eng. 108 (2016) 140–149. [37] J. Crank, Free and Moving Boundary Problems, Calverdon Press, Oxford, 1984. [38] R.W. Lewis, K. Ravindran, Finite element simulation of metal casting, Int. J. Numer. Meth. Eng. 47 (2000) 29–59. [39] C. Bonacina, G. Comini, A. Fasano, M. Primicerio, Numerical solution of phasechange problems, Int. J. Heat Mass Transf. 16 (10) (1973) 1825–1832. [40] R.W. Lewis, K. Morgan, H.R. Thomas, K.N. Seetharamu, The Finite Element Method in Heat Transfer Analysis, John Wiley & Sons, 1996. [41] X.C. Wang, Finite Element Method, first ed., Tsinghua University Press, 2003. [42] V.R. Voller, C.R. Swaminathan, General source-based method for solidification phase change, Numer. Heat Transf. 19 (1991) 175–189. [43] B.W. Zhang, J. Mei, C.Y. Zhang, H.F. Peng, M. Cui, X.W. Gao, A new method for solving two-dimensional inverse heat conduction problems for inverting hybrid multi parameters, in: Proceedings of the CSEE (in press).

6. Conclusion In this paper, the bank thickness of a smelting furnace is predicted by solving a two-dimensional transient nonlinear inverse heat conduction problem. It is the first time to introduce the CVDM into the inverse problem with phase change. The following conclusions can be drawn: (1) The newly developed complex-variable element is accurate for solving the direct problem with phase change, which provides basis for predicting the bank thickness of the smelting furnace. Moreover, the sensitivity coefficients can be also calculated with high accuracy through this element. (2) The present methodology has good accuracy and stability, and initial guessed values have weak effects on inverted bank thickness of a smelting furnace. (3) The numerical example is given to verify the robustness of this method, and the results show that the bank thickness can be predicted with high accuracy, even with certain measurement errors. Acknowledgments Financial support of this work by the National Nature Science Foundation of China (51576026) and the China Postdoctoral Science Foundation (2016M601305) is gratefully acknowledged. References [1] O. Tadrari, M. Lacroix, Prediction of protective banks in high-temperature smelting furnaces by inverse heat transfer, Int. J. Heat Mass Transf. 49 (13–14) (2006) 2180–2189. [2] M. Hafid, M. Lacroix, Inverse heat transfer prediction of the state of the brick wall of a melting furnace, Appl. Therm. Eng. 110 (2017) 265–274. [3] M. Hafid, L. Marcel, Inverse method for simultaneously estimating multi-parameters of heat flux and of temperature-dependent thermal conductivities inside melting furnaces, Appl. Therm. Eng. 141 (2018) 981–989. [4] Y. Gu, L. Wang, W. Chen, C.Z. Zhang, X.Q. He, Application of the meshless generalized finite difference method to inverse heat source problems, Int. J. Heat Mass Transf. 108 (2017) 721–729. [5] M. LeBreux, M. Désilets, M. Lacroix, An unscented Kalman filter inverse heat transfer method for the prediction of the ledge thickness internal high-temperature metallurgical reactors, Int. J. Heat Mass Transf. 57 (1) (2013) 265–273. [6] K. Veerabhadrappa, G. Vinayakaraddy, K.N. Seetharamu, P.G. Hegde, V. Krishna, Transient behavior of three-fluid exchanger with three thermal communications under step change in inlet temperature of fluids using finite element method, Appl. Therm. Eng. 108 (2016) 1390–1400. [7] M. Cui, J. Mei, B.W. Zhang, B.B. Xu, L. Zhou, Y.W. Zhang, Inverse identification of boundary conditions in a scramjet combustor with a regenerative cooling system, Appl. Therm. Eng. 134 (2018) 555–563. [8] E. Divo, A.J. Kassab, Boundary Element Method for Heat Conduction: With Applications in Non-Homogeneous Media, WIT Press, Southampton, 2003. [9] N. Iwona, A.J. Nowak, L.C. Wrobel, Identification of phase change fronts by Bezier splines and BEM, Int. J. Therm. Sci. 41 (6) (2002) 492–499. [10] N.P. Karagiannakis, G.C. Bourantas, A.N. Kalarakis, E.D. Skouras, V.N. Burganos, Transient thermal conduction with variable conductivity using the Meshless Local Petrov-Galerkin method, Appl. Math. Comput. 272 (2016) 676–686. [11] H.T. Yang, Y.Q. He, Solving heat transfer problems with phase change via smoothed effective heat capacity and element-free Galerkin methods, Int. Commun. Heat Mass Transf. 37 (4) (2010) 385–392. [12] S. Han, Finite volume solution of two-step hyperbolic conduction in casting sand, Int. J. Heat Mass Transf. 93 (2016) 1116–1123. [13] Y.C. Yang, W.L. Chen, A nonlinear inverse problem in estimating the heat flux of the

8