COMPARISON OF METHODS FOR SOLVING INVERSE PROBLEMS TO ESTIMATE THE WEAR LINE IN A BLAST FURNACE HEARTH

A key part to the iron-making blast furnace is the hearth. How well it functions significantly affects the performance of the whole furnace and, ultimately, the hot metal quality. The hostile conditions of the hearth lead to the erosion of the refractory lining; the extent of the erosion is known as the wear line. The wear line is usually the main factor limiting the length of a furnace’s campaign. To measure the wear line while the furnace is in regular operation is of course impossible, yet such information would be useful in operation planning. Estimating the wear line in a blast furnace hearth is as close as we can get to obtaining this information, and this paper presents a model for doing just that. The procedure is based on the solution of a nonlinear inverse heat transfer problem. The observations are temperature measurements at points inside the blast furnace hearth and the unknown is the location of the 1150 oC isotherm. This work considers: (i) the finite elements to solve the direct problem, (ii) the fixed-point iteration method to solve the non-linearity of the direct problem due the thermal conductivity depending on the temperature, (iii) three methods to solve the inverse problem. We estimate, using all three methods, the wear line of a hearth with characteristics similar to those of Blast Furnace 3 at steelmaker Arcelor Mittal Tubar ̃ o1. The methods, the performance and accuracy of which this paper compares, are the Levenberg-Marquardt, the iteratively regularized Gauss-Newton, and the conjugate gradient. We validate the solution using simulated measurements with different noise levels.


Introduction
When we determine unknown causes by observing their effects, we are solving an inverse problem.In inverse heat transfer problems (IHTP), some of the unknown causes are initial conditions, thermo-physical properties, and heat flux.Their effects are usually temperature measurements at points along a domain [18].
IHTP belongs to a class of mathematical problems known as ill-posed problems.This means that small perturbations in the observed data can lead to large errors in the solution.Such problems are in direct contrast to Hadamard's well-posed problems [10].Well-posed problems require three conditions: (i) a solution must exist; (ii) the solution must be unique; (iii) given small changes to input data, the solution must be stable (i.e., stability condition).
Such an inverse heat transfer problem is solved in this paper where the wear on a blast furnace hearth, a key part to iron-making blast furnaces, is estimated.The hostile conditions inside a hearth erode the refractory lining, which is known as the wear line.Such erosion is usually the main factor limiting the length of the furnace campaign.While impossible to obtain during regular operation, direct measures of the lining's erosion would be useful in the planning of strategic operations.If the hearth relining is scheduled too late, there is a risk of an outbreak with potentially catastrophic consequences.On the other hand, a too early relining of the hearth causes unnecessary loss of production and capital [4].This paper thus presents a model to better estimate the wear line of a blast furnace hearth.
The lowest temperature at which iron can exist in liquid form, 1150 o C, is known as the eutectic temperature of carbon-saturated iron [22].Finding the 1150 o C isotherm is the same as finding the wear line.In the present non-linear IHTP, the unknown is the location of the 1150 o C isotherm.The observations are temperature measurements at points inside the blast furnace hearth.These measures are taken by thermocouples located inside the refractory walls of the hearth.
A number of papers deal with estimating the 1150 o C isotherm.For example, Gonzalez and Goldschmit in [9], solved the inverse geometry problem based on a radial basis functions (RBF) geometry representation.Zhao et al. [24] compared different characteristics of the allcarbon brick, such as the wear line, and the ceramic cup synthetic hearth bottoms.Brannbacka and Saxén [4] focused on the mathematical formulation of the problem, yielding a general model optimized for fast computation.Zagaria et al. [23] present a model they call MUSA that combines online readings of temperatures and historical data to predict hearth erosion and the skull profile.
This work solves the associated direct problem considering the Galerkin formulation for the finite element method.The thermal conductivity of the blocks making up the hearth's refractory walls are temperature-dependent, causing the direct problem to present a non-linearity.This non-linearity is handled with the fixed-point iteration method.The IHTP is solved using three methods: the Levenberg-Marquardt (LM), the iteratively regularized Gauss-Newton (IRGN), and the conjugate gradient (CG).The three methods are compared in performance and accuracy.The numerical experiments to find the 1150 o C isotherm were realized using a hearth similar to that of Blast Furnace 3 of steelmaker Arcelor Mittal Tubarão.The solution is validated through simulated measurements with different noise levels.
The paper is organized as follows.Section 2 presents details about the direct problem.Section 3 describes the main ideas about inverse heat transfer problems; it also delineates three methods to solve them -the Levenberg-Marquardt, the iteratively regularized Gauss-Newton, and the conjugate gradient methods.Section 4 presents the numerical experiments that consider the simulated measures of temperature in a hearth similar to that of Blast Furnace 3 (BF3).The paper closes with a summary of our main conclusions.

Direct problem
The direct problem associated with estimating a wear line consists of calculating the temperature distribution along the hearth.To solve the direct problem we use the Galerkin formulation for the finite element method [12].Only relatively small irregularities exist in the geometry, material properties, and boundary conditions in the angular direction of the hearth.Hence, the direct problem can be represented, as depicted in Fig. 1, as a two-dimensional axisymmetric problem.Moreover, since it has the characteristics of a thermal equilibrium problem, we can assume the problem to be a steady-state problem.Thus, the heat transfer problem of the hearth is governed by the Eq.(1): where k = k(T ) is the temperature-dependent thermal conductivity.
where T isotherm is the imposed temperature.This is the boundary condition to be estimated in the inverse problem.The boundaries 2 and 5 presents the Neumann boundary condition, where the heat flux is null due to the continuity of the domain on their respective normal directions: Boundaries 3 and 4 present the Robin boundary condition, which occurs due to the external cooling in these regions: where h is the convective heat transfer coefficient and T ∞ is the ambient temperature.

Galerkin finite element formulation
Considering the discretization of Ω in sub-domains Ω e , where e = 1, ..., n el , and n el is the number of elements.Γ e denotes the boundary of Ω e , assuming that ∪ e Ω e = Ω and ∩ e Ω e = ∅.
After doing the appropriate finite element approximation considering the linear triangular element, a linear system, shown in Eq. ( 5), is obtained from the Galerkin formulation, which is solved by the Diagonal Preconditioned Conjugate Gradient method [15]. where, ( where A is the assembling operator, u is the vector of nodal temperatures, N i is the finite element interpolation function for linear triangular elements, n e bc is the number of nodes that the Dirichlet boundary condition is considered and g j is the prescribed value to each node j.

Non-linearity of the thermal conductivity
The thermal conductivity of the refractory blocks from which the hearth is made is temperature-dependent.Thus, Eq.( 1) is actually: This non-linearity is handled with the fixed-point iteration method [15].Eq. ( 9) considers a non-linear iterative process where the thermal conductivity that will be used in the iteration iter + 1 is the one that was obtained in the iteration iter.
In the Blast Furnace 3, the thermal conductivity of the materials that compose the hearth has an approximately linear variation with respect to the temperature.Thus, in this work we consider: where T iter b is the temperature at the barycenter of the element e, obtained from the nodal temperatures of the element; k e is the thermal conductivity with linear behavior at the element e, and the parameters α and β are calculated from the thermal conductivity data.The initial condition T 0 considered is the direct problem approximation with a fixed value for the thermal conductivity for each material.The stopping criterion adopted for the fixed-point iteration method was the Euclidian norm T iter+1 − T iter with a tolerance of 10 −4 .

Inverse problem
The IHTP of estimating the hearth wear line consists of estimating the 1150 o C isotherm.This means it is necessary to estimate the value of the T isotherm that is the prescribed value shown in Eq. ( 2).To estimate this boundary condition, we compare three methods of solving inverse problems: the Levenberg-Marquardt (LM), the iteratively regularized Gauss-Newton (IRGN) and the conjugate gradient (CG).All three consist of minimizing the least squares norm: where S(P ) is the sum of squares to be minimized, called the objective function; Y is the vector of observed temperatures; T (P ) is the vector of estimated temperatures, and I is the number of measurements.The estimated temperatures T i (P ) are obtained from the solution of the direct problem at the domain's measurements positions (x meas , y meas ), using the vector P as the current value of the unknown parameters being estimated (T isotherm , in this problem).

Levenberg-Marquardt method
Levenberg [16] introduced his technique in 1944.In 1963, Marquardt [17] introduced a new approach to what was basically the same technique.Marquardt was trying to obtain a method that, in the neighborhood of the minimum of the ordinary least squares norm, would tend to the Gauss method.However, in the neighborhood of the initial guess used for the iterative procedure, Marquardt wanted the method to tend to the steepest descent method.The least squares norm used for this method is presented in Eq. (11).The iterative solution of the Levenberg-Marquardt method is given by: where P is the vector of parameters that are being estimated; J is the sensitivity matrix (see Section 4, for details); Y is the vector of observed temperatures; T (P ) is the vector of estimated temperatures; µ k is a positive scalar called the damping parameter, and Ω k is a diagonal matrix.Thus, the method consists of solving the following linear system obtained from Eq. ( 12): where P k+1 = P k + ∆P k .In our implementation, Ω k is calculated as follows: The damping parameter µ k is initially considered as µ 0 = 0.001.It is used to improve, in each iteration, the update of the parameters.If the S(P k+1 ) calculated in the current iteration is such that S(P k+1 ) ≥ S(P k ), the vector of parameters P (k+1) is not updated, and µ k is considered to be equal to 10µ k .If the S(P k+1 ) calculated in the current iteration is such that S(P k+1 ) < S(P k ), the vector of parameters P k+1 is updated to that calculated in the current iteration, and µ k+1 is considered be equal to 0.1µ k .As inverse problems are typically ill-posed, the term µ k Ω k in Eq. ( 12) helps with instabilities in cases where det(J T J) ≈ 0. The damping parameter (µ k ) is usually larger at the beginning of the iterations.This is because the problem, in the region around the initial guess used for the iterative procedure, is generally ill conditioned.With larger values for the damping parameter, the LM method tends to the steepest descent method.With smaller such values, the LM method tends to the Gauss-Newton method.

Iteratively regularized Gauss-Newton method
The iteratively regularized Gauss-Newton (IRGN) method was proposed by [3].The IRGN is basically the Gauss-Newton method with regularization techniques.Tests and proposals using the IRGN algorithm can be found in [2,6,9,13,20,21].The least squares norm used for this method is: where L T L is a regularization matrix; α is a regularization parameter; P 0 is the initial suggested value for the vector of estimating parameters P .Both being based on the Gauss-Newton method, the IRGN and LM methods resemble one another.The former, however, considers a different regularization method.The iterative solution of the IRGN method is given by: The solution P k+1 calculated with the IRGN method is used to update P k as follows: where β k is a step size that makes S(P k+1 ) < S(P k ).In our implementation we used β 0 = 1.0 and for every iteration where the condition S(P k+1 ) < S(P k ) was not accepted, we used β k+1 = 0.5β k .The regularization matrix L T L can be calculated as [6,9]: where: where w k ≥ 0 are weighting factors such that 2 k=0 w k = 1.The regularization parameter α is chosen such that α iter > 0 and: with r < 1.The value of α decreases with each iteration, and the first term can be determined as the optimal regularization parameter of the Tikhonov regularization method [7]: where δ is the noise level of the temperature measurements.

Conjugate gradient method
The conjugate gradient method is a powerful iterative technique for solving linear and non-linear inverse problems of parameter estimation.When it is used with appropriated stopping criterion, the method belongs to the class of iterative regularization techniques.Details about the convergence of the method can be found in [5,8,11,14].The least squares norm used for this method is done according to Eq. (11).The iterative solution of the CG method to minimize the function S(P ) can be defined by: where β k is the step size and d k is the direction of descent.The direction of descent can be calculated as: There are different expressions for calculating the coefficient γ k .In this paper we used the Polak-Ribiere [1, 5] definition, since it has improved convergence in nonlinear estimation problems [5,19]: where γ 0 = 0 for k = 0.The term [∇S(P k )] j can be calculated by differentiating Eq. ( 11) with respect to the vector of unknown parameters P : The search step size, β k , is calculated by minimizing S(P k+1 ) with respect to β k and can be given in the matrix form as:

Sensitivity matrix
The sensitivity matrix is given by: where N is the number of unknown parameters and I is the number of measurements.Each element of the sensitivity matrix is called a sensitivity coefficient.The sensitivity coefficient J ij is given by the first derivate of the temperature with respect to the unknown parameter P j : The derivate that appears in Eq.( 30), can be approximated by finite differences.Using a forward difference, the sensitivity coefficients can be approximated by: . ., P j + εP j , . . ., P N ) − T i (P 1 , P 2 , . . ., P j , . . ., P N ) where ε is used as ε = 10 −4 .

Stopping criterion
The following stopping criterion was used for all methods implemented in this work: where ε 1 and ε 2 are prescribed tolerances and • is the Euclidean norm.Equation (32) tests whether the least squares norm is small enough, which is expected for the solution of an inverse problem.The stopping criterion given by Eq.( 33) means that the iterative process should stop when the changes in the parameter vector are quite small.In this work we considered ε 1 = ε 2 = 10 −4 .

Errors
The errors in the estimated temperatures (Eq.( 34)) and in the estimated parameters (Eq.( 35)) are calculated, respectively, as follows: where ε T (P ) is the error in the estimated temperatures; ε P is the error in the estimated parameters; T (P ) is the vector of estimated temperatures; Y is the vector of observed temperatures; P is the vector of estimated parameters, and P real is the vector of real parameters.The values of the errors are given as percentages.

Numerical experiments
The numerical experiments of finding the 1150 o C isotherm are carried out by using simulated measures of temperature in a hearth (located at Arcelor Mittal Tubarão) similar to that of Blast Furnace 3 (BF3).Figure 2(a) shows the points A, B, C, D, E and F on the prescribed boundary, the hearth materials, the thermocouples positions, and the cooling characteristics of the convective boundaries.Figure 2(b) shows the mesh used, generated by Easymesh2 , with 5,881 nodes and 11,382 elements.Table 1 shows the convective parameters of the cooling boundaries and Table 2 shows the thermal conductivity of the hearth materials.
Different values were considered for the weighting factors in the regularization method of the IRGN.All of them, however, led to similar behaviors, such as the number of iterations, computational times and convergence.Since they were so similar, we showed only the results of one combination of these parameters (w 0 = 0.4, w 1 = 0.3, w 2 = 0.3).
Two experiments were done: (i) Experiment 1 estimates the T AF (see, Fig. 2    To generate the simulated measurements, the direct problem is solved considering that all parameters involved with the problem are known, including the T isotherm , the parameter to be estimated by the inverse method.After solving the direct problem with the parameters being considered real, the simulated measurements are generated as follows: where: . ., Y n measuments ] is the vector with real measurements and is the input data of the inverse problem; . ., T n measuments ] is the solution vector obtained from the estimated parameters.δ is the random noise level. We have set value ranges to generate the noise levels.For example, 5% noise is generated by using values between 0 and 5%; 10% noise is generated by using values between 5 and 10%; 15% noise is generated by using values between 10 and 15%, and 20% noise is generated by using values between 15 and 20%.We chose 32 positions for the simulated measurements.Fig. 2(a) shows the thermocouple configuration.
In both experiments, we solved the inverse problem five times for each range of noise level and calculated the average of these five solutions.For each time, we consider the same set of random noise levels for all algorithms and for both experiments.Thus, we can compare the convergence of the estimated parameters for all algorithms and compare the quality of the solution between both experiments.
In both experiments, we considered the temperature of 1400 o C for the prescribed boundary (Eq.( 2 Table 3 compares the three algorithms for Experiment 1.The three algorithms were very similar for this first experiment.They all converged practically to the same value.The number of iterations and the computational time were similar, but in terms of computational time the CG method achieved, at almost all noise levels, a better result.Figure 3 shows the estimated wear for Experiment 1.The wear for each noise level is related to the results in Tab. 3. As we can see, the estimated wear line for different noise levels is almost the same.Table 4 compares the three algorithms for Experiment 2. In this experiment, we can see that the CG method requires more computational time and more iterations.Once again, all three methods converged practically to the same value.We can see that, the LM method achieved results similar to those of the IRGN method.However, the IRGN's results were better with fewer iterations and less computational time.
Figure 4 shows the estimated wear for Experiment 2. The wear for each noise level is related to the results in the Table 4. Once again, the estimated wear line was similar to all the noise levels used, but in Experiment 1 the errors in the estimated parameters were smaller.
The results of Experiment 2, were more sensitive than those of Experiment 1 to the noise in the simulated measurements.Since the second experiment had a larger degree of freedom, that behavior was expected.As both experiments used the same value of noise levels for each related test, we can tell that the second test led to a solution with smaller error in the estimated   temperatures (ε T (P ) ).However, as evidenced in Tabs. 3 and 4, smaller errors in the estimated temperatures do not always mean smaller errors in the estimated parameters.Figures 5,6,7 show, respectively, the behaviors of LM, CG, and IRGN errors for Experiment 2. The value of the vector ε T (P ) decreases, for all algorithms, at each iteration.This is expected, since all of them guarantee that the objective function decreases in each iteration.The reduction in the value of ε T (P ) does not necessarily mean a reduction of the error in the estimated parameters.This behavior appears in accordance with the increase in noise level.Since the domain is a radial section of the blast furnace hearth, we can generate a tridimensional view of the wear line.This is done by rotating the bi-dimensional domain around the center of the hearth, even though it is known that the wear varies across radial sections.

Conclusions
We presented a model that estimates, using inverse problems, the wear line in a blast furnace hearth.On the temperature map along the hearth, the wear line is considered the isotherm of 1150 o C, the lowest temperature at which iron can exist in liquid form.The inverse problem was solved with three methods: the Levenberg-Marquardt, the conjugate gradient, and the iteratively regularized Gauss Newton.
Two experiments were done.In the first one, we estimated only one parameter, the temperature along the prescribed In the second one, we divided the prescribed boundary into five parts, and estimated those parameters.The experiment to estimating one parameter was less sensitive to variations in noise level than was the experiment estimating five parameters.This was due to the latter having more degrees of freedom.Nevertheless, since the prescribed boundary is an approximation of the real problem, it was interesting to divide it to enable future tests with real temperature measures.
For the experiment estimating a single parameter, the three methods presented quite similar behavior -similar number of iterations and similar computational time.For the experiment estimating five parameters, the CG method needed more iterations, and thus more computational time, to converge.As expected, since they have similar formulations, the IRGN and the LM methods were similar in their number of iterations and computational time.

Figure 1 .
Figure 1.Simplified domain (a)), which is considered the prescribed value through Segment AF ; (ii) Experiment 2 estimates the T AB , T BC , T CD , T DE and T EF which are the five parts of the prescribed boundary (segments AB, BC, CD, DE and EF ).

Figure 2 .
Figure 2. Materials and mesh of the BF3.

Figure 3 .
Figure 3.Estimated wear line for different noise levels for the Experiment 1.

Figure 4 .
Figure 4.Estimated wear line for different noise levels for the Experiment 2