1. 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 space-dependent 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 time-dependent 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.
2. Governing Equation
The heat-conducting body shown in
Figure 1a is made of a material of constant thermal conductivity
, density
, and specific heat
, and is initially at the temperature
. For the time
, a timewise varying heat flux
is applied at the boundary surface
. Convective heat transfer is imposed on boundary surfaces
with corresponding heat transfer coefficients
and surrounding temperatures
.
For this problem, the two-dimensional transient heat conduction equation with no heat generation is
The boundary and initial conditions are
where
is the time. Since the heat conduction problem is concerned with an irregular body, the irregular shape (the
and
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 (
) to the (
) variables [
2,
13]. The transformation gives
where
and
are grid control functions. If
, a smooth grid over the physical domain is generated. Thus, Equation (5) becomes
where
are the coefficients obtained from the elliptic grid generation method. The transformed boundary and initial conditions becomes
where
is the initial condition
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
where
. 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
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
,
, as follows
4. 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 times (). 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:
- (2)
A sinusoidal function in Test case 2:
- (3)
A triangular function in Test case 3:
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
Table 1.
The grid size is
, the sensor is placed at the node
(
Figure 2), the initial temperature is
, the final time is
, and the time step is
. Thus, the number of transient readings of the single sensor
is
. This implies that the number of unknown parameters is
. 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.
In this study, three different measurement errors of
are considered. The stopping criteria (Equation (31)) for the test cases with these measurement errors are
Since the size of the Jacobian matrix is
, we will deal here with a Jacobian matrix of size
. As soon as the nodal temperature at the sensor place is calculated at each time
, 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
,
and
,
are computed from Equations (20) and (22), respectively, and then the sensitivity coefficients can be obtained from
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 () and time is equal to , then the computational cost of the entire transient heat conduction solution for nodes at the final time is approximately equal to . On the other hand, the computational cost of Equation (23) for a single sensor placed at only one node, , at the final time is approximately equal to . Therefore, the computational cost of the computation of Jacobian matrix at each iteration is 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 (
) 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
, and the time step
are shown in
Figure 3a,b, respectively, and the temperature distribution obtained by the finite element software COMSOL is depicted in
Figure 3c. The temperature history of the place of the sensor,
, obtained by the explicit solver and the software COMSOL is depicted in
Figure 4. Moreover, the temperature distribution at the final time
at a given line, here along the nodes
,
,
, obtained by the explicit solver and the software COMSOL is shown in
Figure 5. The comparison between the results shows very good agreement, thereby confirming the correct implementation of the explicit solver.
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
Figure 6a,
Figure 7a,
Figure 8a (for the case of no measurement error,
),
Figure 6c,
Figure 7c,
Figure 8c (for the measurement error of
),
Figure 6e,
Figure 7e,
Figure 8e (for the measurement error of
), and
Figure 6g,
Figure 7g,
Figure 8g (for the measurement error of
). 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.
The convergence histories of the objective function for Test cases 1–3 and cases of no measurement error and the measurement errors of
,
, and
are shown in
Figure 6b,
Figure 7b,
Figure 8b (for the case of no measurement error,
),
Figure 6d,
Figure 7d,
Figure 8d (for the measurement error of
),
Figure 6f,
Figure 7f,
Figure 8f (for the measurement error of
), and
Figure 6h,
Figure 7h,
Figure 8h (for the measurement error of
), 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
in Test case 2 in which the number of iterations is 467 (larger than 401 for
). The reason is that the iterative process for the case of the measurement error of
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 ill-posed 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
Table 2 (for cases of no measurement error and the measurement errors of
,
, and
). As can be seen in
Table 2, in spite of large unknown variables (10000 in these test cases) and large final time, (
), the short computation time in the test cases 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.
It can also be seen that in a neighborhood of
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
as
By approaching the final time
, 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,
, is
which is a very small number. Substituting a very small value for
into Equation (27),
, gives a very small value for
. Likewise, substituting a very small number for
into Equation (26),
, results in
, as observed.