Analysis of a Fundamental Procedure for Solving the Inverse Heat Conduction Problem Applied to Simple Systems

Abstract

Heat transfer phenomena develop in various natural and artificially created processes. Fundamental laws of physics allow the transfer mechanisms to be classified; however, describing the phenomena is relatively complex, even if the analysis is limited to conduction. In particular, to determine the temperature distribution in a solid body, the definition of the boundary conditions that perturb it is required. Such conditions mathematically model a fluid’s hydrodynamic and thermodynamic behavior, and when the temperature differences are significantly high, the flow by radiation. It is then complex to define the functions of the thermal boundary conditions and solve the well-posed problem. Naturally, nonlinear system results and applying numerical methods are constant in the analysis. However, a unique solution for the thermal field in a solid does ensure. Alternatively, the scheme of the discretization of the system allows us to propose that through the knowledge of a fraction of the thermal field, the boundary condition is quantified independently of its nature. Such a procedure is called inverse analysis and has the characteristic of not satisfying the single solution criterion. However, some cases of interest can treat, and the estimate is guaranteed to be highly accurate.

1. Introduction

The researchers Dowding and Beck [1] define that inverse heat conduction problems (IHCPs) can be categorized in various ways, (1) by the solution technique or algorithm (Function Specification, Regularization, Laplace Transform, Conjugate Gradient, Mollification, etc.), (2) by the solution method (Duhamel’s Theorem, Finite Difference Method, Finite Element Method, etc.), and (3) by the time domain (Stoltz Method, Sequential, and Complete Domain Method). Beck et al. [2] developed a whole theory, and other researchers have used it to adapt it to particular objectives [3, 4, 5]. Therefore, the characteristic of the inverse problems of not meeting the criteria of existence and uniqueness of the solution has promoted the implementation of more robust mathematical methods in each classification of the IHCP to optimize its stability since minor errors in the measured data can induce significant inaccuracies in the estimated variable [6, 7, 8, 9, 10, 11, 12]. The analyses have oriented the study with an analytical mathematical approach [13, 14, 15] and to applying numerical solutions in various physical phenomena, particularly thermal phenomena [16, 17, 18, 19].

2. Inverse heat conduction problem

Transient thermal phenomena are relatively complex to be described through mathematical formulations. However, it is possible to obtain acceptable approximations in systems that can be bounded. In this regard, the problems posed by Beck et al. [2] for a solid with cylindrical geometry will be formulated and developed. With this objective, the numerical method based on heat flow balances results in sufficient precision and easy programming.

It begins with defining the heat conduction equation in its differential form for the cylindrical coordinate system represented by the variables r, θ, and z. Eq. (1) specifies the material’s properties through k, ρ, and Cp; thermal conductivity; density; and thermal capacity, respectively. The field or dependent variable T is temperature, and t represents the variable time. Finally, $\dot{q}$ defines the heat generated or consumed per unit volume.

$$-\left[\frac{1}{r}\frac{\partial }{\partial r}\left(-\mathit{kr}\frac{\partial T}{\partial r}\right)+\frac{1}{r}\frac{\partial }{\partial \theta }\left(-k\frac{1}{r}\frac{\partial T}{\partial \theta }\right)+\frac{\partial }{\partial z}\left(-k\frac{\partial T}{\partial z}\right)\right]+\dot{q}=\frac{\partial \rho C_{p}T}{\partial t}$$E1

Mainly, to clarify the process of estimating thermal boundary conditions through inverse analysis, interference from internal heat “sources” or “wells” is not desired, so the $\dot{q}$ term did cancel. From this new form of Eq. (1), the formulations of the specific problems will be developed, which will address the inverse analysis.

2.1 One-dimensional analysis in the radial direction of heat transfer—An ill-posed problem for estimating the function in the active boundary

Considering the dimensional restriction to a flow direction simplifies the problem. Additionally, although the physical properties associated with a material depend on temperature, due to the sequential numerical solution strategy that will be adopted, constants can be defined at the instant of calculation time. That is, it is evident that the temperature distribution is unknown at the moment of calculation, but is know at the last instant. Therefore, if it does consider that there will be no significant variations in the properties due to relatively small temperature changes, then they are estimated and projected to the calculation time step. In this way, the well-posed or direct heat conduction problem (DHCP) is formulated according to the information in Figure 1.

Figure 1.

The mathematical formulation of the DHCP in the radial direction and transitory state.

The problem is mathematically nonlinear; therefore, numerical analysis is required to obtain a solution for the temperature variable. The finite difference method was selected based on the systems’ geometric characteristics and the developments’ intelligibility.

In the methodology, the system domain is subdivided into adjacent sections creating so-called control volumes (CVs). Figure 2 shows a simplified scheme to indicate the representative basic rings according to the definition of the geometry and direction of heat flow. As can be deduced, the notation ∆r represents a fraction of the radius of the cylinder.

Figure 2.

Discretization of the cylindrical system considering a one-dimensional analysis in the radial direction.

Equation (2), called the general equation of heat flux balance Q, is equivalent to the mathematical model presented in Figure 1 and will be applied to each CV in the scheme.

$$Q_{\mathit{\text{In}}}-Q_{\mathit{\text{Out}}}=Q_{\mathit{\text{Accumulation}}}$$E2

By definition of Fourier’s law for the heat flow entering Q_In and leaving Q_Out, in addition to the heat flow per unit volume that accumulates in a transitory state Q_Accumulation, the system of equations corresponds to the number of CV of meshed space. If a Crank–Nicolson technique approach does consider calculating the field variable concerning the transient state, the expressions in Table 1 are obtained.

VC	General heat balance equation	$i$
1	$q_i^t{A}_{\mathit{\text{In}},1}-\left(-\overline{k_{i\to i+1}}{A}_{\mathit{\text{Out}},i}\frac{T_{i+1}^t-T_i^t}{∆\mathrm{r}}\right)={\rho }_i{C}_{\mathit{pi}}{V}_1\frac{T_i^{t+∆t}-T_i^t}{∆t}$	$1$
2	$-\overline{k_{i-1\to i}}{A}_{\mathit{\text{In}},i}\frac{T_i^t-T_{i-1}^t}{∆\mathrm{r}}-\left(-\overline{k_{i\to i+1}}{A}_{\mathit{\text{Out}},i}\frac{T_{i+1}^t-T_i^t}{∆\mathrm{r}}\right)={\rho }_i{C}_{\mathit{pi}}{V}_i\frac{T_i^{t+∆t}-T_i^t}{∆t}$	$2,\dots ,n-1$
3	$-\overline{k_{n-1\to n}}{A}_{\mathit{\text{In}},i}\frac{T_n^t-T_{n-1}^t}{∆\mathrm{r}}-0={\rho }_n{C}_{\mathit{pn}}{V}_n\frac{T_n^{t+∆t}-T_n^t}{∆t}$	$n$

Table 1.

Summary of the development of the general equation of balance of heat flow Eq. (2) on the CVs in the discretization of the cylindrical system in the radial direction.

In each system considered, a mesh sensitivity analysis should do perform. It can deduce that obtaining results with greater precision of the temperature distribution requires defining a more significant number of CVs. Therefore, it is required to determine the minimum number of CVs in which no differences are detected in the calculation result. On the other hand, the parameter $∆\mathrm{r}$ can be heterogeneous when system sections are defined with a different degree of mesh. The strategy helps to focus sensitivity on important sections and to make computational time efficient.

The function $\mathbf{F}$ specified in Figure 1 can be defined as Dirichlet, Neumann, or Robinson type. However, the general heat flux balance equation determines a Neumann-type function. Therefore, the general form of the equations can be stated as follows:

$$d_i{T}_{i-1}^t+e_i{T}_i^t+f_i{T}_{i+1}^t=T_i^{t+∆t}-g_i,i=1,2,\dots ,n$$E3

Where d, e, f, and g represent the respective terms product of the algebraic development of the corresponding equations in Table 1. In this way, the system of equations of the type does generate:

$$\left[\begin{array}{ccccc}{e}_1& f_1& \cdots & 0& 0\\ d_2& e_2& f_2& 0& 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& 0& d_{n-1}& e_{n-1}& f_{n-1}\\ 0& 0& \cdots & d_n& e_n\end{array}\right]\left[\begin{array}{c}{T}_1^t\\ T_2^t\\ ⋮\\ T_{n-1}^t\\ T_n^t\end{array}\right]=\left[\begin{array}{c}{T}_1^{t+∆t}-g_1\\ T_2^{t+∆t}-g_2\\ ⋮\\ T_{n-1}^{t+∆t}-g_{n-1}\\ T_n^{t+∆t}-g_n\end{array}\right]$$E4

Implementing solution strategies for matrix systems will depend on their structure, but the Thomas method, successive over-relaxation, and LU decomposition do suggest. A solution to the well-posed or direct problem is relevant because it is required to solve the ill-posed or inverse problem. The algorithm for solving the direct heat conduction problem establishes the core of the solution of the inverse analysis since it is used repeatedly in each cycle of estimating the thermal boundary conditions depending on the value of the stabilization parameter.

In turn, the wrongly stated or inverse problem defines that the boundary condition at r = R is unknown; however, if there is knowledge of the thermal history obtained during a monitoring event at least one point within the solid. The formulation is summarized in Figure 3.

Figure 3.

Formulation of the ill-posed problem or IHCP for the same system referenced in Figure 1.

The implementation of Beck’s theory, based on the sequential method of specification of the function, refers to the determination in each calculation step of heat flux density q(t). Estimating the value of data of the function depends exclusively on the internal gradients in the solid and not on the potential extraction or supply of heat by an external medium (as foreseen, a Neumann-type boundary condition does define as an unknown, the reason over the Dirichlet and Robinson types is the natural condition of their specification through the heat flux general balance equation).

It is more common and with greater dominance that in the field of applicability of boundary conditions in heat transfer problems, Robinson-type conditions identify where a heat transfer coefficient “h” explicitly expresses. The magnitude of the values of a function h represents complex behavior conditions of media whose primary heat transfer mechanism is convection. However, some data may be associated with a combination with a radiation mechanism. In addition, they are related to the temperature of the surface of the solid and to the temperature of the fluid sine through Newton’s law of cooling, defined by Eq. (5).

$$q=h\left(T_s-T_{\infty }\right)$$E5

Particularly, it is possible to determine a function of the Robinson-type boundary condition with the solution methodology of the IHCP. However, Newton’s law of cooling approach determines that once the difference between the surface and bulk temperatures of the fluid tends to zero, the function does not adequately satisfy the thermal boundary condition that causes small gradients in the solid. This condition exposes a field of research of great interest. The solution by inverse analysis could support the modification of Newton’s cooling law, being used as a reference to model a mathematical complement that allows determining the appropriate value of the heat transfer coefficient at surface temperatures close to the temperature of the fluid bulk.

Returning to the specification of a Neumann-type condition is justified to describe the complete transitory heat transfer process. For this purpose, applying the norm of least squares defined by Eq. (6) as an objective function is proposed.

$$S=\sum _{l=1}^{\ell }\sum _{j=1}^n{\left(Y_{j,M+l-1}-T_{j,M+l-1}\right)}^2$$E6

The objective function establishes that the square of the difference between the registered and calculated temperatures, proposing a value of the border condition, which will be minimized. Since more recorded histories can integrate, then the minimization of the sum of the squares of all the differences is considered. In addition, since there is the possibility of knowing the temperature data after a reference instant, as much data as necessary do use to obtain weights of the immediate “future” and is also considered in the minimization. This last definition does use to stabilize the determination of the magnitude of the value of the boundary condition. As it is associated with the sequence of calculation time steps, the amount of data considered is called the number of future time steps (ℓ). Naturally, to achieve stabilization, the minimum value of the number must be two, and for each calculation, the matching value of the proposed boundary condition does consider (q_M-1 = 0 is recommended). Following the above, the mathematical development for determining the boundary condition’s value does obtain by the expression of Eq. (7).

$$\begin{array}{l}{q}_M-q_{M-1}=\frac{1}{∆M}\sum _{l=1}^{\ell }\sum _{j=1}^m\left(Y_{j,M+l-1}-T_{j,M+l-1}^{\ast }\right)\left(T_{j,M+l-1;q}^{\ast }\right),\\ ∆M=\sum _{l=1}^l\sum _{j=1}^m{\left(T_{j,M+l-1;q}^{\ast }\right)}^2,\\ T_{j,M+l-1;q}^{\ast }=\frac{\partial T_{j,M+l-1}}{\partial q_M}\end{array}$$E7

As can be deduced, considering ℓ = 2 or greater implies solving the same number of times the formulation of the DHCP. In addition, as seen in Eq. (7), according to the term called sensitivity coefficient, it is necessary to solve again the same number of times the result of the derivation of the same formulation. Figure 4 shows the compendium of the corresponding development.

Figure 4.

Mathematical formulation for the calculation of the sensitivity coefficients for the solution of the IHCP of the cylindrical system with radial unidirectional flow and a single boundary condition.

Since a large amount of information must be processed to specify a single piece of data for the function q(t), it is convenient to focus efforts on developing a computer application. Figure 5 contains the flowchart of the pre-explanation calculation algorithm to obtain the sequential specification of the complete transient event function.

$$S_Y={\left[\frac{N}{N-1}\sum _{\mathrm{\aleph }=1}^N{\left(Y_{\mathrm{\aleph }}\left(t\right)-{\widehat{Y}}_{\mathrm{\aleph }}\left(t\right)\right)}^2\right]}^{\frac{1}{2}}$$E8

Figure 5.

Flowchart of the algorithm for solving the IHCP.

Currently, there is the advantage of having efficient computer systems. In this sense, the algorithm can modify to include decision-making parameters to optimize the solution’s stability without influencing the estimation’s precision. For example, a weighted error parameter at each calculation step is defined by Eq. (8). The equation represents the calculation of the standard deviation by relating the recorded and calculated thermal histories to the function of the estimated boundary condition.

The degree of error will be associated with the “quality” of the information of each thermal history fed to the calculation algorithm and the designation of the number of future time steps (l). A high value of l helps to stabilize the IHCP solution when the quality of the thermal history is low; that is, the imprecision due to the effect of noise introduced in the data does minimize. However, when the quality of the thermal history is high, the magnitude of the function q(t) values is underestimated. The underestimation can be noted graphically by the trend of a smooth transition of the function q(t) curve, above all, in moments where the temperature could change suddenly. In this way, with a previous analysis of the quality of the thermal history, the error parameter can be used to include an adjustment process of l that does not imply a critical underestimation of the function q(t) concerning the thermal information used.

Fundamentally, the previous implies that before analyzing natural systems, exercises are carried out to verify the formulations and the coding of algorithms through submission to benchmark cases. The fundamental methodology fulfills at least the following steps:

• Establish a system with a known thermal boundary condition and calculate thermal histories solving the well-posed or direct heat conduction problem. Based on a thermodynamic analysis, negative data indicates heat extraction from the solid.

• Choose at least one thermal history calculated at a specific position, preferably the location close to the border of interest.

• Provide at least one thermal history regarding the point before the application that solves the wrongly stated or inverse problem.

• Compare the estimated thermal boundary condition function, as well as the thermal history calculated with it, with the corresponding data from steps 1 and 2. A reliable IHCP solution tool obtains if there are no significant differences. Otherwise, one must try to stabilize the solution through the number of future time steps or modify the mathematical formulation.

In Figure 6, the graphs obtained from the previous steps are summarized. It can be noted that abrupt magnitude changes in the boundary condition represent extreme test points for the IHCP solution method. The bias presented by the estimated function establishes a research topic that should be dealt with separately from the discussion here. However, the method can project to natural systems with the results obtained.

Figure 6.

Benchmark boundary conditions q(t) with (a) triangular and (b) pulse shapes, used to assess the accuracy of the estimate $\widehat{q}\left(t\right)$ by solving the IHCP.

The sensors and data acquisition systems must calibrate before a metrology certifying unit. On the other hand, the systems coupled to the measurement target must be shielded as well as possible to avoid environmental disturbances and human error. In essence, for diverse systems, the quality of the measurements of the field variable used to estimate boundary conditions represents the paradigm of inverse analysis. In addition, there is the component of access to the system of interest, often made impossible by the operating conditions, the dimensions, or the specific location. However, with the current advancement of filters in electronic and computer systems, it has not been necessary to increase the complexity of mathematical methods to stabilize the solution of ill-posed problems.

An example of a direct, unfiltered measurement of thermal history on a cylindrical specimen is presented in Figure 7. The trend from the beginning (t = 0 s) to end (t = 34 s) in the measurement represents a decrease in temperature. Therefore, it is defined as a cooling curve. The measurement was obtained in a cylinder that meets unidirectional flow conditions by keeping an equivalent Height-Diameter ratio of 4:1. The AISI 304-type stainless steel material guarantees that there will be no exothermic or endothermic internal reactions per unit volume equivalent to the $\dot{q}$ term. Installing a type K bayonet thermocouple for temperature sensing required drilling with a 1 mm diameter drill bit 1.6 mm from the curved surface from the top to half the height of the specimen. The diameter of the analyzed specimen is 12.7 mm.

Figure 7.

Thermal history recorded inside a cylindrical specimen where heat flow is obtained by cooling from a high temperature. In this case, the thermal history $Y\left(t\right)$ is called the “cooling curve”.

The results of estimating the boundary condition using the thermal history and solving the IHCP are expressed as the function $\widehat{q}\left(t\right)$ in the graph of Figure 8a. In turn, with the history heat flux function, the temperature field in the cross-section of the cylinder was determined (solving the DHCP), and the thermal history of the CV where the thermocouple does virtually locate were identified. Figure 8b graph plots the measured Y(t) (empty circular marker) and calculated $\widehat{Y}\left(t\right)$ thermal histories (continuous line). As can be seen, no significant differences are observed between the thermal histories, indicating that the boundary condition has done the estimate with good precision.

Figure 8.

Results of the application of the IHCP solution algorithm for the determination of (a) the function of the heat flux history $\widehat{q}\left(t\right)$ and (b) the thermal history in the virtual position of the thermocouple $\widehat{Y}\left(t\right)$ (solid line).

However, when analyzing the standard deviation weighting, an optimization modification did implement in the IHCP solution algorithm, specifically, the option to auto-adjust the number of future time steps between 2 and a user-defined upper limit, comparing the results to the standard deviation. In this process, the value of l is increased by one unit in each loop cycle. If the magnitude of the deviation is less than the previous one, then a new increment is performed. Otherwise, the result does keep with the last best-estimated value. Figure 9 shows the standard deviation calculation called the residual at each calculation step through the dashed curve. When implementing the recursive process for determining ℓ considering a maximum value of 10, it can be seen through the dotted line that the residual decreased. The most notable difference does find between 23 and 30 s. The history of ℓ is also recorded and does show by the solid line with filled square symbols.

Figure 9.

Calculate the (residual) standard deviation used to determine the value of the number of future time steps ℓ (solid line with filled square marker). The curves (a) dashed line and (b) dotted line show the differences in the magnitude of the residual by using a constant and adjustable value of ℓ, respectively.

The result of the IHCP solution with ℓ adjustable and ℓ = 2 does compare in the curves of Figure 10. As seen in the heat flux histories curves of Figure 10a, a value of ℓ greater than 2 tends to underestimate. However, it also indicates a tendency toward stabilizing the solution if an alteration occurs due to disturbances in the experimental thermal history. In this case, the comparison of the deviation percentages of the thermal accounts calculated in the virtual position of the thermocouple in Figure 10b shows that the increase of ℓ improved the estimation of the function.

Figure 10.

Differences in (a) the estimation of the heat flux function with a constant value of ℓ = 2 (continuous line) and adjustable to a maximum of 10 (open circles markers), evidenced by the comparison of (b) the percentage of deviation between the experimental temperature history and that calculated with ℓ = 2 (solid line) and with adjustable ℓ (dotted line).

According to the above, training in understanding the phenomena of heat transfer by the conduction mechanism and the numerical analysis allows us to identify the integrity and precision of the estimation of boundary conditions by solving the problems posed in the reverse form. In addition to the materials’ constitutive properties, the analysis’s complexity increases with the number of active boundaries and dimensional coordinates, as seen in the studies in the following sections.

2.2 One-dimensional analysis of heat flow in radial direction: Ill-posed problem for estimating one function in each frontier

Modifying the system of the previous section by making a circular cut centered on the axial axis generates a new thermal boundary without insulation. Therefore, it is necessary to formulate the new IHCP considering a dual estimation procedure. Figure 11 outlines the case study with the following characteristics:

• The concentric conditions of the surfaces in the radial coordinate determine that the gradients in the angular direction of the system are not significant. It is also established that the height-diameter ratio meets the minimum requirement of 4:1.

• The heat conduction in the solid does determine by ignoring the discontinuity in the material due to the insertion of temperature sensors (instrumentation with thermocouples), and the absence of sources or internal heat sinks in the solid does consider.

• The material is considered isotropic, and the thermophysical properties (k, ρ, and Cp) are a function of temperature.

• A homogeneous heat flux history q1(t) can specify the boundary condition on the inner surface.

• A homogeneous heat flux history q2(t) can specify the boundary condition on the external surface.

• The initial temperature field is known.

• At least two points, information is available on the temperature distribution inside the body to adequately estimate the unknown thermal boundary conditions.

Figure 11.

Formulation of the IHCP for a cylindrical system with the unidirectional flow in the radial direction with two active boundaries.

The solution of the mathematical formulation is obtained again by applying numerical methods mainly due to:

1. The nonlinearity of the partial differential equation in the transient state and thermophysical properties dependent on the temperature that governs the heat conduction problem,

2. The nonlinearity of the border conditions, as well as,

3. The vast magnitude of the time domain (depending on the density of experimental information) in which the phenomenon does analyze.

In this way, the sequential technique of estimation of the function at intervals of analysis delimited in time (fractions of the total time of the observed phenomenon) allows transforming the nonlinear problem into a “temporarily” linear one with a solution of the system for each step of time through the approach of an “objective function” and the introduction of the concept of “number of future time steps” used for the estimation calculation.

According to Beck’s approach [2], (Eq. (9)), the result of a Taylor series expansion, determines the matrix operation for calculating the temperature distribution concerning the operational boundary conditions.

$$\mathbf{T}={\mathbf{T}}_{\mathbf{q}=0}+\mathbf{Xq}$$E9

Carrying out a slight adaptation of the notation to include a positioning matrix, Eq. (10) does obtain.

$$\mathbf{q}=\mathbf{A\beta },\mathbf{A}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],\mathbf{\beta }=\left[\begin{array}{c}{q}_1\left(t\right)\\ q_2\left(t\right)\end{array}\right]$$E10

The Eq. (9) is rewritten to get Eq. (11):

$$\mathbf{T}={\mathbf{T}}_{\mathbf{\beta }=0}+\mathbf{Z\beta },\mathit{Z}=\mathit{XA}$$E11

Therefore, the target equation was modified as follows to obtain Eq. (12):

$$S={\left(\mathbf{Y}-\mathbf{T}\right)}^T\left(\mathbf{Y}-\mathbf{T}\right)$$E12

Substituting Eq. (11) into Eq. (12) results in the Eq. (13):

$$S={\left(\mathbf{Y}-{\mathbf{T}}_{\mathbf{\beta }=0}+\mathbf{Z\beta }\right)}^T\left(\mathbf{Y}-{\mathbf{T}}_{\mathbf{\beta }=0}+\mathbf{Z\beta }\right)$$E13

Therefore, when performing the minimization operation of Eq. (13) concerning β, the relationship defined by Eq. (14) is deduced to estimate the unknown boundary conditions:

$$\widehat{\mathbf{\beta }}={\left({\mathbf{X}}^T\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^T\left(\mathbf{Y}-{\mathbf{T}}_{\mathbf{\beta }=0}\right)$$E14

The development of the algebraic operations indicated by Eq. (14) determines Eqs. (15)–(17) as the applicable relations for any set of thermal histories and the parameter value of the number of future time steps for stabilizing the solution.

$$\widehat{\mathbf{\beta }}=\left[\begin{array}{c}{\widehat{q}}_1\left(M\right)\\ {\widehat{q}}_2\left(M\right)\end{array}\right]=\frac{1}{\left(C_{11}{C}_{22}-C_{12}{C}_{12}\right)}\left[\begin{array}{c}{C}_{22}{D}_1-C_{12}{D}_2\\ -C_{12}{D}_1+C_{11}{D}_2\end{array}\right]$$E15

$$C_{\mathit{uv}}=\sum _{l=1}^{\ell }\sum _{j=1}^m\left[\sum _{p=1}^l{a}_{\mathit{ju}}\left(p\right)\right]\left[\sum _{p=1}^l{a}_{\mathit{jv}}\left(p\right)\right],u=1,2;v=1,2$$E16

$$D_w=\sum _{l=1}^{\ell }\sum _{j=1}^m\left[\prod _{p=1}^l{a}_{\mathit{jw}}\left(p\right)\right]\left[{\mathbf{Y}}_{j,M+p-1}-{\mathbf{T}}_{j,\mathbf{\beta }\left(M+p-1\right)=0}\right],w=1,2$$E17

Where the terms $a_{\mathit{ju}},a_{\mathit{jv}},a_{\mathit{jw}}$ are components of the matrices, which in turn are the components of the sensitivity coefficient matrix, which is defined by Eq. (18):

$$a_{\mathit{jK}}\left(p\right)=\frac{\partial T\left(\mathbf{x}j,t_p\right)}{\partial \left[q_K\left(1\right)\right]}$$E18

The field of sensitivity coefficients does calculate by the same procedure used for the example in the previous section. Therefore, it is necessary to formulate the equivalent problem expressed in Figure 4. From Eq. (18), it can be deduced that the term x_j refers to the histories of the sensitivity coefficients corresponding to the positions of the sensors used to record the thermal histories Y_j.

According to the above, Figure 12 summarizes the graphs of the analysis of heat extraction from a 50.8 mm external diameter cylinder with a 12.7 mm bore. The hollow cylinder was heated to a temperature of 150°C and exposed to room temperature airflow through the central duct and to still air on the outer surface. The holes to instrument the hollow cylinder near its thermally active borders met the same conditions as in the case presented in Section 2.1 (1 mm diameter at 1.6 mm distance from the surfaces).

Figure 12.

Comparison of recorded thermal histories $Y_j\left(t\right)$ and calculated ${\widehat{Y}}_j\left(t\right)$ at (a) inner boundary and (b) outer boundary of the hollow cylinder system, and estimated thermal boundary conditions defined as (c) heat flux $\widehat{q}\left(t\right)$ and (d) heat flow $Q\left(t\right)$, respectively. Figure prepared by the authors.

Figure 12a and b show the recorded cooling curves $Y_j\left(t\right)$ (open circle markers) and the calculated ${\widehat{Y}}_j\left(t\right)$ (solid lines) with the DHCP solution using the boundary conditions thermals estimated by the IHCP solution. Figure 12c and d show the functions $q_1\left(t\right)$ and $q_2\left(t\right)$ in terms of heat flux and heat flow to note that the magnitude of the surface may indicate a misunderstanding on which boundary is more heat extracted. Initially, the most significant amount of heat does remove through the inner surface, where a forced convective phenomenon occurs experimentally.

2.3 Bidimensional analysis of heat flow in radial and axial directions: Ill-posed problem for estimating multiple functions in one frontier

Another case study from which it can be deduced that applying the IHCP solution for estimating boundary conditions can extrapolate to more complex systems will present below. Again, the characteristics of the system are cited to define the governing heat transfer equation:

• Consider a disk of small thickness where the characteristics of symmetric conditions that determine the gradients in the angular direction of the system are not significant. Nevertheless, a heterogeneously distributed boundary condition does determine from its geometric center on one of the faces with the largest surface. In contrast, on the opposite face, it is considered homogeneous and known.

• The heat conduction in the solid does determine by ignoring the discontinuity in the material due to the insertion of temperature sensors (instrumentation with thermocouples), and the absence of sources or internal heat sinks in the solid does consider.

• The material is considered isotropic, and the thermophysical properties (k, ρ, and Cp) are a function of temperature.

• The heterogeneously distributed boundary condition can be specified by a series of discrete heat flux histories bounded by specific adjacent annular regions of defined width.

• A homogeneous heat flux history can specify the boundary condition on the opposite surface.

• The initial temperature field is known.

• At least one point each region information is available on the temperature distribution inside the body to adequately estimate the unknown thermal boundary conditions.

Figure 13 shows the schematization of the system of interest and the corresponding mathematical formulation.

Figure 13.

Formulation of the IHCP for a cylindrical system with the bidirectional flow in the radial and axial directions with two active boundaries.

In this analysis, the mathematical treatment of one of the two previous cases can be implemented to obtain the solution of the IHCP through the objective function. In this case, the methodology followed in Section 2.1 will be chosen. It is worth mentioning that in the case of Section 2.2, the distribution matrix of the domain of the functions on the boundaries of the system must formulate. Figure 14 and Eq. (19) determine the discretization of the disc in two dimensions and the type of penta-diagonal matrix of the system of equations, respectively. However, the system can increase according to the product (number of nodes in the r-direction) * (number of nodes in the z-direction). Again, matrix solution techniques are recommended by LU decomposition or successive over-relaxation (SOR).

Figure 14.

Simplified scheme of the discretization in two dimensions (r, z) to define the CVs in the disk system.

$$\left[\begin{array}{cccccccccccccc}{b}_{11}& c_{12}& 0& 0& 0& 0& e_{21}& 0& 0& 0& 0& 0& 0& \dots \\ a_{11}& b_{12}& c_{13}& 0& 0& 0& 0& e_{22}& 0& 0& 0& 0& 0& \dots \\ 0& a_{12}& b_{13}& c_{14}& 0& 0& 0& 0& e_{23}& 0& 0& 0& 0& \dots \\ 0& 0& a_{13}& b_{14}& c_{15}& 0& 0& 0& 0& e_{24}& 0& 0& 0& \dots \\ 0& 0& 0& a_{14}& b_{15}& c_{16}& 0& 0& 0& 0& e_{25}& 0& 0& \dots \\ 0& 0& 0& 0& a_{15}& b_{16}& 0& 0& 0& 0& 0& e_{26}& 0& \dots \\ d_{11}& 0& 0& 0& 0& 0& b_{21}& c_{22}& 0& 0& 0& 0& e_{31}& \dots \\ 0& \ddots & 0& 0& 0& 0& \ddots & \ddots & \ddots & 0& 0& 0& 0& \ddots \end{array}\right]$$

$$\left[\begin{array}{c}{T}_{11}^t\\ T_{12}^t\\ T_{13}^t\\ T_{14}^t\\ T_{15}^t\\ T_{16}^t\\ T_{21}^t\\ T_{22}^t\\ T_{23}^t\\ T_{24}^t\\ T_{25}^t\\ T_{26}^t\\ T_{31}^t\\ ⋮\end{array}\right]=\left[\begin{array}{c}{T}_{11}^{t+∆t}-f_{11}\\ T_{12}^{t+∆t}-f_{12}\\ T_{13}^{t+∆t}-f_{13}\\ T_{14}^{t+∆t}-f_{14}\\ T_{15}^{t+∆t}-f_{15}\\ T_{16}^{t+∆t}-f_{16}\\ T_{21}^{t+∆t}-f_{21}\\ T_{22}^{t+∆t}-f_{22}\\ T_{23}^{t+∆t}-f_{23}\\ T_{24}^{t+∆t}-f_{24}\\ T_{25}^{t+∆t}-f_{25}\\ T_{26}^{t+∆t}-f_{26}\\ T_{31}^{t+∆t}-f_{31}\\ ⋮\end{array}\right]$$E19

As previously defined, the objective function does minimize concerning the set of estimable heat flux density functions q, such that Eq. (20) is determined:

$$\frac{\partial S}{\partial \mathbf{q}}=\frac{\partial }{\partial \mathbf{q}}\left[\sum _{l=1}^{\ell }\sum _{j=1}^m{\left(Y_{j,M+l-1}-T_{j,M+l-1}\right)}^2\right]=0$$E20

The procedure leads to the expression of the Eq. (21):

$$0=\sum _{l,j}\left(Y_{j,M+l-1}-T_{j,M+l-1}\right)\frac{{\partial T}_{j,M+l-1}}{{\partial q}_{\beta ,M}},\beta =1,2,\dots ,5$$E21

A general discrete definition of the Taylor series expansion for the determination of the thermal field does express as (Eq. (22)):

$$T_M=T_{M-1}+f\left[T_{M-1},q_{M-1}\right]∆q$$E22

The differential form of Eq. (22) applied to the case study can be expressed as Eq. (23):

$$T_{j,M+l-1}=T_{j,M+l-1}^{\ast }+\sum _{l,j}\left(\frac{{\partial T}_{j,M+l-1}^{\ast }}{{\partial q}_{\beta ,M}}\right)\left(q_{\beta ,M}-q_{\beta ,M-1}\right),\beta =1,2,\dots ,5$$E23

By substituting Eq. (22) in Eq. (20), the expression for the solution of the IHCP obtain Eq. (24):

$$\sum _{\alpha =1}^n\left(q_{\beta ,M}-q_{\beta ,M-1}\right)\sum _{l,j}\left(\frac{{\partial T}_{j,M+l-1}^{\ast }}{{\partial q}_{\beta ,M}}\bullet \frac{{\partial T}_{j,M+l-1}^{\ast }}{{\partial q}_{\beta ,M}}\right)=$$

$$\sum _{l,j}\left(Y_{j,M+l-1}-T_{j,M+l-1}^{\ast }\right)\frac{{\partial T}_{j,M+l-1}^{\ast }}{{\partial q}_{\beta ,M}},\beta =1,2,\dots ,5$$E24

Eqs. (22) and (23) again show the definition of the sensitivity coefficient. The data set corresponding to the histories at the position of the j^th sensor due to a change in the heat flux in the β^th domain of the boundary at the computation time M defines the equivalent of the cooling rate histories recorded at those points of the solid. As can be deduced, each set is independent and isolated in the boundary condition over the domain defined in front of the sensor at the time of calculation without considering changes in the other boundary conditions.

Figure 15 shows the results of a benchmark case to complete the exposition. Tables 2–6 summarize the information on the direct and inverse problem statements in Figure 13. The material is considered AISI 304 stainless steel; therefore, the average values of the thermophysical properties of this material did use. The analyzed disk has dimensions of 100 mm in diameter and a thickness of 6 mm.

Figure 15.

Benchmark case for the analysis of the multiple estimations of heterogeneously distributed thermal boundary conditions over the larger-magnitude surfaces of a disk; (a) the proposed $q_i\left(t\right)$ functions do virtually distribute in ω_i sections, and DHCP is solved to obtain (b) the thermal histories $Y_i\left(t\right)$ at virtual sensors positions. In turn, these do use to solve the IHCP estimating the thermal boundary conditions ${\widehat{q}}_i\left(t\right)$, and (c) they are compared with the origin functions, and finally (d) the thermal histories ${\widehat{Y}}_i\left(t\right)$ originated by a new DHCP solution are compared in this cycle.

ω1	ω2	ω3	ω4	ω5
0 ≤ r < 0.005	0.005 ≤ r < 0.015	0.015 ≤ r < 0.025	0.025 ≤ r < 0.035	0.035 ≤ r ≤ 0.05

Table 2.

Dimensions of the bounded regions for the specification of the heterogeneous boundary condition (0 ≤ r < R).

Y1	Y2	Y3	Y4	Y5
(0, 0.0018)	(0.01, 0.0018)	(0.02, 0.0018)	(0.03, 0.0018)	(0.04, 0.0018)

Table 3.

Ordered pair (r, z) to define the location of the temperature sensors Y_j(t).

Material properties	k (Wm^−1°C^−1)	ρ∙Cp (Jm^−3°C^−1)
	22.0	5.125 × 10^6

Table 4.

Average thermo-physical properties of AISI 304 type stainless steel.

CVr	CVz	SOR parameter	ℓ
50	10	0.6	3

Table 5.

Parameters for the numerical solution.

ω1	q1(t)	−7.0E5 0.0 ≤ t ≤ 0.5; −5.6E6 t + 2.1E6 0.5 < t < 1.0; 5.6E6 t – 9.1E6 1.0 ≤ t ≤ 1.5; −7.0E5 1.5 ≤ t ≤ 4.0
ω2	q2(t)	−7.0E5 0.0 ≤ t ≤ 1.0; −5.6E6 t + 4.9E6 1.0 < t < 1.5; 5.6E6 t – 1.19E6 1.5 ≤ t ≤ 2.0; −7.0E5 2.0 ≤ t ≤ 4.0
ω3	q3(t)	−7.0E5 0.0 ≤ t ≤ 1.5; −5.6E6 t + 7.7E6 1.5 < t < 2.0; 5.6E6 t – 1.47E6 2.0 ≤ t ≤ 2.5; −7.0E5 2.5 ≤ t ≤ 4.0
ω4	q4(t)	−7.0E5 0.0 ≤ t ≤ 2.0; −5.6E6 t + 1.05E6 2.0 < t < 2.5; 5.6E6 t – 1.75E6 2.5 ≤ t ≤ 3.0; −7.0E5 3.0 ≤ t ≤ 4.0
ω5	q5(t)	−7.0E5 0.0 ≤ t ≤ 2.5; −5.6E6 t + 1.33E6 2.5 < t < 3.0; 5.6E6 t – 2.03E6 3.0 ≤ t ≤ 3.5; −7.0E5 3.5 ≤ t ≤ 4.0

Table 6.

Parameters of the distributed heterogeneous boundary condition.

The five functions that will distribute in the annular regions of the disk surface defined in Table 2 do consider to develop the case of analysis. The functions are determined as histories of heat flux density and initially present a constant period, later to reach a maximum value with a steep slope and then return to the same state of constant magnitude. Each function presents an initial period of different tendency such that the curves $q_i\left(t\right)$ overlap at some instant, as shown in Figure 15a and the relationships in Table 6. Formulating the DHCP concerning the system in Figure 13 defines the boundary condition of the opposite with a value $q\left(t\right)=0$ and an initial field of constant temperatures of 920°C.

The solution with the discretization proposed in Table 5 and with the position of the sensors defined in the places indicated in Table 3 is shown graphically in Figure 15b with cooling curves $Y_i\left(t\right)$. It can be noted that the system cools a maximum of 170°C and that the gradients generated in the mass of each region define the subsequent cooling trend of the entire system. It can also be noted that neither the sharp change points of the $q_i\left(t\right)$ original functions nor the maximums of each curve are accurately estimated with the corresponding ${\widehat{q}}_i\left(t\right)$ (see Figure 15c). However, when solving the DHCP with these new boundary conditions, a very acceptable adjustment of the ${\widehat{Y}}_i\left(t\right)$ cooling curves is observed (see Figure 15d), so the differences between these systems with the approximation may be due to the approximation criteria implied by the numerical methods for both the DHCP and IHCP solution.

2.4 Applications

One of the most impressive fields of application of inverse analysis in heat conduction problems does link to the aerospace industry, specifically in the race for mastery of trips to outer space. The knowledge of the heating conditions of the surfaces of the shuttles exposed to the situation of reentry into the atmosphere contributed to the design of coatings, selection of suitable materials, and cooling systems to counteract the extreme heating and consequent failures. In this way, the same concept is applied to systems at the level of the earth’s surface when considering nuclear fission and fusion reactors. In other processes of interest, the analysis is projected to the design and monitoring of reactors for extracting mineral substances and their subsequent chemical refining, as well as in the design processes of materials with optimal allotropic properties, either by sintering or thermosetting, where the control of the rate of temperature change is depended on to obtain optimal microstructural characteristics. In general, one can think of any other system where thermal boundary conditions are complex to determine by fluid dynamics analysis, radiant systems, or a combination thereof. However, a limitation occurs with the impossibility of instrumenting with sensors and implementing data acquisition systems in the advice of interest. Another stemmed from the ability of computer systems to process a large amount of data, which also determined the number of dimensions implicit in the analyses. However, at present, it has been minimized. The selection of the number of analysis dimensions will depend on the degree of detail that defines the characteristics of the boundary condition and the conditions of the phenomenon studied.

3. Conclusions

It is shown in the case studies that it is possible to analyze the phenomena of heat transfer by conduction in solids by formulating inverse problems using essential foundations of mathematical concepts. As it has been talked about in the course of the discussion, treating these problems requires an appropriate quality of information that will process to determine thermal boundary conditions as a primary objective. Other complications associated with the stability of the solution to the ill-posed problem can adopt different theories. However, the sequential process of estimating the function by implementing the least square’s objective function guarantees an efficient technique, such that only at least one recorded thermal history near the boundary of interest is needed to obtain an estimate. More invasive instrumentation would significantly affect the continuity of the materials. Therefore, with the knowledge bases about the phenomena that occur, correct criteria can be taken regarding the precision of the IHCP solution. In such a way that, for example, it can deduce that the analyses are highly reliable due to finding magnitudes of relatively high gradients in the systems under study. Estimating thermal boundary conditions through the IHCP solution helps to weigh the magnitudes associated with interaction with highly complex external conditions. Additionally, it is practical to define the boundary conditions of transient systems in terms of heat flux histories due to the state of the mathematical formulation that governs the problems. The information can be used as a validation parameter for modeling convective, radiant, or a combination of both media transfer phenomena. The amount of data currently represents a manageable setback due to the high capacity of computer systems.

Acknowledgments

Thanks are due to the Tecnológico Nacional de México, which, through the Instituto Tecnológico de Morelia, provided financial support for the publication of this material.

Conflict of interest

The authors declare no conflict of interest.

Nomenclature

$\beta$	index of the heat flux function
$\mathbf{\beta }$	heat flux function matrix
$\widehat{\mathit{\beta }}$	matrix of estimated functions of heat flux
$\partial$	partial derivative
$∆$	magnitude of increase
$\aleph$	thermal histories data index
$\rho$	density
$\theta$	angular coordinate
${\omega }_i$	radial domain magnitude
$a_i$	coefficient of systems of algebraic equations
$A$	area perpendicular to the direction of unidirectional flow
$\mathbf{A}$	boundary conditions individuality matrix
$a_{j\mathrm{u}}$	component of the sensitivity coefficient matrix matrices
$a_{j\mathrm{v}}$	component of the sensitivity coefficient matrix matrices
$a_{j\mathrm{w}}$	component of the sensitivity coefficient matrix matrices
$b_i$	coefficient of systems of algebraic equations
$c_i$	coefficient of systems of algebraic equations
$C_{\mathrm{p}}$	specific heat capacity
$C_{\mathrm{uv}}$	IHCP solution component
$d_i$	coefficient of systems of algebraic equations
$D_{\mathrm{w}}$	IHCP solution component
$e_i$	coefficient of systems of algebraic equations
$f_i$	coefficient of systems of algebraic equations
$\mathbf{F}$	function of the Dirichlet, Neumann or Robinson type boundary condition
$g_i$	coefficient of systems of algebraic equations
$h$	convective heat transfer coefficient
$i$	geometric discretization index/index of mathematical operators
$j$	sensor index/index of mathematical operators
$k$	thermal conductivity
$l$	index of the number of future time steps
ℓ	number of future time steps
$m$	number of sensors
$M$	reference to the calculation time instant
$n$	number of control volumes in radial direction
$N$	number of data from thermal histories
$p$	index of mathematical operators
$q$	heat flux function
$\widehat{q}$	estimated heat flux function
$\mathbf{q}$	heat flux function matrix
$\dot{q}$	heat generated or consumed per unit volume
$Q$	heat flow
$r$	radial coordinate
$R$	magnitude of the radius of cylindrical systems
$S$	least squares function (objective function)
$S_{\mathrm{Y}}$	standard deviation from thermal histories
$t$	variable time
$T$	variable temperature
$T_{\mathrm{s}}$	surface temperature
$T_{\infty }$	bulk fluid temperature
$T^{\ast }$	abbreviation for sensitivity coefficient definition
$\mathbf{T}$	temperature field
$u$	index of mathematical operators
$v$	index of mathematical operators
$V_i$	amount of control volume
$w$	geometric discretization index/index of mathematical operators
${\mathbf{x}}_j$	histories of sensitivity coefficients
$\mathbf{X}$	sensitivity coefficient matrix
$Y$	thermal history
$\mathbf{Y}$	thermal history matrix
$z$	axial coordinate