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

Abstract

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.

1. 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.

2. 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 $\dot{q}(t)$ at boundary surface ${\Gamma }_1$ and convective heat transfer on boundary surfaces ${\Gamma }_i,i=2,3,4$ with corresponding heat transfer coefficients $h_2(t)$, $h_3$, and $h_4$ and surrounding temperatures $T_{{\infty }_i},i=2,3,4$. The thermal conductivity, density, and specific heat of the body are $k_T$, $\rho$, and $c$, respectively.

The governing equation for a two-dimensional transient heat conduction problem with no heat generation can be expressed as [17,18]

$$k_T\left(\frac{{\partial }^{2}T(x,y,t)}{\partial x^2}+\frac{{\partial }^{2}T(x,y,t)}{\partial y^2}\right)=\rho c\frac{\partial T(x,y,t)}{\partial t}\mathrm{in}\mathrm{physical}\mathrm{domain}\Omega (x,y)$$

(1)

with the boundary and initial conditions

$$\frac{\partial T(x,y,t)}{\partial n_1}=\frac{\dot{q}(t)}{k_T}\mathrm{on}\mathrm{boundary}\mathrm{surface}{\Gamma }_1(x,y)$$

(2)

$$\frac{\partial T(x,y,t)}{\partial n_2}=-\frac{h_2(t)}{k_T}\left(T_{{\Gamma }_2}(x,y,t)-T_{{\infty }_2}\right)\mathrm{on}\mathrm{boundary}\mathrm{surface}{\Gamma }_2(x,y)$$

(3)

$$\frac{\partial T(x,y,t)}{\partial n_i}=-\frac{h_i}{k_T}\left(T_{{\Gamma }_i}(x,y,t)-T_{{\infty }_i}\right)\mathrm{on}\mathrm{boundary}\mathrm{surface}{\Gamma }_i(x,y),i=3,4$$

(4)

$$T(x,y,0)=T_0(x,y)\mathrm{in}\mathrm{physical}\mathrm{domain}\Omega (x,y)$$

(5)

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 $\xi$ and $\eta$ 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 ($\xi ,\eta ,t$) variables [12,18]. The transformation results in

$$\left(\frac{\alpha \frac{{\partial }^{2}T(\xi ,\eta ,t)}{\partial {\xi }^2}-2\beta \frac{{\partial }^{2}T(\xi ,\eta ,t)}{\partial \xi \partial \eta }+\gamma \frac{{\partial }^{2}T(\xi ,\eta ,t)}{\partial {\eta }^2}}{J^2}+({\nabla }^2\xi )\frac{\partial T(\xi ,\eta ,t)}{\partial \xi }+({\nabla }^2\eta )\frac{\partial T(\xi ,\eta ,t)}{\partial \eta }\right)=\frac{\rho c}{k_T}\frac{\partial T(\xi ,\eta ,t)}{\partial t}$$

(6)

where ${\nabla }^2\xi =P(\xi ,\eta )$ and ${\nabla }^2\eta =Q(\xi ,\eta )$ are grid control functions. If $P(\xi ,\eta )=Q(\xi ,\eta )=0$, then a smooth grid over the physical domain is obtained. Therefore, Equation (6) becomes

$$\left(\frac{\alpha \frac{{\partial }^{2}T(\xi ,\eta ,t)}{\partial {\xi }^2}-2\beta \frac{{\partial }^{2}T(\xi ,\eta ,t)}{\partial \xi \partial \eta }+\gamma \frac{{\partial }^{2}T(\xi ,\eta ,t)}{\partial {\eta }^2}}{J^2}\right)=\frac{\rho c}{k_T}\frac{\partial T(\xi ,\eta ,t)}{\partial t}\mathrm{in}1<\xi <M,1<\eta <N,\mathrm{for}t>0$$

(7)

where

$$\alpha =x_{\eta }^2+y_{\eta }^2$$

$$\beta =x_{\xi }{x}_{\eta }+y_{\xi }{y}_{\eta }$$

$$\gamma =x_{\xi }^2+y_{\xi }^2$$

$$J=x_{\xi }{y}_{\eta }-x_{\eta }{y}_{\xi }\left(\mathrm{Jacobian}\mathrm{of}\mathrm{transformation}\right)$$

(8)

The transformed boundary and initial conditions can be expressed as

$${\left(\frac{-1}{J\sqrt{\gamma }}\left(\gamma \frac{\partial T(\xi ,\eta ,t)}{\partial \eta }-\beta \frac{\partial T(\xi ,\eta ,t)}{\partial \xi }\right)\right)}_{{\Gamma }_1}=\frac{\dot{q}(t)}{k_T}\mathrm{at}1<\xi <M,\eta =1,\mathrm{for}t>0$$

(9)

$${\left(\frac{1}{J\sqrt{\gamma }}\left(\gamma \frac{\partial T(\xi ,\eta ,t)}{\partial \eta }-\beta \frac{\partial T(\xi ,\eta ,t)}{\partial \xi }\right)\right)}_{{\Gamma }_2}=-\frac{h_2(t)}{k_T}\left(T(\xi ,\eta ,t)-T_{{\infty }_2}\right)\mathrm{at}1<\xi <M,\eta =N,\mathrm{for}t>0$$

(10)

$${\left(\frac{-1}{J\sqrt{\alpha }}\left(\alpha \frac{\partial T(\xi ,\eta ,t)}{\partial \xi }-\beta \frac{\partial T(\xi ,\eta ,t)}{\partial \eta }\right)\right)}_{{\Gamma }_3}=-\frac{h_3}{k_T}\left(T(\xi ,\eta ,t)-T_{{\infty }_3}\right)\mathrm{at}1<\eta <N,\xi =1,\mathrm{for}t>0$$

(11)

$${\left(\frac{1}{J\sqrt{\alpha }}\left(\alpha \frac{\partial T(\xi ,\eta ,t)}{\partial \xi }-\beta \frac{\partial T(\xi ,\eta ,t)}{\partial \eta }\right)\right)}_{{\Gamma }_4}=-\frac{h_4}{k_T}\left(T(\xi ,\eta ,t)-T_{{\infty }_4}\right)\mathrm{at}1<\eta <N,\xi =M,\mathrm{for}t>0$$

(12)

$$T(\xi ,\eta ,0)=T_0^{*}(\xi ,\eta )\mathrm{in}1<\xi <M,1<\eta <N,\mathrm{for}t=0$$

(13)

where the initial condition $T_0(x,y)$ is rewritten as $T_0^{*}(\xi ,\eta )$ in terms of the variables $\xi$ and $\eta$. 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 $\Delta \xi =\Delta \eta =1$)

$$f_{\xi }=\frac{1}{2}(f_{i+1,j}-f_{i-1,j})$$

$$f_{\eta }=\frac{1}{2}(f_{i,j+1}-f_{i,j-1})$$

$$f_{\xi \xi }=f_{i+1,j}-2f_{i,j}+f_{i-1,j}$$

$$f_{\eta \eta }=f_{i,j+1}-2f_{i,j}+f_{i,j-1}$$

$$f_{\xi \eta }=\frac{1}{4}(f_{i+1,j+1}-f_{i-1,j+1}-f_{i+1,j-1}+f_{i-1,j-1})$$

(14)

where $f\equiv 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-central-space (FTCS) discretization and the relations in Equation (14), we get

$$\begin{gathered} \frac{1}{J^2}\left(\alpha (T_{i+1,j}^n-2T_{i,j}^n+T_{i-1,j}^n)-2\beta \frac{1}{4}(T_{i+1,j+1}^n-T_{i-1,j+1}^n-T_{i+1,j-1}^n+T_{i-1,j-1}^n)+\gamma (T_{i,j+1}^n-2T_{i,j}^n+T_{i,j-1}^n)\right) \\ =\frac{\rho c}{k_T}\frac{T_{i,j}^{n+1}-T_{i,j}^n}{\Delta t},i=2,\dots ,M-1,j=2,\dots ,N-1\mathrm{for}t>0 \end{gathered}$$

(15)

where $\Delta 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_{i,j}^{n+1}$. That is, the nodal temperatures at the time level $n+1$, $T_{i,j}^{n+1}$, can be determined from the knowledge of nodal temperatures at the previous time level $n$, $T_{i,j}^n$, as follows

$$\begin{gathered} T_{i,j}^{n+1}=T_{i,j}^n+ \\ \frac{k_T\Delta t}{\rho c{J}^2}\left(\alpha (T_{i+1,j}^n-2T_{i,j}^n+T_{i-1,j}^n)-2\beta \frac{1}{4}(T_{i+1,j+1}^n-T_{i-1,j+1}^n-T_{i+1,j-1}^n+T_{i-1,j-1}^n)+\gamma (T_{i,j+1}^n-2T_{i,j}^n+T_{i,j-1}^n)\right) \end{gathered}$$

(16)

3. The Inverse Analysis

3.1. Objective Function

The inverse heat transfer problem of interest deals with the estimation of the time-dependent heat transfer coefficient $h_2(t)$ applied at the time $t_i$ and at the surface ${\Gamma }_2$, $h_2(t_i)$, $i=1,\dots ,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

$$\min \left\{\underset{{\mathbf{h}}_2\mathrm{at}{\Gamma }_2}{\mathrm{J}}:={\Vert {\mathbf{T}}_e(Si,Sj,t)-{\mathbf{T}}_m(Si,Sj,t)\Vert }^2:\mathrm{Equation}\hspace{0.17em}(1)\mathrm{in}\Omega ,\mathrm{BCs}\mathrm{and}\mathrm{IC}\mathrm{in}\mathrm{Equations}\left(2\right)–\left(5\right)\right\}$$

(17)

where ${\mathbf{h}}_2={\left[h_2(t_1),h_2(t_2),h_2(t_3),\dots ,h_2(t_r)\right]}^{\mathbf{T}}$. Therefore, the objective function can be expressed as

$$\mathrm{J}={\displaystyle \sum _{i=1}^r{\left[T_e(Si,Sj,t_i)-T_m(Si,Sj,t_i)\right]}^2}$$

(18)

3.2. Sensitivity Analysis

The calculation of the gradient of the objective function $\mathrm{J}$ defined by Equation (18) with respect to $h_2(t_i)$, $i=1,\dots ,r$ is required in gradient-based minimization methods. Therefore, it can be written

$$\frac{\partial \mathrm{J}}{\partial h_2(t_{i^{\prime }})}=2{\displaystyle \sum _{i=1}^r\left[T_e(Si,Sj,t_i)-T_m(Si,Sj,t_i)\right]}\frac{\partial T_e(Si,Sj,t_i)}{\partial h_2(t_{i^{\prime }})}$$

(19)

The sensitivity coefficients $\frac{\partial T_e(Si,Sj,t_i)}{\partial h_2(t_{i^{\prime }})}$ ($i=1,\dots ,r$, $i^{\prime }=1,\dots ,r$) can be explicitly expressed using the chain rule (using the constant thermal conductivity $k_T$) as follows

$$\frac{\partial T_e(Si,Sj,t_i)}{\partial h_2(t_{i^{\prime }})}=\frac{\frac{\partial T_e(Si,Sj,t_i)}{\partial k_T}}{\frac{\partial h_2(t_{i^{\prime }})}{\partial k_T}}$$

(20)

The expression in the numerator of Equation (20), $\frac{\partial T_e(Si,Sj,t_i)}{\partial k_T}$, can be obtained by taking derivative of $T_{i,j}^{n+1}$ in Equation (16) with respect to $k_T$, as follows

$$\begin{gathered} \frac{\partial T_e^{n+1}(Si,Sj,t_i)}{\partial k_T}=\frac{\Delta t}{\rho c{J}^2}(\alpha (T_{Si+1,Sj}^n-2T_{Si,Sj}^n+T_{Si-1,Sj}^n) \\ -2\beta \frac{1}{4}(T_{Si+1,Sj+1}^n-T_{Si-1,Sj+1}^n-T_{Si+1,Sj-1}^n+T_{Si-1,Sj-1}^n)+\gamma (T_{Si,Sj+1}^n-2T_{Si,Sj}^n+T_{Si,Sj-1}^n)) \end{gathered}$$

(21)

The expression in the denominator of Equation (20), $\frac{\partial T_e(Si,Sj,t_i)}{\partial k_T}$, can be obtained from the boundary condition involving the heat transfer coefficient $h_2(t)$, Equation (9), as follows

$${\left(\frac{1}{J\sqrt{\gamma }}\left(\gamma \frac{\partial T(\xi ,\eta ,t)}{\partial \eta }-\beta \frac{\partial T(\xi ,\eta ,t)}{\partial \xi }\right)\right)}_{{\Gamma }_2}=-\frac{h_2(t)}{k_T}\left(T(\xi ,\eta ,t)-T_{{\infty }_2}\right)$$

(22)

therefore, we get

$$\frac{\partial h_2(t_{i^{\prime }})}{\partial k_T}=\frac{-1}{T(\xi ,\eta ,t)-T_{{\infty }_2}}{\left(\frac{1}{J\sqrt{\gamma }}\left(\gamma \frac{\partial T(\xi ,\eta ,t)}{\partial \eta }-\beta \frac{\partial T(\xi ,\eta ,t)}{\partial \xi }\right)\right)}_{{\Gamma }_2}$$

(23)

where the terms $T_{\xi }$, $T_{\eta }$, $J$, $\gamma$ and $\beta$ are computed using the finite-difference expressions associated with the surface ${\Gamma }_2$. That is,

$$f_{\xi }=\frac{1}{2}(f_{i+1,N}-f_{i-1,N})$$

$$f_{\eta }=\frac{1}{2}(3f_{i,N}-4f_{i,N-1}+f_{i,N-2})$$

where $f\equiv 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).

$$\begin{gathered} \frac{\partial T_e(Si,Sj,t_i)}{\partial h_2(t_{i^{\prime }})}={\frac{1}{\left(\frac{-1}{T(\xi ,\eta ,t)-T_{{\infty }_2}}{\left(\frac{1}{J\sqrt{\gamma }}\left(\gamma \frac{\partial T(\xi ,\eta ,t)}{\partial \eta }-\beta \frac{\partial T(\xi ,\eta ,t)}{\partial \xi }\right)\right)}_{{\Gamma }_2}\right)}}_{t_{i^{\prime }}}\frac{\Delta t}{\rho c{J}^2}(\alpha (T_{Si+1,Sj}^{t_i}-2T_{Si,Sj}^{t_i}+T_{Si-1,Sj}^{t_i}) \\ -2\beta \frac{1}{4}(T_{Si+1,Sj+1}^{t_i}-T_{Si-1,Sj+1}^{t_i}-T_{Si+1,Sj-1}^{t_i}+T_{Si-1,Sj-1}^{t_i})+\gamma (T_{Si,Sj+1}^{t_i}-2T_{Si,Sj}^{t_i}+T_{Si,Sj-1}^{t_i})) \end{gathered}$$

(24)

It can be seen that all sensitivity coefficients $\frac{\partial T_e(Si,Sj,t_i)}{\partial h_2(t_{i^{\prime }})}$ can be computed during the transient solution of the direct heat conduction equation. The sensitivity matrix $\mathbf{Ja}$ can be explicitly written as

$$\mathbf{J}{\mathbf{a}}_{h_2(t)}={\left[\begin{array}{cccc}\frac{\partial T_e(Si,Sj,t_1)}{\partial h_2(t_1)}& \frac{\partial T_e(Si,Sj,t_1)}{\partial h_2(t_2)}& \frac{\partial T_e(Si,Sj,t_1)}{\partial h_2(t_3)}& \begin{array}{cc}\dots & \frac{\partial T_e(Si,Sj,t_1)}{\partial h_2(t_r)}\end{array}\\ \frac{\partial T_e(Si,Sj,t_2)}{\partial h_2(t_1)}& \frac{\partial T_e(Si,Sj,t_2)}{\partial h_2(t_2)}& \frac{\partial T_e(Si,Sj,t_2)}{\partial h_2(t_3)}& \begin{array}{cc}\dots & \frac{\partial T_e(Si,Sj,t_2)}{\partial h_2(t_r)}\end{array}\\ \frac{\partial T_e(Si,Sj,t_3)}{\partial h_2(t_1)}& \frac{\partial T_e(Si,Sj,t_3)}{\partial h_2(t_2)}& \frac{\partial T_e(Si,Sj,t_3)}{\partial h_2(t_3)}& \begin{array}{cc}\dots & \frac{\partial T_e(Si,Sj,t_3)}{\partial h_2(t_r)}\end{array}\\ \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_1)}\end{array}& \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_2)}\end{array}& \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_3)}\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \dots \end{array}& \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_r)}\end{array}\end{array}\end{array}\right]}_{r\times r}$$

(25)

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^{\prime }>i$ (the terms above the main diagonal of the sensitivity matrix), $\frac{\partial T_e(Si,Sj,t_i)}{\partial h_2(t_{i^{\prime }})}=0$. Thus, we get

$$\mathbf{J}{\mathbf{a}}_{h_2(t)}={\left[\begin{array}{cccc}\frac{\partial T_e(Si,Sj,t_1)}{\partial h_2(t_1)}& 0& 0& \begin{array}{cc}\dots & 0\end{array}\\ \frac{\partial T_e(Si,Sj,t_2)}{\partial h_2(t_1)}& \frac{\partial T_e(Si,Sj,t_2)}{\partial h_2(t_2)}& 0& \begin{array}{cc}\dots & 0\end{array}\\ \frac{\partial T_e(Si,Sj,t_3)}{\partial h_2(t_1)}& \frac{\partial T_e(Si,Sj,t_3)}{\partial h_2(t_2)}& \frac{\partial T_e(Si,Sj,t_3)}{\partial h_2(t_3)}& \begin{array}{cc}\dots & 0\end{array}\\ \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_1)}\end{array}& \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_2)}\end{array}& \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_3)}\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \dots \end{array}& \begin{array}{c}⋮\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_r)}\end{array}\end{array}\end{array}\right]}_{r\times r}$$

(26)

3.3. 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 ${\mathbf{d}}^{\left(\mathrm{k}\right)}$ using a search step length ${\beta }^{\left(\mathrm{k}\right)}$.

$${\mathbf{h}}_2{}^{(\mathrm{k}+1)}={\mathbf{h}}_2{}^{\left(\mathrm{k}\right)}+{\beta }^{\left(\mathrm{k}\right)}{\mathbf{d}}^{\left(\mathrm{k}\right)}$$

(27)

The negative of the gradient direction $\nabla {\mathrm{J}}^{\left(\mathrm{k}\right)}$ denotes the direction of steepest descent ${\mathbf{d}}^{\left(\mathrm{k}\right)}$. Thus,

$${\mathbf{d}}^{\left(\mathrm{k}\right)}=-\nabla {\mathrm{J}}^{\left(\mathrm{k}\right)}$$

(28)

In the expression ${\mathbf{h}}_2{}^{(\mathrm{k}+1)}={\mathbf{h}}_2{}^{\left(\mathrm{k}\right)}+{\beta }^{\left(\mathrm{k}\right)}{\mathbf{d}}^{\left(\mathrm{k}\right)}={\mathbf{h}}_2{}^{\left(\mathrm{k}\right)}-{\beta }^{\left(\mathrm{k}\right)}\nabla {\mathrm{J}}^{\left(\mathrm{k}\right)}$, the search step-length ${\beta }^{\left(\mathrm{k}\right)}>0$ is given as follows [12]

$${\beta }^{\left(\mathrm{k}\right)}=-\frac{{[\mathbf{J}{\mathbf{a}}^{\left(\mathrm{k}\right)}{\mathbf{d}}^{\left(\mathrm{k}\right)}]}^{\mathrm{T}}[{\mathbf{T}}_e-{\mathbf{T}}_m]}{{[\mathbf{J}{\mathbf{a}}^{\left(\mathrm{k}\right)}{\mathbf{d}}^{\left(\mathrm{k}\right)}]}^{\mathrm{T}}[\mathbf{J}{\mathbf{a}}^{\left(\mathrm{k}\right)}{\mathbf{d}}^{\left(\mathrm{k}\right)}]}$$

(29)

Optimization Algorithm

The proposed numerical procedure to estimate the time-dependent heat transfer coefficient at the boundary surface ${\Gamma }_2$, $h_2(t)$, can be summarized as follows:

• Measure the temperatures at the sensor place $S_{Si,Sj}$ and the time $t_i$ ($i=1,\dots ,r$), $T_m(Si,Sj,t_i)$.
• Solve the direct problem to obtain the temperature values at the sensor place and the time $t_i$ ($i=1,\dots ,r$), $T_e(Si,Sj,t_i)$, through solving Equations (7)–(13).
• Compute the objective function value (${\mathrm{J}}^{\left(\mathrm{k}\right)}$) using Equation (18).
• 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.
• Compute the sensitivity matrix $\mathbf{J}{\mathbf{a}}_{h_2(t)}$ from Equation (26), the gradient direction $\nabla {\mathrm{J}}^{(\mathrm{k})}$ from Equation (19), the direction of descent ${\mathbf{d}}^{\left(\mathrm{k}\right)}$ from Equation (28), and the search step-length ${\beta }^{\left(\mathrm{k}\right)}$ from Equation (29).
• Update ${\mathbf{h}}_2$ from Equation (27) and set the next iteration ($\mathrm{k}=\mathrm{k}+1$) and return to the step 2.

3.4. Stopping Criterion

The minimization process can be terminated if

$${\mathrm{J}}^{(\mathrm{k})}<\epsilon$$

(30)

where $\epsilon$ is chosen based on obtaining stable and appropriate results. In this study, when there are no measurement errors, $\epsilon ={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,

$$\epsilon =r{\sigma }^2$$

(31)

where $\sigma$ is the standard deviation of the measurement errors [17] and is assumed constant in this study ($\sigma =0.001$ and $\sigma =0.005$). Here, the measured temperatures containing random errors, $T_{\mathrm{meas}}(Si,Sj,t_i)$, $(i=1,\dots ,r)$, are generated by adding an error term $\omega \sigma$ to the exact temperatures $T_{\mathrm{exact}}(Si,Sj,t_i)$ to give

$$T_{\mathrm{meas}}(Si,Sj,t_i)=T_{\mathrm{exact}}(Si,Sj,t_i)+\omega \sigma$$

(32)

where $\omega$ is a random variable with normal (Gaussian) distribution, zero mean, and unitary standard deviation. Assuming 99% confidence for the measured temperature, $\omega$ lies in the range $-2.576\le \omega \le 2.576$ and it is randomly generated by using MATLAB.

4. 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,\dots ,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 (Figure 2, Figure 3 and Figure 4)

$$h_2^{(1)}(t)=\{\begin{array}{ll}12,& 0\le t<100\hspace{0.17em}\mathrm{s}\\ 8,& 100\le t<200\hspace{0.17em}\mathrm{s}\\ \frac{16-8}{300-200}(t-200)+8,& 200\le t\le 300\hspace{0.17em}\mathrm{s}\\ \frac{10-16}{400-300}(t-300)+16,& 300<t\le 400\hspace{0.17em}\mathrm{s}\\ 10,& 400<t\le 500\hspace{0.17em}\mathrm{s}\\ 11,& 500<t\le 600\hspace{0.17em}\mathrm{s}\end{array}$$

$$h_2^{(2)}(t)=10+\log (t^{0.5})\sin (\frac{t}{6})+t^{0.45}$$

and $h_2^{(3)}(t)$ which is an arbitrary waveform generated in MATLAB used here to model the variation of the heat transfer coefficient with time.

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. Data used for the heat conduction problem (the body is made of stainless steel (type 304) [20]).

$\dot{q}(\frac{\mathrm{W}}{{\mathrm{m}}^2})$	$k_T(\frac{\mathrm{W}}{\mathrm{m}·°\mathrm{C}})$	$\rho (\frac{\mathrm{kg}}{{\mathrm{m}}^3})$	$c(\frac{\mathrm{J}}{\mathrm{kg}·°\mathrm{C}})$	$h_2(\frac{\mathrm{W}}{{\mathrm{m}}^2·°\mathrm{C}})$	$\begin{array}{l}{h}_i(\frac{\mathrm{W}}{{\mathrm{m}}^2·°\mathrm{C}}),\\ i=3,4\end{array}$	$\begin{array}{l}{T}_{{\infty }_i}(°\mathrm{C}),\\ i=2,3,4\end{array}$
$2000+1000\sin (\frac{\pi t}{180})$	$14.9$	$7900$	$477$	$h_2^{(1)}(t),h_2^{(2)}(t),h_2^{(3)}(t)$	$5$	$30$

.

In all simulations in this study, the heat conducting body is meshed using a grid size of $M\times N=40\times 40$, the temperature measurement sensor is placed at node $(Si,Sj)=(\frac{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 }^{\circ }\mathrm{C}$, the final time is $t_r=600\hspace{0.17em}\mathrm{s}$, and the time step is $\Delta t=0.1\hspace{0.17em}\mathrm{s}$. Thus, the number of transient readings of the single sensor $S$ is $r=\frac{t_r}{\Delta t}=\frac{600\hspace{0.17em}\mathrm{s}}{0.1\hspace{0.17em}\mathrm{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 $\sigma =0.001$ and $\sigma =0.005$ are considered. The stopping criteria for the test case with the following measurement errors are

$$\sigma =0.001\Rightarrow \epsilon ={\sigma }^{2}r={0.001}^2(6000)=0.006$$

$$\sigma =0.005\Rightarrow \epsilon ={\sigma }^{2}r={0.005}^2(6000)=0.15$$

As the size of Jacobian matrix is $r\times r$, we will deal in this study with a Jacobian matrix of size $6000\times 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)=\frac{\partial T_e(Si,Sj,t_i)}{\partial k_T}$, $i=1,\dots ,r$ and $h_2{}_{{}_{k_T}}(j,1)=\frac{\partial h_2(t_j)}{\partial k_T}$, $j=1,\dots ,r$ are computed from Equations (21) and (23), respectively, and then the sensitivity coefficients can be obtained using the following pseudo-code

$$\begin{array}{l}\mathrm{do}\mathrm{i}=1,\mathrm{r}\\ \mathrm{do}\mathrm{j}=1,\mathrm{r}\\ \hspace{1em}\mathrm{if}(\mathrm{i}.\mathrm{LT}.\mathrm{j})\mathrm{then}\\ \hspace{1em}\mathrm{Ja}(\mathrm{i},\mathrm{j})=0.0\\ \hspace{1em}\mathrm{else}\\ \hspace{1em}\mathrm{Ja}(\mathrm{i},\mathrm{j})=\frac{{\mathrm{T}}_{{\mathrm{k}}_{\mathrm{T}}}(\mathrm{i},1)}{{\mathrm{h}}_2{}_{{}_{{\mathrm{k}}_{\mathrm{T}}}}(\mathrm{j},1)}\\ \hspace{1em}\mathrm{endif}\\ \mathrm{enddo}\\ \mathrm{enddo}\end{array}$$

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, $S(\frac{M}{2},N-3)$, 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.

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 Figure 8a, Figure 9a and Figure 10a (for the case of no measurement error, $\sigma =0$), Figure 8c, Figure 9c and Figure 10c (for the measurement error of $\sigma =0.001$), and Figure 8e, Figure 9e and Figure 10e (for the measurement error of $\sigma =0.005$). By comparing the desired and final functional forms shown in the above figures, it can be seen that the desired functional forms are recovered accurately which implies that the inverse analysis is not strongly affected by the errors involved in the temperature measurements due to the accuracy of the proposed sensitivity analysis scheme. When measurement errors exist, some oscillatory behaviors are observed around the exact values due to the ill-posed nature of the inverse heat transfer problem. The convergence histories of the objective function for the three functional forms of interest are shown in Figure 8b, Figure 9b and Figure 10b (for the case of no measurement error, $\sigma =0$), Figure 8d, Figure 9d and Figure 10d (for the measurement error of $\sigma =0.001$), and Figure 8f, Figure 9f and Figure 10f (for the measurement error of $\sigma =0.005$). The details of the results, including the initial and desired values for the unknown time-dependent heat transfer coefficient, the initial and final values of the objective function, and the number of iterations required to reach the solutions are given in

Table 2. A summary of results for the estimation of the time-dependent heat transfer coefficient.

Functional Form	Initial Guess, ${\mathit{h}}_{2_{\mathbf{initial}}}$	Temperature Measurement Error, $\mathit{\sigma }$	Initial Value of Objective Function	Final Value of Objective Function	Number of Iterations
$h_2^{(1)}$	$10$	$\sigma =0.0$	$3.50$	$9.99\times {10}^{-7}$	$3654$ (~244 min)
$h_2^{(1)}$	$10$	$\sigma =0.001$	$3.51$	$0.006$	$490$ (~33 min)
$h_2^{(1)}$	$10$	$\sigma =0.005$	$3.65$	$0.15$	$217$ (~14 min)
$h_2^{(2)}$	$25$	$\sigma =0.0$	$43.14$	$7.90\times {10}^{-7}$	$1037$ (~69 min)
$h_2^{(2)}$	$25$	$\sigma =0.001$	$43.14$	$0.006$	$487$ (~32 min)
$h_2^{(2)}$	$25$	$\sigma =0.005$	$43.26$	$0.15$	$371$ (~25 min)
$h_2^{(3)}$	$12$	$\sigma =0.0$	$34.28$	$6.95\times {10}^{-7}$	$948$ (~63 min)
$h_2^{(3)}$	$12$	$\sigma =0.001$	$34.28$	$0.006$	$387$ (~26 min)
$h_2^{(3)}$	$12$	$\sigma =0.005$	$34.42$	$0.15$	$197$ (~13 min)

. The computation time for each iteration (the direct and inverse solutions) is about 4 s. In spite of large unknown variables (6000 in the test case) and large final time, $t_r=600\hspace{0.17em}\mathrm{s}$, this short computation time confirms that the employed inverse analysis based on the proposed sensitivity analysis is very efficient. The results are obtained by a FORTRAN compiler and computations are run on a PC with Intel Core i5 and 6G RAM.

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 $\mathbf{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

$$\left[\begin{array}{c}0\\ 0\\ 0\\ \begin{array}{c}\begin{array}{c}⋮\\ 0\end{array}\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_r)}\end{array}\end{array}\right]$$

(33)

From Equation (19), ($\frac{\partial \mathrm{J}}{\partial {\mathbf{h}}_2}=2\mathbf{J}{\mathbf{a}}^{\mathrm{T}}\left[{\mathbf{T}}_e-{\mathbf{T}}_m\right]$), we can write the gradient of the objective function $\mathrm{J}$ with respect to ${\mathbf{h}}_2$ at the final time $t_r$ as

$$\frac{\partial \mathrm{J}}{\partial h_2(t_r)}=2{\left[\begin{array}{c}0\\ 0\\ 0\\ \begin{array}{c}\begin{array}{c}⋮\\ 0\end{array}\\ \frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_r)}\end{array}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}{T}_e(Si,Sj,t_1)-T_m(Si,Sj,t_1)\\ T_e(Si,Sj,t_2)-T_m(Si,Sj,t_2)\\ T_e(Si,Sj,t_3)-T_m(Si,Sj,t_3)\\ \begin{array}{c}\begin{array}{c}⋮\\ T_e(Si,Sj,t_{r-1})-T_m(Si,Sj,t_{r-1})\end{array}\\ T_e(Si,Sj,t_r)-T_m(Si,Sj,t_r)\end{array}\end{array}\right]=2\frac{\partial T_e(Si,Sj,t_r)}{\partial h_2(t_r)}\left[T_e(Si,Sj,t_r)-T_m(Si,Sj,t_r)\right]$$

(34)

which is a very small number. Substituting a very small value for $\nabla {\mathrm{J}}_{h_2(t_r)}^{(\mathrm{k})}$ into Equation (28), $d_{h_2(t_r)}^{(\mathrm{k})}=-\nabla {\mathrm{J}}_{h_2(t_r)}^{(\mathrm{k})}$, results in a very small value for $d_{h_2(t_r)}^{(\mathrm{k})}$. Likewise, substituting a very small number for $d_{h_2(t_r)}^{(\mathrm{k})}$ into Equation (27), $h_2{(t_r)}^{(\mathrm{k}+1)}=h_2{(t_r)}^{\left(\mathrm{k}\right)}+{\beta }^{\left(\mathrm{k}\right)}{d}_{h_2(t_r)}^{(\mathrm{k})}$, results in $h_2{(t_r)}^{(\mathrm{k}+1)}\approx h_2{(t_r)}^{\left(\mathrm{k}\right)}$, as observed. In other words, by approaching the final time $t_r$, 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.

5. 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 temperature-dependent 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.