Estimation of Functional Form of Time-Dependent Heat Transfer Coefficient Using an Accurate and Robust Parameter Estimation Approach: An Inverse Analysis

This paper presents a numerical method to address function estimation problems in inverse heat transfer problems using parameter estimation approach without prior information on the functional form of the variable to be estimated. Using an inverse analysis, the functional form of a time-dependent heat transfer coefficient is estimated efficiently and accurately. The functional form of the heat transfer coefficient is assumed unknown and the inverse heat transfer problem should be treated using a function estimation approach by solving sensitivity and adjoint problems during the minimization process. Based on proposing a new sensitivity matrix, however, the functional form can be estimated in an accurate and very efficient manner using a parameter estimation approach without the need for solving the sensitivity and adjoint problems and imposing extra computational cost, mathematical complexity, and implementation efforts. In the proposed sensitivity analysis scheme, all sensitivity coefficients can be computed in only one direct problem solution at each iteration. In this inverse heat transfer problem, the body shape is irregular and meshed using a body-fitted grid generation method. The direct heat conduction problem is solved using the finite-difference method. The steepest-descent method is used as a minimization algorithm to minimize the defined objective function and the termination of the minimization process is carried out based on the discrepancy principle. A test case with three different functional forms and two different measurement errors is considered to show the accuracy and efficiency of the used inverse analysis.


Introduction
In a heat transfer problem, the accuracy of thermophysical properties and boundary conditions is critical to obtain an accurate numerical simulation. As a boundary condition, convective heat transfer depends on different parameters such as time, surface geometry, and surface temperature, to name a few. The accurate determination of the convective heat transfer coefficient is a difficult task as convection is a very complicated phenomenon and expensive experiments with sophisticated instruments are required to appropriately unravel its dynamics [1].
The advent of high-speed and high-capacity computers and the development of different regularization methods over the past decades have played significant roles in successful applications of numerical inverse methods, as inexpensive alternatives to costly and time-consuming experiments with sophisticated instruments, to appropriately estimate unknown heat transfer quantities such as the heat transfer coefficient [2][3][4][5][6][7][8][9][10]. Inverse heat transfer problems are mathematically challenging problems because they are ill-posed and the accuracy of the estimation of an unknown quantity is very sensitive to measurement errors [11][12][13]. If the unknown quantity to be estimated (in this study, the heat transfer coefficient) can be expressed as a constant parameter [14,15] or represented by a few parameters [16], a parameter estimation approach may be used to estimate the parameters thereby estimating the unknown quantity. However, a function estimation approach should be used to estimate the unknown functional form of the unknown quantity when there is no information available on the functional form of the unknown quantity. In the function estimation approach, the sensitivity and adjoint problems are required to obtain the gradient of objective function with respect to unknown functional form which impose additional mathematical developments and computational costs on the inverse analysis. In this study, based on the numerical procedure employed in [17], a two-dimensional transient inverse heat conduction problem is considered. The thermophysical properties are assumed constant, the geometry of heat-conducting body is irregular, and the body is subject to Neumann and Robin boundary conditions at its boundary surface parts. Using the parameter estimation approach initially developed in [17] for the estimation of unknown functional form of a time-dependent heat flux (a boundary condition of second kind) imposed at a boundary surface, here the unknown functional form of a time-dependent heat transfer coefficient(a third-kind boundary condition) is estimated efficiently and accurately without involving the solution of the sensitivity and adjoint problems. Thus, the mathematical development effort and the computational cost are reduced significantly.
To do so, the heat-conducting body (the physical domain) is mapped onto a regular computational domain in order to take advantage of the ease of implementation of the finite-difference method to explicitly solve the transient heat conduction equation and the associated boundary conditions. As the body shape is irregular, a body-fitted (elliptic) grid generation method is used to mesh the irregular domain which makes the proposed method general and applicable to any irregular domain as long as it can be mapped onto a regular computational domain. Using the chain rule to relate the temperature at sensor place and the time-dependent heat transfer coefficient applied on the part of the body boundary, explicit expressions are derived to compute sensitivity coefficients during the solution of the transient heat conduction equation without the need for solving the sensitivity and adjoint equations. The steepest-descent method, as an iterative regularization method, with a stopping criterion specified by discrepancy principle is used to minimize the objective function and reach the solution accurately. A test case with three complicated functional forms of timewise variation of the heat transfer coefficient is presented to reveal the accuracy, efficiency, and robustness of the inverse analysis. Moreover, two different measurement errors are considered. It is shown that the inverse analysis is not strongly affected by the errors involved in the temperature measurements and the unknown functional forms of the timewise variation of the heat transfer coefficient can be recovered with excellent accuracy. As stated before, the objective of this study is to present a parameter estimation approach to estimate the unknown functional form of a time-dependent heat transfer coefficient efficiently and accurately.

Governing Equation
The body shown in Figure 1a is initially at the temperature T 0 . At time t > 0, it is exposed to a time-dependent heat flux . q(t) at boundary surface Γ 1 and convective heat transfer on boundary surfaces Γ i , i = 2, 3, 4 with corresponding heat transfer coefficients h 2 (t), h 3 , and h 4 and surrounding temperatures T ∞ i , i = 2, 3, 4. The thermal conductivity, density, and specific heat of the body are k T , ρ, and c, respectively.  The governing equation for a two-dimensional transient heat conduction problem with no heat generation can be expressed as [17,18] ae with the boundary and initial conditions where t is the time. Since the heat-conducting body is irregular, it (the x and y physical domain) can be mapped onto a regular one (the ξ and η computational domain). The elliptic grid generation method is employed to generate a grid over the physical 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 [12,18]. The transformation results in , then a smooth grid over the physical domain is obtained. Therefore, Equation (6)  The governing equation for a two-dimensional transient heat conduction problem with no heat generation can be expressed as [17,18] with the boundary and initial conditions ∂T(x, y, t) ∂T(x, y, t) T(x, y, 0) = T 0 (x, y) in physical domain Ω(x, y) where t is the time. Since the heat-conducting body is irregular, it (the x and y physical domain) can be mapped onto a regular one (the ξ and η computational domain). The elliptic grid generation method is employed to generate a grid over the physical 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 [12,18]. The transformation results in where ∇ 2 ξ = P(ξ, η) and ∇ 2 η = Q(ξ, η) are grid control functions. If P(ξ, η) = Q(ξ, η) = 0, then a smooth grid over the physical domain is obtained. Therefore, Equation (6) becomes where The transformed boundary and initial conditions can be expressed as where the initial condition T 0 (x, y) is rewritten as T * 0 (ξ, η) in terms of the variables ξ and η. Now the finite-difference method can be employed to discretize the derivatives present in the above equations in the regular computational domain, as follows (assuming ∆ξ = ∆η = 1) where f ≡ x, y, T. One-sided forward and one-sided backward relations are used to discretize the boundary condition equations. The explicit method can be used to solve the resulting transient heat conduction equation, Equation (7). Using forward-time-centralspace (FTCS) discretization and the relations in Equation (14), we get where ∆t is the time step. Taking into account the stability criterion, the time-marching procedure can be used to solve Equation (15) and obtain T n+1 i,j . That is, 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 Energies 2021, 14, 5073 5 of 18

Objective Function
The inverse heat transfer problem of interest deals with the estimation of the timedependent heat transfer coefficient h 2 (t) applied at the time t i and at the surface Γ 2 , h 2 (t i ), i = 1, . . . , r (r is the number of time steps) using the transient readings of a single sensor S placed at the point (Si, Sj) inside the heat conducting body. In inverse analysis, the aim is to minimize the mismatch between the estimated temperatures at the sensor place, T e (Si, Sj, t i ), computed from the solution of the direct transient heat conduction problem using the estimated heat transfer coefficient h 2 (t i ) and the measured temperatures T m (Si, Sj, t i ) over the time domain 0 < t < t r . This can be mathematically expressed as a least-squares minimization as follows where Therefore, the objective function can be expressed as

Sensitivity Analysis
The calculation of the gradient of the objective function J defined by Equation (18) with respect to h 2 (t i ), i = 1, . . . , r is required in gradient-based minimization methods. Therefore, it can be written The sensitivity coefficients . . , r, i = 1, . . . , r) can be explicitly expressed using the chain rule (using the constant thermal conductivity k T ) as follows The expression in the numerator of Equation (20), , can be obtained by taking in Equation (16) with respect to k T , as follows The expression in the denominator of Equation (20), , can be obtained from the boundary condition involving the heat transfer coefficient h 2 (t), Equation (9), as follows therefore, we get where the terms T ξ , T η , J, γ and β are computed using the finite-difference expressions associated with the surface Γ 2 . That is, where f ≡ x, y, T. Therefore, the sensitivity coefficients in Equation (20) can be calculated by dividing the term in Equation (21) by the one in Equation (23).
It can be seen that all sensitivity coefficients can be computed during the transient solution of the direct heat conduction equation. The sensitivity matrix Ja can be explicitly written as . . .
We know that the temperature estimated at any time is independent of a yet-to-occur future heat transfer coefficient component [11,19] which gives rise to a lower-triangular sensitivity matrix. That is, for i > i (the terms above the main diagonal of the sensitivity matrix), . . .

The Steepest-Descent Method
The steepest-descent optimization method is used in this study to minimize the objective function given by Equation (18) by searching along the direction of steepest descent d (k) using a search step length β (k) .
The negative of the gradient direction ∇J (k) denotes the direction of steepest descent d (k) . Thus, , the search steplength β (k) > 0 is given as follows [12] Optimization Algorithm The proposed numerical procedure to estimate the time-dependent heat transfer coefficient at the boundary surface Γ 2 , h 2 (t), can be summarized as follows: 1.
Measure the temperatures at the sensor place S Si,Sj and the time t i (i = 1, . . . , r), Solve the direct problem to obtain the temperature values at the sensor place and the time t i (i = 1, . . . , r), T e (Si, Sj, t i ), through solving Equations (7)-(13). 3.

4.
If value of the objective function obtained in step 3 is less than the specified stopping criterion, the optimization is finished. Otherwise, go to step 5.
Update h 2 from Equation (27) and set the next iteration (k = k + 1) and return to the step 2.

Stopping Criterion
The minimization process can be terminated if where ε is chosen based on obtaining stable and appropriate results. In this study, when there are no measurement errors, ε = 10 −6 . However, if the temperature measurements contain errors, then the discrepancy principle is used to stop the iterative procedure and obtain stable results. In this case, where σ is the standard deviation of the measurement errors [17] and is assumed constant in this study (σ = 0.001 and σ = 0.005). Here, the measured temperatures containing random errors, T meas (Si, Sj, t i ), (i = 1, . . . , r), are generated by adding an error term ωσ to the exact temperatures T exact (Si, Sj, t i ) to give 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
A test case with three different complicated functional forms of variation of the heat transfer coefficient with time is presented to investigate the accuracy, efficiency, and robustness of the proposed sensitivity analysis method to estimate the time-dependent heat transfer coefficient on part of the boundary of a heat conducting body. Initially the heat transfer coefficient is assumed to be known, the transient heat conduction problem is then solved to calculate the temperature at the sensor place at times t i (i = 1, . . . , r). Then, the calculated temperatures are used as simulated measured ones to recover the initially used heat transfer coefficient. The three different forms of timewise variation of the heat transfer coefficient considered are as follows (Figures 2-4) h (1) ì ï ï ï ï ï ï ï ï ï ï ï ï í ï ï ï ï ï ï ï < £ ï ï ï ï < £ ï î = + + Assuming the heat conducting body is made of stainless steel (type 304), the numerical values of the coefficients involved in the test case are listed in Table 1.   ì ï ï ï ï ï ï ï ï ï ï ï ï í ï ï ï ï ï ï ï < £ ï ï ï ï < £ ï î = + + Assuming the heat conducting body is made of stainless steel (type 304), the numerical values of the coefficients involved in the test case are listed in Table 1.     Assuming the heat conducting body is made of stainless steel (type 304), the numerical values of the coefficients involved in the test case are listed in Table 1. Table 1. Data used for the heat conduction problem (the body is made of stainless steel (type 304) [20]).
2 (t) 5 30 In all simulations in this study, the heat conducting body is meshed using a grid size of M × N = 40 × 40, the temperature measurement sensor is placed at node (Si, Sj) = ( M 2 , N − 3) = (20, 37) (close to the boundary subject to convective heat transfer with the convective heat transfer coefficient h 2 to obtain sensible sensitivity coefficients) ( Figure 5), the initial temperature is T(x, y, 0) = T 0 (x, y) = 20 • C, the final time is t r = 600 s, and the time step is ∆t = 0.1 s. Thus, the number of transient readings of the single sensor S is r = t r ∆t = 600 s 0.1 s = 6000. This means that that the number of unknown parameters is 6000. Thus, the estimation of such a large number of unknown parameters using the parameter estimation approach commonly used in the literature is not feasible. However, using the proposed sensitivity analysis, one can handle the estimation of the large number of unknown parameters accurately and efficiently. In this study, two different measurement errors of σ = 0.001 and σ = 0.005 are considered. The stopping criteria for the test case with the following measurement errors are σ = 0.001 ⇒ ε = σ 2 r = 0.001 2 (6000) = 0.006 σ = 0.005 ⇒ ε = σ 2 r = 0.005 2 (6000) = 0.15 As the size of Jacobian matrix is r × r, we will deal in this study with a Jacobian matrix of size 6000 × 6000. Once the temperature at the sensor place is obtained at the time t i , the elements of the Jacobian matrix can be calculated during the transient solution using the obtained expression for the sensitivity coefficients; that is, during the solution of the direct problem, the terms T k T (i, 1) = ∂T e (Si, Sj, t i ) ∂k T , i = 1, . . . , r and h 2 k T (j, 1) = ∂h 2 (t j ) ∂k T , j = 1, . . . , r are computed from Equations (21) and (23), respectively, and then the sensitivity coefficients can be obtained using the following pseudo-code Initially, the implementation of the direct problem solver is validated with the results obtained from the commercial finite element software COMSOL. To do so, using the data given in Table 1, h (2) 2 , and the body shown in Figure 5, the temperature distribution in the body is calculated by the two methods (our finite-difference explicit code, Equation (16), using two different time steps of 0.1 and 0.001 s and the finite element software COMSOL) which is shown in Figure 6. Moreover, the temperature history of the place of the sensor, , obtained by both methods is shown in Figure 7. The comparison between the results reveals a very good agreement hereby verifying the correct implementation of the explicit solver. ¶ = ¶ Initially, the implementation of the direct problem solver is validated with the re sults obtained from the commercial finite element software COMSOL. To do so, using th data given in Table 1, (2) 2 h , and the body shown in Figure 5, the temperature distribution in the body is calculated by the two methods (our finite-difference explicit code Equation (16)   In this inverse heat conduction problem, three different and complicated functional forms of timewise variations (including the variations which are difficult to be recovered by an inverse analysis such as discontinuities and sharp corners [12]) for the heat transfer coefficient are chosen to examine the accuracy, efficiency, and robustness of the inverse analysis presented in this study. A comparison of the initial (guessed), final, and desired heat transfer coefficients is shown in Figures 8a-10a   In this inverse heat conduction problem, three different and complicated functional forms of timewise variations (including the variations which are difficult to be recovered by an inverse analysis such as discontinuities and sharp corners [12]) for the heat transfer coefficient are chosen to examine the accuracy, efficiency, and robustness of the inverse analysis presented in this study. A comparison of the initial (guessed), final, and desired heat transfer coefficients is shown in Figures 8a, 9a and 10a      From the above figures, we can also see that the estimated heat transfer coefficient deviates from the exact one in a neighborhood of t r and approaches the initially guessed heat transfer coefficient. The mathematical reason is that by approaching the final time t r , the number of the zero elements in the column vectors of the sensitivity matrix Ja also increases so that there exists only one nonzero element in the last column vector because the sensitivity matrix is a lower-triangular matrix (see Equation (26)). Thus, the last column vector can be written as  From Equation (19), , we can write the gradient of the objective function J with respect to h 2 at the final time t r as which is a very small number. Substituting a very small value for ∇J there is no significant modification in the value of h 2 during the minimization process and the heat transfer coefficient retains its initially guessed value until the end of the minimization process.

Conclusions
Based on explicit expressions, an accurate and efficient sensitivity analysis scheme was proposed to estimate the unknown functional form of a time-dependent heat transfer coefficient applied on part of the boundary of an irregular heat conducting body subjected to specified initial and boundary conditions through transient readings of a single sensor placed inside the irregular body. Since there was no prior information available on the functional form of the variable to be estimated, the commonly used method to address this inverse heat conduction problem was based on the function estimation approach. However, here a parameter estimation approach was proposed to estimate the unknown functional form accurately and efficiently. As the body geometry was general (irregular), the physical domain was mapped onto a regular computational domain in order to take advantage of the ease of implementation of the finite-difference method to explicitly solve the transient heat conduction equation and associated boundary conditions. The chain rule was used to relate the sensor temperature and the time-dependent heat transfer coefficient applied on the part of the body boundary, the two ingredients of the sensitivity coefficients. Formulating this way, all sensitivity coefficients could be computed during the transient solution of the direct heat conduction problem without the need for solving the sensitivity and adjoint equations. The steepest-descent method with a stopping criterion specified by the discrepancy principle was used to minimize the objective function and reach the solution accurately. The accuracy, efficiency, and robustness of the inverse analysis were presented through considering three different complicated functional forms. As a future study, more challenging problems of heat conduction in materials with temperaturedependent or space-dependent thermal conductivity (in functionally graded materials, for instance) may be considered. In this study, the transient heat transfer equation (direct solution) was solved by the explicit method; the feasibility of derivation of the sensitivity coefficients using an implicit method needs to be investigated.