Function Estimation in Inverse Heat Transfer Problems Based on Parameter Estimation Approach

: A new sensitivity analysis scheme is presented based on explicit expressions for sensitivity coefficients to estimate timewise varying heat flux in heat conduction problems over irregular geometries using the transient readings of a single sensor. There is no prior information available on the functional form of the unknown heat flux; hence, the inverse problem is regarded as a function estimation problem and sensitivity and adjoint problems are involved in the solution of the inverse problem to recover the unknown heat flux. However, using the proposed sensitivity analysis scheme, one can compute all sensitivity coefficients explicitly in only one direct problem solution at each iteration without the need for solving the sensitivity and adjoint problems. In other words, the functional form of the unknown heat flux can be numerically estimated by using the parameter estimation approach. In this method, the irregular shape of heat-conducting body is meshed using the boundary-fitted grid generation (elliptic) method. Explicit expressions are given to compute the sensitivity coefficients efficiently and the steepest-descent method is used as the minimization method to minimize the objective function and reach the solution. Three test cases are presented to confirm the accuracy and efficiency of the proposed inverse analysis.


Introduction
Direct heat transfer problems are concerned with the determination of temperature distribution in a heat-conducting body from known values for thermo-physical properties, geometrical configuration, boundary conditions, and heat flux.In contract, inverse heat transfer problems (IHTPs) deal with the determination of the thermo-physical properties, the geometrical configuration, the boundary conditions, and the heat flux from the knowledge of the temperature distribution within or on some part of the boundary of the heat-conducting body.Unlike the direct heat transfer problems which are well-posed, the inverse heat transfer problems are ill-posed and inherently unstable and extremely sensitive to measurement errors.Due to the ill-posed nature of IHTPs, these problems are widely regarded as mathematically challenging problems and many methods have been proposed to address the underlying challenges.Over the past decades, the numerical methods for solving the inverse heat transfer problems have received much attention due to the ever-increasing power of computers.Among the numerical methods used to overcome the instabilities in inverse heat transfer problems are the iterative regularization methods.In these gradient-based methods, the inverse problem solution is improved sequentially, and the number of iterations required to obtain a stable solution is set based on the measurement errors using the discrepancy principle [1][2][3].
During the past decades, inverse heat transfer analysis has been extensively used to estimate the constant and variable heat flux using temperature measurement taken at some points inside the heat conducting body or on some part of the boundary of the body.In [4], a three-dimensional inverse heat conduction problem is solved to determine the unknown transient boundary heat flux in an irregular domain using simulated temperature measurements taken at some appropriate locations and times on different parts of the boundary.Then, the least-squares functional is minimized using the conjugate gradient method.In [5], the inverse problem deals with the estimation of the spacedependent heat flux in two-dimensional steady-state heat conduction problems in the presence of variable thermal conductivity.The conjugate gradient method is used as a minimization method to minimize the least-squares objective function.In [6], the inverse problem is concerned with the simultaneous estimation of heat transfer coefficient and heat flux imposed on different parts of the boundary of an eccentric hollow cylinder.The thermal conductivity of the cylinder varies with space and is not considered constant.In [7], an inverse analysis is used to estimate an imposed heat flux on the inner surface of a long cylinder in a transient heat conduction problem using the measured temperature on the outer surface of the cylinder.In [8], a time-dependent heat flux applied at the inner surface of a functionally graded hollow cylinder is estimated via inverse analysis.The simulated temperature measurements are taken within the cylinder and the minimization of the functional is performed by the conjugate gradient method.In [9], an inverse analysis is employed to estimate the unknown base heat flux in an irregular fin made of a material with space-dependent thermal conductivity.The temperature measurements are taken within the fin and the conjugate gradient is used to minimize the objective function.In [10], using an inverse analysis, the heat flux applied on the high-temperature wall of an engine is estimated.The thermal conductivity is a function of temperature and the imposed heat flux is a function of time and position.In [11], the surface heat flux is estimated from two temperature measurements taken inside a finite domain with temperature-dependent thermal properties using a sequential inverse method.In [12], a timedependent wall heat flux applied on the surface of a solid medium is estimated using an inverse analysis in which the minimization of the functional is performed by the Levenberg-Marquardt method.The temperature measurement data are obtained from thermocouples embedded inside the solid medium.
Inverse problems can be regarded either as a parameter estimation or as a function estimation approach [2].The parameter estimation approach can be used whenever the unknown quantity can be expressed by a few parameters.For example, the recovery of a constant thermal conductivity can be viewed as a parameter estimation approach.However, if no prior information is available on the functional form of the unknown quantity, then the function estimation approach should be used to estimate the unknown functional form of the unknown quantity.To the best of the author's knowledge, this study is the first one in which the unknown functional form of a variable heat flux is estimated by the parameter estimation approach through the simultaneous estimation of so many parameters efficiently and accurately.
The function estimation approach to estimate timewise varying surface heat flux, when no a priori information is available on the functional form of the variable heat flux, involves the solution of sensitivity and adjoint problems which require additional mathematical developments and have computational costs equivalent to that of the direct problem.In this study, the inverse problem is formulated based on the parameter estimation approach to solve the function estimation problem of the recovering of the timewise varying heat flux with no a priori information of the functional form of the unknown variable heat flux.
To do so, a two-dimensional transient heat conduction equation subjected to appropriate initial and boundary conditions (with Neumann and Robin conditions at boundaries of a general irregular heat conducting body) is solved using an explicit scheme.Initially, the two-dimensional irregular domain is transformed into a regular computational domain, and all computations needed to solve the direct and inverse problem equations are performed in the regular computational domain.A body-fitted grid generation method (here the elliptic grid generation method) is used to mesh the physical domain, and the finite-difference method, a method chosen for its easy-to-implement feature, is used to approximate the derivatives of the field variable (temperature) at the grid nodes by algebraic expressions.A novel, efficient, accurate, and very easy to implement sensitivity analysis scheme is developed to compute the sensitivity (Jacobian) matrix by employing the chain rule to relate the temperature at the sensor place and the timewise varying wall heat flux applied on the part of the domain boundary.As mentioned above, the inverse analysis using the proposed sensitivity analysis scheme is not involved with the sensitivity and the adjoint equations, thereby reducing the associated development efforts and the computational costs significantly.The main novelty of the proposed sensitivity analysis is that all sensitivity coefficients can be computed efficiently in only one direct problem solution (during the transient solution), with no need for the solution of the sensitivity and adjoint equations (to compute the gradient of the objective function with respect to the variables), irrespective of the number of unknown parameters, which is extremely large in this study.A nonlinear least-square formulation is used to define the objective function and the steepest-descent method with an appropriate stopping criterion specified by the discrepancy principle is employed to minimize the objective function and recover the unknown functional form accurately. Three test cases are presented to reveal the accuracy and efficiency of the employed inverse analysis.In this study, three different measurement errors will be considered.As will be shown, the proposed inverse analysis is not strongly affected by the errors involved in the temperature measurements and the unknown timewise varying heat flux can be recovered with good accuracy.
The inverse analysis for the transient heat conduction problem presented in this study is sufficiently general; that is, it can be used to estimate the timewise varying wall heat flux applied on part of the boundary of a general irregular two-dimensional region (a heat-conducting body) with both Neumann and Robin conditions at the boundaries as long as the general two-dimensional region can be mapped onto a regular computational domain.

Governing Equation
The heat-conducting body shown in Figure 1a  For this problem, the two-dimensional transient heat conduction equation with no heat generation is (1) The boundary and initial conditions are ( , , )-on boundary surface Γ ( , ), 2, 3, 4 where t is the time.Since the heat conduction problem is concerned with an irregular body, the irregular shape (the x and y physical domain) is mapped onto a regular one (the ξ and η computational domain).The elliptic grid generation method is used to mesh the domain.Then, the heat conduction equation and its associated boundary and initial conditions can be transformed from the ( , , x y t ) to the ( , , ξ η t ) variables [2,13].The transformation gives smooth grid over the physical domain is generated.Thus, Equation (5) becomes , for 0 where are the coefficients obtained from the elliptic grid generation method.The transformed boundary and initial conditions becomes where * 0 ( , ) T ξ η is the initial condition 0 ( , ) T x y rewritten in terms of the variables ξ and η .The derivatives appearing in the transformed heat conduction equation, Equation( 6), can be discretized in the regular computational domain using the finite-difference method, as follows (assuming where  , , f x y T .For the discretization of the boundary condition equations, one-sided forward and one-sided backward relations should be used.The transformed differential equation, Equation (6), can be approximated by the finite-difference method using the explicit method.Using forward-time-central-space (FTCS) discretization and the relations in Equation ( 13), we get where t is the time step.Considering stability criterion, Equation ( 14) can be solved using the time-marching procedure to obtain . In other words, the nodal temperatures at the time level , can be determined from the knowledge of nodal temperatures at the previous time level n , ,

Objective Function
The aim of this study is to estimate the timewise varying heat flux applied at the time i t and the surface 1 Γ , ( ) i q t ,  1, i f using the transient readings of a single sensor S placed at the point ( , ) Si Sj inside the body.To do so, an inverse analysis is used so that the square of the difference between the estimated temperatures at the sensor place, ( , , ) T Si Sj t , computed from the solution of the direct heat conduction problem using the estimated heat flux ( ) i q t and the measured temperatures ( , , ) T Si Sj t over the time domain   0 f t t is minimized.This can be mathematically expressed as at Γ min : ( , , ) ( , , ) : Equation ( 1) in Ω, BCs and IC in Equations ( 2)-( 4) where q ( ), ( ), ( ),..., ( ) f q t q t q t q t T .The inverse analysis is used to minimize the following objective function expression:

Sensitivity Analysis
The inverse problem is concerned with the calculation of the gradient of the objective function J defined by Equation ( 17) with respect to ( ) i q t ,  1, i f .Thus, we can write In Equation ( 18), ( ) ) are called the sensitivity coefficients and can be explicitly expressed using the chain rule (using the constant thermal conductivity T k ) as follows The term in the numerator of Equation ( 19 , with respect to T k ,as follows The term in the denominator of Equation ( 19), , can be obtained from the boundary condition involving the applied heat flux, Equation ( 8), as follows where the terms ξ T , η T , J , γ , and β are computed using the finite-difference expressions associated with the surface 1 Γ .
Therefore, the expression in Equation ( 19) can be computed by dividing the expression in Equation (20) by the one in Equation ( 22).
Therefore, the sensitivity coefficients ( ) can be computed in only one single direct problem solution (during the transient solution) without the need for solving sensitivity and adjoint problems.The sensitivity matrix Ja can be explicitly written as T Si Sj t T Si Sj t T Si Sj t T Si Sj t q t q t q t q t T Si Sj t q t q t q t q t (24) However, the temperature estimated at any time is independent of a yet-to-occur future heat flux component [1,14] and the sensitivity matrix is a lower-triangular one.In other words, for   i i (the terms above the main diagonal of the sensitivity matrix), . Thus, we get T Si Sj t (25)

The Steepest-Descent Method
In this study, the steepest-descent optimization method is used to solve the inverse heat transfer problem.In this method, the objective function given by Equation ( 17) is minimized by searching along the direction of descent (k)  d using a search step length (k)  β .
The direction of descent is obtained from the gradient direction  (k) J , as follows The search step-length is given as follows Optimization Algorithm The following algorithm presents the direct and inverse analysis steps used to estimate the timewise varying heat flux applied at the boundary surface 1 Γ : 1. Specify the physical domain, the boundary and initial conditions, and the measured temperatures at the sensor place , Si Sj

S
and the time i t T Si Sj t .
2. 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 i t ( ( , , ) T Si Sj t , through solving Equations ( 6)-( 12).
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 ( ) q t Ja from Equation (25).10.From Equation (26), evaluate the new values for  q , namely   (k 1) q .
11. Set the next iteration ( k = k + 1 ) and return to step 2.

Stopping Criterion
Without measurement errors, the inverse problem can be terminated if 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 ε .For the temperature measurements with error, the discrepancy principle is used to terminate the iterative procedure.In the discrepancy principle, if the difference between estimated and measured temperatures is of the order of magnitude of the measurement errors, then the solution is assumed to be sufficiently accurate, that is, where σ is the standard deviation of the measurement errors and is assumed constant in this study (  0.1 σ ,  0.5 σ , and  1.0 σ ).We can obtain the following value for ε by substituting Equation (30) into Equation (17) (objective function definition) Then, the iterative procedure is terminated when the following criterion is satisfied here the measured temperatures containing random errors, meas ( , , ) where ω is a random variable with normal (Gaussian) 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.

Results
Three test cases are presented to investigate the accuracy and efficiency of the proposed sensitivity analysis method to estimate the timewise varying heat flux applied on part of the boundary of a heat conducting body.Initially, it is assumed that the heat flux is known, and the transient heat conduction problem is then solved to estimate the temperature at the sensor place at Then, the estimated temperatures are used as simulated measured ones to recover the initially used heat flux.Three different forms of timewise variation are considered for the heat flux, as follows 1) A step-function in Test case 1: 3) A triangular function in Test case 3:  .This implies that the number of unknown parameters is 10000 .Thus, the parameter estimation approach used so far in the literature is not feasible to estimate this extremely large number of unknown parameters.However, it will be shown that using the proposed sensitivity analysis, the estimation of such a large number of unknown parameters is feasible in an accurate and efficient manner.  ( , 5 2

M S
, used to measure the temperature at time i t .
In this study, three different measurement errors of  0.1,0.5,1.0 σ are considered.The stopping criteria (Equation (31)) for the test cases with these measurement errors are 1.0 1.0 (10000) 10000 σ ε σ f Since the size of the Jacobian matrix is  f f , we will deal here with a Jacobian matrix of size  10000 10000 .As soon as the nodal temperature at the sensor place is calculated at each time i t , the Jacobian matrix coefficients can be immediately calculated during the transient solution using the explicit expression for the sensitivity coefficients; that is, during the transient solution of the direct problem, the terms , )  ( ,1)

Ja i j T i Ja i j q j
This means that the proposed sensitivity analysis scheme is very efficient and does not contribute significantly to the computational cost of the solution.If we assume that the computational cost of the expression given by Equation ( 15) at any node ( , i j ) and time i t is equal to a , then the computational cost of the entire transient heat conduction solution for  M N nodes at the final time f t is approximately equal to    M N f a .On the other hand, the computational cost of Equation ( 23) for a single sensor placed at only one node, ( , ) S Si Sj , at the final time f t is approximately equal to   1 f a .Therefore, the computational cost of the computation of Jacobian matrix at each iteration of the computational cost of the direct problem solution.In this study, since , the computational cost of the computation of Jacobian matrix is of the computational cost of the direct problem solution.
The initial guesses used for the estimation of the unknown timewise varying heat fluxes in the test cases are as follows: Initially, the temperature distribution for the heat conducting body used in one of the test cases (Test case 2) is compared to the temperature distribution obtained from the commercial finite element software COMSOL to validate the implementation of the direct problem solution.To do so, the timewise varying heat flux given in Test case 2 (   5000 400 sin( ),0< 2000s 180 ) is considered.The grid and the temperature distribution obtained by the explicit solver, Equation (15), using the problem data given in Table 1, a grid size of  20 20 , and the time step   0.2s t are shown in Figure 3a,b  Three different functions (step, sinusoidal, and triangular function) are considered for the timewise variation of the applied heat flux to investigate the inverse analysis presented here.These different functions are selected so that they can reflect the accuracy and efficiency of the inverse analysis as they involve sharp comers and discontinuities, which are the most difficult to be recovered by an inverse analysis [2].Considering the three different functions, a comparison of the initial (guessed), optimal, and desired heat fluxes is shown in Figures 6a-8a  ).It can be seen that very accurate results are obtained for both cases of no measurement error and the measurement error, even for functions containing sharp comers and discontinuities.), respectively.Without the measurement error, a 100% reduction in the objective function and complete recovering of the unknown heat flux are achieved in all test cases, which shows the accuracy of the proposed sensitivity analysis scheme.Moreover, an acceptable recovering of unknown heat flux even in the presence of large measurement errors and a significant reduction in the objective function value implies that the proposed inverse analysis is not strongly affected by the errors involved in the temperature measurements.In addition, in all test cases with and without measurement error, a high convergence rate can be seen and an approximately 99% reduction in the objective function value takes place in the first ten iterations.As can be seen in the figures, the number of required iterations to obtain a stable solution using the discrepancy principle in recovering the unknown heat flux is decreased by increasing the measurement error.The only exception is the case of the measurement error of  0.1 σ in Test case 2 in which the number of iterations is 467 (larger than 401 for  0.0 σ ).The reason is that the iterative process for the case of the measurement error of  0.1 σ could be terminated much sooner because the discrepancy principle, Equation (30), is an approximation expression.The value of the objective function at iteration 314 is 101.007 and the value of the objective function at iteration 467 is 99.999.The objective function is decreased from 101.007 to 99.999 in 154 iterations.In other words, one could terminate the iterative process at iteration 314 instead of iteration 467 to get a stable solution.In the case of the measurement error, some oscillatory behaviors are observed around the exact values due to the illposed nature of the inverse heat transfer problem.The details of the results, including the initial and desired values for the variable heat flux, the initial and final values of the objective function, the computation time, the number of iterations, and the percentage of the decrease in the objective function, are given in  It can also be seen that in a neighborhood of f t the estimated heat flux deviates from the exact one and approaches the initial heat flux.One method to overcome this drawback and reduce the effects of the initial heat flux on the solution in the time interval of interest is to consider a final time larger than that of interest.From Equation (18) we can write the gradient of the objective function J with respect to  q as Ja T T q

J
(34) By approaching the final time f t , the elements with zero value in the column vectors of the Jacobian matrix Ja increase and there is only one nonzero element in the last column vector (see Equation ( 25)) because the Jacobian matrix is a lower-triangular matrix.For example, the gradient of the objective function J with respect to the heat flux at the final time, ( )  q t q t , as observed.

Conclusions
An explicit sensitivity analysis scheme was developed in this study to estimate the unknown functional form of the timewise varying heat flux applied on part of the boundary of a heatconducting body subjected to specified initial and boundary conditions through transient readings of a single sensor placed inside the body.There was no prior information available on the functional form of the unknown heat flux.Unlike the function estimation approach which is the widespread approach to address this problem, in this study a parameter estimation approach was used for the first time to estimate the unknown functional form accurately and efficiently through demonstrating three test cases.The transient heat conduction equation was solved using an explicit scheme.Initially, the two-dimensional irregular domain was transformed into a regular computational domain and all computations needed to solve the direct and inverse problem equations were carried out in the regular computational domain.The elliptic grid generation method was used to mesh the physical domain and the finite-difference method was used to approximate the derivatives of the field variable (temperature) at the grid nodes by algebraic ones.A novel, efficient, accurate, and very easy to implement sensitivity analysis scheme was developed by using the chain rule to obtain an explicit relation between the temperature at the sensor place and the timewise varying wall heat flux applied on the part of the domain boundary which allowed for the calculation of the sensitivity coefficients at each time step during the transient solution.Thus, the inverse analysis using the proposed sensitivity analysis scheme is not involved with sensitivity and adjoint equations.The main novelty is that all sensitivity coefficients can be computed efficiently in only one direct problem solution (during the transient solution), without the need for the solution of the sensitivity and adjoint equations (to compute the gradient of the objective function with respect to the variables).The steepest-descent method with an appropriate stopping criterion specified by the discrepancy principle was used to minimize the objective function and recover the unknown functional form accurately.The accuracy and efficiency of the proposed numerical procedure were demonstrated through presenting three test cases.

Figure 2 .
Figure 2. The place of the sensor S ,( , 5)  2 ,..., j f are computed from Equations (20) and (22), respectively, and then the sensitivity coefficients can be obtained from

Figure 3 .Figure 4 .
Figure 3. Validation of the direct problem solver using the finite element software COMSOL.The result obtained by using our explicit code (a,b), and the result obtained by using COMSOL (c).

Figure 6 .
Figure 6.Estimation of the timewise varying heat flux  step q

Figure 7 .Figure 8 .
Figure 7. Estimation of the timewise varying heat flux  sinusoidal q

Table 1 .
The heat-conducting body is made of stainless steel (type 304).The numerical values of the coefficients involved in the test cases are listed in Table1.Data used for the heat conduction problem (the body is made of stainless steel (type 304)).

Table 2 .
A summary of results for the estimation of the timewise varying heat flux