Explicit Sensitivity Coe ﬃ cients for Estimation of Temperature-Dependent Thermophysical Properties in Inverse Transient Heat Conduction Problems

: Explicit expressions are obtained for sensitivity coe ﬃ cients to separately estimate temperature-dependent thermophysical properties, such as speciﬁc heat and thermal conductivity, in two-dimensional inverse transient heat conduction problems for bodies with irregular shape from temperature measurement readings of a single sensor inside the body. The proposed sensitivity analysis scheme allows for the computation of all sensitivity coe ﬃ cients in only one direct problem solution at each iteration with no need to solve the sensitivity and adjoint problems. In this method, a boundary-ﬁtted grid generation (elliptic) method is used to mesh the irregular shape of the heat conducting body. Explicit expressions are obtained to calculate the sensitivity coe ﬃ cients e ﬃ ciently and the conjugate gradient method as an iterative gradient-based optimization method is used to minimize the objective function and reach the solution. A test case with di ﬀ erent initial guesses and sensor locations is presented to investigate the proposed inverse analysis.


Introduction
The accuracy of the numerical simulation of heat transfer problems relies significantly on the accuracy of data, including, among others, the thermophysical properties such as the thermal conductivity and the specific heat of heat-conducting body. If, no a priori information on the thermophysical properties is available, inverse methods, as inexpensive alternatives to expensive experiments with sophisticated instruments, may be employed to estimate the unknown properties accurately. The inverse methods are widely used to estimate the thermal conductivity and the specific heat. Sawaf et al. [1] estimated linearly temperature-dependent thermal conductivity and specific heat capacity for an orthotropic solid using an inverse analysis. The Levenberg-Marquardt iterative procedure is used to minimize the objective function. Flach and Özişik [2] employed an inverse analysis to estimate spatially varying thermal conductivity and heat capacity per unit volume of an one-dimensional slab using the Levenberg-Marquardt method. Talukdar et al. [3] used ant colony optimization algorithm for the estimation of temperature-dependent thermal conductivity and specific heat. Mohebbi et al. [4] derived explicit sensitivity coefficients to estimate the linearly temperature-dependent thermal conductivity in steady-state heat conduction problems using an inverse analysis. The conjugate-gradient method is used as a minimization method to minimize the objective function and estimate the unknown parameters. Liu [5] used a hybrid method based on a combination of the modified genetic algorithm and Levenberg-Marquardt method to identify simultaneously the fluid thermal conductivity and heat capacity for a transient inverse heat transfer - The steady-state or transient heat conduction problems are concerned with regular bodies only (inability to consider the irregular bodies) and the heat conduction equation is solved using the traditional finite-difference method. -Most of the boundary conditions considered are associated with either a constant temperature (Dirichlet boundary condition) or an insulated surface. Hence, the Neumann and Robin boundary condition types are not addressed. -Most of the earlier works have been limited to the one-dimensional heat conduction problems.
Therefore, a sufficiently accurate general methodology for solving transient heat conduction problems in the presence of temperature-dependent thermophysical properties, which can handle a general 2D domain and a variety of boundary conditions, is required. This study is concerned with direct and inverse transient heat conduction problems in general 2D domains with different boundary conditions.
The inverse heat transfer problems are mathematically challenging problems due to their ill-posed nature. They are inherently unstable and very sensitive to noise and special methods are required to treat them. Among such methods are iterative regularization methods [7,8]. In these methods, the original objective function expression, defined as a nonlinear least-square formulation, is not modified and the minimization of the objective function is performed by a gradient-based minimization method such as conjugate-gradient method or steepest descent method and the discrepancy principle can be used as a criterion to stop the iteration process and obtain a reasonably stable solution. As this study is concerned with a parameter estimation problem, the gradient of the objective function with respect to the unknown parameters (as needed in gradient-based minimization methods) can be computed by using the finite-difference method (by additional direct problem solutions and forming the sensitivity matrix) or adjoint method. Both methods involve an increase in computational cost. Moreover, the adjoint method has its own mathematical complexity. In this study, however, using the obtained explicit sensitivity coefficients, the gradient of the objective function can be computed in only one direct problem solution at each iteration, thereby decreasing the computational cost significantly.

Governing Equation
The heat conducting body shown in Figure 1a is made of a material whose thermal conductivity and specific heat are linearly temperature-dependent variables, namely k T = k 1 + k 2 T and c = c 1 + c 2 T where k 1 , k 2 , c 1 , and c 2 are constant. The density of the body is ρ and the body is initially at the temperature T 0 . For the time t > 0, a time wise varying heat flux . q(t) is applied at the boundary surface Γ 1 . Convective heat transfer is imposed on boundary surfaces Γ i , i = 2, 3, 4 with corresponding heat transfer coefficients h i , i = 2, 3, 4 and surrounding temperatures T ∞ i , i = 2, 3, 4. For this problem, the two-dimensional transient heat conduction equation with no heat generation is for the linearly temperature-dependent thermal conductivity and the specific heat, namely k T = k 1 + k 2 T and c = c 1 + c 2 T, Equation (1) is expressed as ∂T(x,y,t) ∂x ∂T(x,y,t) ∂y ∂T(x,y,t) ∂t (2) ∂T(x,y,t) ∂t For this problem, the two-dimensional transient heat conduction equation with no heat generation is for the linearly temperature-dependent thermal conductivity and the specific heat, namely = + by simplifying Equation (2), we get The boundary and initial conditions are ¶ = ¶ where t is the time. Since the x and y physical domain (the geometry of the heat-conducting body) is irregular, it is mapped onto a regular one (the ξ and η computational domain). The elliptic The boundary and initial conditions are ∂T(x, y, t) T(x, y, 0) = T 0 (x, y) in physical domain Ω(x, y) where t is the time. Since the x and y physical domain (the geometry of the heat-conducting body) is irregular, it is mapped onto a regular one (the ξ and η computational domain). The elliptic grid generation method (by solving two Laplace equations ξ xx + ξ yy = 0 and η xx + η yy = 0) is used to mesh the physical domain. Then the transient heat conduction equation and its associated boundary and initial conditions can be transformed from the (x, y, t) to the (ξ, η, t) variables [7, [9][10][11]. The transformation gives where are the coefficients obtained from the elliptic grid generation method. By simplifying Equation (7), we get The transformed boundary and initial conditions become where T * 0 (ξ, η) is the initial condition T 0 (x, y) rewritten in terms of the variables ξ and η. The derivatives appearing in the transformed heat conduction equation, Equation (9), can be discretized in the regular computational domain using the finite-difference method, as follows (assuming ∆ξ = ∆η = 1) where f ≡ x, y, T. One-sided forward and one-sided backward relations should be used to discretize the boundary condition equations. Using forward-time-central-space (FTCS) discretization and the relations in Equation (15), and employing an explicit approach, we can approximate the differential Equation (9) as follows where ∆t and n are the time step and the time level, respectively. By considering the stability criterion for the solution of Equation (16), a time-marching procedure can be used to obtain T n+1 i,j . In other words, the nodal temperatures at the time level n + 1, T n+1 i,j , can be determined from the knowledge of nodal temperatures at the previous time level n, T n i,j , as follows

Objective Function
The objective of this study is to separately estimate the linearly temperature-dependent thermophysical properties such as the thermal conductivity and the specific heat using the transient readings of a single sensor S placed at the point (Si, S j) inside the body. To do so, using an inverse analysis the square of the difference between the estimated temperatures at the sensor place,T e (Si, S j, t i ), computed from the solution of the direct heat conduction problem using the estimated thermophysical properties and the measured temperatures T m (Si, S j, t i ) over the time domain 0 < t < t f is minimized. This can be mathematically expressed as where C is a positive constant and can be considered as C = 10 n , n = 1, 2, 3, . . . The inverse analysis is used to minimize the following objective function expression:

Sensitivity Analysis
The calculation of the gradient of the objective function J defined by Equation (19) with respect to the unknown variables, k 1 , k 2 , c 1 , and c 2 is needed in the inverse problem. Thus, we can write where P ≡ k 1 , k 2 , c 1 , c 2 . In Equation (20), C ∂T e (Si,Sj,t i ) ∂P (i = 1, f ) are called the sensitivity coefficients and can be obtained by differentiating the obtained expression for T n+1 e (Si, S j, t i ) in Equation (17) with respect to k 1 , k 2 , c 1 , and c 2 , as follows where α, β, γ, J, T ξ , T η , T ξξ , T ξη , and T ηη are computed using the finite-difference expressions associated with the sensor place, S(Si, S j) at the time t i , i = 1, f . By considering Equations (21)-(24), it can be seen that the inverse problem of estimating k 1 and k 2 is linear (as the sensitivity coefficient expressions are independent of k 1 and k 2 ) and the inverse problem of estimating c 1 and c 2 is nonlinear (as the sensitivity coefficient expressions depend on c 1 and c 2 ).
Moreover, it can be noticed that all sensitivity coefficients C can be computed in only one single direct problem solution without the need for solving sensitivity and adjoint problems. The sensitivity matrix Ja can be explicitly written as and

The Conjugate Gradient Method (CGM)
In this study, the conjugate gradient optimization method is used to solve the inverse heat transfer problem. The objective function given by Equation (19) is minimized by searching along the direction of descent d (k) using a search step length β (k) .
where P ≡ k 1 , k 2 , c 1 , c 2 are parameters to be estimated. The direction of descent of the current iteration is obtained as a linear combination of the direction of descent of the previous iteration and the gradient direction ∇J (k) . Therefore, The Polak-Ribiere formula [12] is used to compute the conjugation coefficient: The search step-length is given as follows [7] The parameter update in the conjugate-gradient or steepest-descent methods can be written as < 0 which may result in a negative thermophysical property which is unacceptable. So, the objective function term (least squares norm) is multiplied by the constant C to prevent a negative thermophysical property. The multiplication of the objective function by a constant does not change the optimal solution. Then, from Equation (20), the gradient of the objective function, ∇J (k) , is also multiplied by C. As the sensitivity coefficients all are multiplied by the constant C, from Equation (31) we find that the search-step length β (k) is multiplied by 1 . This way the updated parameter P (k+1) will be positive and the iteration process can continue. we can interpret this matter as follows: The incorporation of the constant C into the objective function definition has no effect on both the optimal solution and the gradient and only decreases the step-length as 1 C of its value when the constant C is not considered (i.e., we initiate the iterative process with a smaller step-length): The value of C = 10 n , n = 1, 2, 3, . . . is a convention based on numerical experience which works well for the parameter estimation problems. In this study, until C = 10 3 , we still get negative thermophysical property and the iteration process is terminated and hence C = 10 4 is chosen to initiate the iterations.

Optimization Algorithm
The following algorithm presents the direct and inverse analysis steps used to estimate the linearly temperature-dependent thermal conductivity k T = k 1 + k 2 T and the specific heat c = c 1 + c 2 T separately. To do so, either two parameters k 1 and k 2 for the thermal conductivity or two parameters c 1 and c 2 for the specific heat are estimated simultaneously. In other words, it is assumed that k T = k 1 + k 2 T is known when estimating c = c 1 + c 2 T, and vice versa. Simultaneous estimation of all parameters k 1 , k 2 , c 1 , and c 2 is not considered in this study. - Specify the physical domain, the boundary and initial conditions, the thermophysical properties, and the measured temperatures at the sensor place S Si,Sj and the time t Generate the boundary-fitted grid over the heat-conducting body using the elliptic grid generation method. 3. Solve the direct problem to obtain the temperature values at the sensor place and the time Using Equation (19), compute the objective function (J (k) ).

5.
If value of the objective function obtained in step 4 is less than the specified stopping criterion, the optimization is finished. Otherwise, go to step 6. 6.
Compute the sensitivity matrix Ja k 1 from Equation (26).

8.
Compute the conjugation coefficients γ 10. Compute the search step-length β 12. With new value for k 1 repeat the steps 6 to 11 for k 2 . 13. Set the next iteration (k = k + 1) and return to the step 2.
-Estimation of c = c 1 + c 2 T (assuming k T = k 1 + k 2 T is known): Repeat the steps 1 to 13 given above for the simultaneous estimation of c 1 and c 2 .

Stopping Criterion
If there are no errors in the temperature measurements, then the iterations can be stopped when where ε is a small specified number chosen based on obtaining stable and appropriate results. In this study, for the case of no measurement error, ε = 0.5. However, if there are errors in the temperature measurements, the discrepancy principle is used to stop the iterative process. This principle states that if the difference between estimated and measured temperatures is of the order of magnitude of the measurement errors, that is, then the solution is assumed to be sufficiently accurate and the iterative process can be terminated. In Equation (33), σ is the standard deviation of the measurement errors and is assumed constant in this study (σ = 0.5). If we substitute Equation (33) into Equation (19) (objective function definition), then the value for ε can be obtained as Then the iterative process can be stopped when the following criterion is satisfied In this study, the measured temperatures containing random errors, T meas (Si, S j, t i ), (i = 1, f ), are generated by adding an error term ωσ to the exact temperatures T exact (Si, S j, t i ) to give where ω is a random variable with normal distribution, zero mean, and unitary standard deviation. Assuming 99% confidence for the measured temperature, ω lies in the range −2.576 ≤ ω ≤ 2.576 and it is randomly generated by using MATLAB.

Simultaneous Estimation of the Parameters
The simultaneous estimation of all four parameters k 1 , k 2 , c 1 , and c 2 traps in local minima (close to the global minima). It is common to combine a global optimization method like genetic algorithm with gradient based ones (such as what is considered in this study) to have both good initial guesses and convergence speed to reach a global solution (as in [5]). In this work, the accurate estimation of any pair of (k 1 , k 2 ), (c 1 , c 2 ) and constant parameters (c 1 , k 1 ) is possible (global solution). Moreover, it should be noted that only one single sensor (inside the body) is used in this work which makes the minimization problem difficult. The investigation of the effect of using multiple sensors (inside or at the boundary surface), the place of sensors, and linear-dependence or independence of all sensitivity coefficients are not included in this work as the objective of this work is to derive explicit expressions for the sensitivity coefficients and verify their accuracy. Moreover, the estimation of the parameters separately does not limit the practical importance of the sensitivity coefficients obtained in this study. The reason is that at first the temperature-dependent thermal conductivity components k 1 and k 2 can be estimated using the method the author proposed in [4]. In this method, using a steady-state heat transfer analysis (as k 1 and k 2 are constant in both steady-state and transient analyses) and temperature measurements at a boundary surface (not inside the body), these two parameters can be estimated. Then, by having k 1 and k 2 , two parameters c 1 and c 2 can be estimated using the method presented here (using a transient heat transfer analysis). It is worth mentioning that the sensitivity coefficients expressions developed here for k 1 and k 2 are applied for the points located inside a heat-conducting body and are different from the ones developed before for k 1 and k 2 in [4] which are applied for the points located at a boundary surface.

Results
A test case is presented here to investigate the accuracy and efficiency of the proposed sensitivity analysis scheme to estimate the linearly temperature-dependent thermal conductivity k T = k 1 + k 2 T and the specific heat c = c 1 + c 2 T separately. Initially it is assumed that the thermophysical property to be estimated is known, the transient heat conduction problem is then solved to estimate the temperature at the sensor place at times t i (i = 1, f ). Then, the estimated temperatures are employed as simulated measured ones to recover the initially used thermophysical property. The numerical values of the coefficients involved in the test case are listed in Table 1.
The grid size used in this study is M × N = 20 × 20 (Figure 2), the initial temperature is Before proceeding further, the accuracy and correct implementation of direct heat conduction solver should be confirmed. In order to validate the results from the direct problem solution, the temperature distribution for the heat conducting body used in this study is compared to the temperature distribution obtained from the commercial finite element software COMSOL. To do so, the timewise varying heat flux 50000 + 40000 sin( πt i 180 ) is considered. The temperature distribution obtained by the explicit solver, Equation (17), using the problem data given in Table 1, a grid size of 20 × 20, and the time step ∆t = 0.1 s is shown in Figure 3a and the temperature distribution obtained by the finite element software COMSOL is depicted in Figure 3b. The temperature histories at the sensor places (i, j) = (3, 12) ( Figure 4a) and (i, j) = (18, 8) (Figure 5a) obtained by both the explicit solver and the software COMSOL are compared and presented in Figure 4b or Figure 5b, respectively. Moreover, the temperature distribution at the final time t f = 500 s at a given line, here along the nodes (i, j), i = 1, . . . , M, j = 10 ( Figure 6a) obtained by the explicit solver and the software COMSOL is shown in Figure 6b. The comparison between the results shows a very good agreement, thereby confirming the correct implementation of the explicit solver.
), m . C 2,3,4 The grid size used in this study is Before proceeding further, the accuracy and correct implementation of direct heat conduction solver should be confirmed. In order to validate the results from the direct problem solution, the temperature distribution for the heat conducting body used in this study is compared to the temperature distribution obtained from the commercial finite element software COMSOL. To do so, the timewise varying heat flux + 50000 40000 sin( ) 180 i πt is considered. The temperature distribution obtained by the explicit solver, Equation (17), using the problem data given in Table 1, a grid size of 20 20 , and the time step is shown in Figure 3a and the temperature distribution obtained by the finite element software COMSOL is depicted in Figure 3b. The temperature histories at the sensor places = ( , ) (3,12) i j ( Figure 4a) and = ( , ) (18, 8) i j (          Two different places are considered for the single sensor S to study the effect of the sensor location on the final results. The sensor is placed at the grid nodes (i, j) = (3,12) and (i, j) = (18, 8) successively. Once the nodal temperature at the sensor place is calculated at each time t i , the sensitivity coefficients can be immediately computed during the transient solution using the explicit expressions obtained previously. This implies that the proposed sensitivity analysis scheme is very efficient and does not contribute significantly to the computational cost of the solution. As stated previously, we assume that c = c 1 + c 2 T is known when estimating k T = k 1 + k 2 T, and vice versa. Initially we assume that the thermal conductivity is known and is equal to k T = 12.5 + 0.05T. The aim of inverse problem is now to recover the initially used specific heat c = 450 + 0.02T using two different initial guesses: These different initial guesses are selected so that they can reflect the accuracy and efficiency of the inverse analysis. The results of the recovery of the specific heat components c 1 and c 2 as well as the history of the objective function during the minimization process are shown in Figure 7 (for the initial guess c initial 1 ) and Figure 8 (for the initial guess c initial 2 ) by using the sensor placed at the node (i, j) = (3, 12) and Figure 9 (for the initial guess c initial 1 ) and Figure 10 (for the initial guess c initial 2 ) by using the sensor placed at the node (i, j) = (18, 8) and Figure 11 (for the initial guess c initial 2 ) by using the sensor placed at the node (i, j) = (18, 8) and the measurement error of σ = 0.5. As shown in Table 2, without the measurement error, a 100% reduction in the objective function value and complete recovery of the unknown specific heat components are achieved by initiating the minimization process from both initial guesses, which shows the accuracy of the proposed sensitivity analysis scheme. Moreover, it can be seen that the proposed inverse analysis is not strongly affected by the errors involved in the temperature measurements, and the unknown parameters are recovered accurately.
Then we assume that the specific heat is known and is equal to c = 450 + 0.02 T. The aim of inverse problem is then to recover the initially used thermal conductivity k T = 12.5 + 0.05 T using two different initial guesses: The results of the recovery of the thermal conductivity components k 1 and k 2 as well as the history of the objective function during the minimization process are shown in Figure 12 (for the initial guess k T initial 1 ) and Figure 13 (for the initial guess k T initial 2 ) by using the sensor placed at the node (i, j) = (3, 12) and Figure 14 (for the initial guess k T initial 1 ) and Figure 15 (for the initial guess k T initial 2 ) by using the sensor placed at the node (i, j) = (18, 8) and Figure 16 (for the initial guess k T initial 1 ) by using the sensor placed at the node (i, j) = (3, 12) and the measurement error of σ = 0.5. The summary of the results given in Table 3 reveals that without the measurement error a 100% reduction in the objective function value is obtained and the unknown thermal conductivity components are recovered accurately by initiating the minimization process from both initial guesses, which again shows the accuracy of the proposed sensitivity analysis scheme. Like the specific heat estimation case, it can be seen that the estimation of thermal conductivity components is not strongly affected by the errors involved in the temperature measurements and the unknown parameters are recovered accurately. by using the sensor placed at the node = ( , ) (18, 8) i j and Figure 16 (for the initial guess initial 1 T k ) by using the sensor placed at the node = ( , ) (3,12) i j and the measurement error of = 0.5 σ . The summary of the results given in Table 3 reveals that without the measurement error a 100% reduction in the objective function value is obtained and the unknown thermal conductivity components are recovered accurately by initiating the minimization process from both initial guesses, which again shows the accuracy of the proposed sensitivity analysis scheme. Like the specific heat estimation case, it can be seen that the estimation of thermal conductivity components is not strongly affected by the errors involved in the temperature measurements and the unknown parameters are recovered accurately.     Figure 16. Estimation of k 1 (a) and k 2 (b), and the objective function versus iteration number (c) using the initial guess k T initial 1 = 30.0 + 0.5 T( W m. • C ) and a single sensor placed at the node (i, j) = (3, 12) by considering the measurement error of σ = 0.5). Table 2. A summary of results for the estimation of the temperature-dependent specific heat components c 1 and c 2 (c = c 1 + c 2 T).  Table 3. A summary of results for the estimation of the temperature-dependent thermal conductivity components k 1 and k 2 (k T = k 1 + k 2 T).

Conclusions
Novel explicit expressions were derived for the sensitivity coefficients to estimate the temperature-dependent thermal conductivity and the temperature-dependent specific heat in general two-dimensional heat-conducting bodies using an inverse transient heat conduction analysis. Due to inability of the traditional finite-difference method to effectively handle the solution of heat transfer problems involving the irregular shapes, the irregular heat-conducting body was transformed into a regular computational domain to perform all computations related to the direct and inverse transient heat conduction solution using the finite-difference method, a method chosen for its widespread use in numerical heat transfer and ease of implementation. The irregular body was meshed using the elliptic grid generation technique because of its capability of producing a smooth grid over the body. Then, the governing equation and the associated boundary conditions were solved in the computational domain to obtain the temperature distribution over the body at each time step. In the inverse analysis, an objective function was defined based on the least squares norm and the explicit expressions were obtained for the sensitivity coefficients by differentiating the obtained expression for the temperature at the sensor location with respect to the unknown parameters. In addition to being accurate expressions, the use of explicit sensitivity coefficients eliminates the need for extra solution of direct problem to calculate the sensitivity coefficients using the finite-difference method or eliminates the need for solving adjoint problem to compute the gradient of the objective function with respect to the unknown parameters. We showed that, using the derived explicit sensitivity coefficients, the sensitivity matrix was computed in only one direct problem solution (at each iteration) and hence the computational cost was decreased significantly. The conjugate gradient method, as an iterative gradient-based optimization method, was used to minimize the objective function and reach the solution. A test case with different initial guesses and the sensor locations was presented to investigate the proposed inverse analysis. The obtained accurate results showed that the estimation of the temperature-dependent thermal conductivity or specific heat using the proposed method is not affected by the initial guess or sensor placement. Likewise, the results revealed that the method is not affected significantly by the error in the temperature measurement. This confirms that the proposed inverse analysis is very accurate, robust, and efficient.