Solving of Two-Dimensional Unsteady Inverse Heat Conduction Problems Based on Boundary Element Method and Sequential Function Specification Method

The boundary element method (BEM) and sequential function specification method (SFSM) are used to research the inverse problem of boundary heat flux identification in the two-dimensional heat conduction system. The future time step in the SFSM is optimized by introducing the residual error principles to get the more accurate inversion results. For the forward problems, the BEM is used to calculate the required temperature value of discrete point; for the inverse problems, the impacts of different future time steps, measuring point position, and measuring error on the inversion results are discussed. Furthermore, the comparison is made for the optimal future time step obtained by introducing the residual error principle and the inherent future time step. The example analysis shows that the method proposed still has higher accuracy when the measuring error exists or the measuring point position is far away from the boundary heat flux.


Introduction
The inverse heat conduction problems (IHCP) are to measure the temperature at the heat conduction system boundary or internal point or points by using the experimental method to obtain partial temperature information and inverse some unknown parameters: the boundary condition, material thermophysical parameter, internal heat source and boundary geometry, and so on [1][2][3][4].The IHCP researches have wide application background and are nearly applied in all fields of science engineering: the aerospace engineering, bioengineering, power engineering, machine manufacturing, chemical engineering, nuclear physics, metallurgy, material processing, equipment geometry optimization, nondestructive testing, and so on [5][6][7][8][9].The domestic and foreign scholars have made many researches on IHCP.Duda identified the heat flux in two-dimensional transient heat conduction and reconstructed the transient temperature field by utilizing the finite element method (FEM) and Levenberg-Marquardt method in ANSYS Multiphysics software.The method mentioned above was applied to the identification of aerodynamic heating on an atmospheric reentry capsule [10,11].Luo et al. proposed the decentralized fuzzy inference method applicable to unsteady IHCP by dispersion and coordination of measurement information on time domain on the basis of researching steady IHCP by using the decentralized fuzzy inference method [12,13].Qian et al. solved the unsteady IHCP by using the SFSM and conjugate gradient method, which sufficiently demonstrated the effectiveness of these two methods, analyzed, and compared the advantages and disadvantages of these two methods [14][15][16].Lin et al. proposed an improved SFSM and researched the IHCP with time-varying internal heat source [17].Cabeza et al. researched the one-dimensional transient inverse problems.It can be found that the time step is the key parameter in the SFSM [18].Shao et al. used the conjugate gradient methods and SFSM for heat flux inversion for the lower surface of the fixed geometric domain in heat flux reversion for the variable geometric domain and verified the stability and effectiveness of such method [19,20].Lesnic et al. identified the thermophysical parameters in one-dimensional transient heat conduction problems by using the BEM [21].Ershova and Sidikova used the residual error principles in the Tikhonov regularization method and completed the crystal phonon spectrum identification tasks [22].Weizhen proposed a subsection identification method, by which the IHCP identified by variable thermophysical parameters was solved.Furthermore, it is proved that the calculation accuracy and efficiency of such method are superior to the common optimization algorithm [23].Li and Liu researched IHCP by using the BEM and identified the irregular boundaries [24,25].Zhou et al. solved the heat conductivity coefficient in the two-dimensional transient inverse problems by using the BEM and gradient regularization method and obtained the relatively accurate inversion results [26].Yaparova solved the inverse heat conduction boundary value problems with stable boundary based on the Laplace and Fourier transformation method [27].
For the boundary heat flux identification problem in the heat conduction system, the BEM is used to solve the two-dimensional unsteady forward problem without internal heat source; the SFSM is used to solve the inverse problem.In the process of solving the inverse problem, the future time step in the inversion process is optimized by introducing the residual error principle to improve the inversion accuracy.

Unsteady Forward Problem
2.1.Boundary Integral Equation.The mathematical model of the two-dimensional unsteady heat conduction problem without internal heat source for isotropic bodies is described below: where Γ 1 is the first boundary condition, Γ 2 is the second boundary condition, Γ 3 is the third boundary condition, and Γ = Γ 1 + Γ 2 + Γ 3 is the boundary of the whole region.a is the heat diffusion rate a = λ/cρ, c is the specific heat capacity of the object, ρ is the density of the object, and λ is heat conductivity coefficient of the object.T is the temperature, q is the normal component of the heat flux vector at the boundary surface, and n is the coordinate along the external normal vector, q ′ = h/λ T f , where T f is the ambient environment temperature and h is the surface heat transfer coefficient.
The weight function T * is introduced by the weighted residual method to get [28] The left side of Formula ( 2) is decomposed into The first item of Formula ( 3) is converted into Use Green's second formula: D v∇ 2 u − u∇ 2 v d D = s v ∂u/∂n − u ∂v/∂n d s , where s is the boundary curve of region D and d s is the arc differential.
Based on Green's formula, Formula (4) is converted as Formula ( 5) is substituted into Formula (3), and Formula (3) is substituted into Formula (2) to get The single time integration by parts is conducted for Formula ( 7) is substituted into Formula (6) to get Because Γ = Γ 1 + Γ 2 + Γ 3 , ∂T/∂n = q, and ∂T * /∂n = q * , Formula (8) can be converted as Formula ( 10) is substituted into Formula (9) to get Further simplify and combine the second item and the third item at the right side of Formula (11) to get The basic solution of such equation is shown below: where d is the space dimension.The two-dimensional problem will be discussed in this paper; d = 2 and r = The basic solution has the following characteristics: Formula ( 13) is derived: where D = ∂r/∂n * r, where D is the vertical distance from origin point i to the boundary.Formulas ( 14) and ( 15) are substituted into Formula ( 12) and simplified to get where

Dispersion of Boundary Integral Equation.
When the time domain is divided, functions T and q will change with 3 Complexity time.But its change is small than the change of T * and q * and can be ignored.Therefore, it can be considered constant within a small time interval.The time integration by subsections can be made for Formula (16).
The internal time level is integrated: where E i b is the exponential integral function which can be calculated by a series: where C is Euler's constant C = 0 57721566 ; when 0 ≤ b ≤ 1, the approximate value of the first five items is taken generally.Based on the above formula, Formula (18) can be written as When the space domain is divided, the domain Ω is divided into M units, and the boundary Γ is divided into N units; when j is within the boundary Γ 1 , Γ 2 , the formula can be written as The linear interpolation is used in this paper.Because the interpolation function of the linear unit is

24
the boundary curve can be approximate to the straight line.In each unit, the value of T and q is taken on the endpoint of the unit and linearly approximated.Formula ( 23) is collated to get where Given Formula (25) is converted as Similarly, when j is within the boundary, Formula ( 22) can be written as Formula ( 29) is converted as Formulas ( 28) and ( 31) are written in a matrix form: The value of T and q at the boundary node can be obtained through Formula (32).Taking C i = 1, the temperature at any internal point can be obtained through Formulas ( 16), (25), and (29).

Mathematical Model of Rectangle Plate Heat Transfer
Process. Figure 1 is the model of a two-dimensional unsteady heat conduction system without internal heat source.The rectangle plate as shown in Figure 1 is used.The boundary D 1 , D 2 , D 3 is heat insulated.The boundary D 4 has the timedependent heat flux q t .By using the relevant mathematical model, Formula (1) can be converted as Among them, Duda has done a lot of experimental and theoretical research on a heat transfer model and heat conduction, which has achieved many valuable results [29][30][31].
Duda identified the transient temperature distribution and the local heat flux on unknown boundary edges with the help of a flat plate heat conduction system shown in Figures 2 and 3 and the inversion method in the literature.The determination of the heating process of a simple twodimensional plate (Figure 2) and the identification of simultaneous transient heat flow by conduction (Figure 3) can demonstrate high accuracy and stability of the proposed algorithm.Figure 2 is the model of square plate heat conduction.Figure 3 is the model of a rectangular cross-section of an infinitely long beam [10,11].

Unsteady Inverse Problem
3.1.Objective Function of Inverse Problem.The inverse problem corresponding to the forward problem is to inverse the unknown boundary heat flux q X at the present based on the temperature information of measuring point S in the region, other known boundary conditions, and thermophysical parameter in the forward problems.At time t X , the heat flux value q 1 , q 2 , q 3 , … , q X−1 at time X − 1 and the measured value T X , T X+1 , … , T X+R−1 of measuring point temperature at later time R are known.
Its corresponding objective function can be defined as below: 5 Complexity where q X is the parameter to be inversed.T X+n−1 q X is the calculated temperature value of measuring point S at time t X+n−1 and obtained by solving the forward problem through the predicted value of q X .T X+n−1 mea is the measured temperature value of measuring point S at time t X+n−1 .

Sequential Function Specification Method for Inverse
Problem.Suppose that the boundary heat flux satisfies the certain function relationship with the time domain t X , t X+R−1 : Given the initial predicted value of heat flux q X to be inversed as q X p , considering that the model described in Formula (33) has linear characteristics, the calculated value T X+n−1 of measuring point is T X+n−1 q X = T X+n−1 q X p + Z X+n−1 q X − q X p , 36 where Z X+n−1 = ∂T X+n−1 /∂q X is the sensitivity coefficient of heat flux to be inversed.q X is derived by the equation in Formula (33) and the corresponding sensitivity equation can be obtained:

37
For the sensitivity equation solving, as Z is unrelated to q X , the BEM for the forward problems is used to obtain the sensitivity coefficient.
Substitute Formula (36) into the objective function (34) and given ∂J q X /∂q X = 0 to get

38
Further simplify to get The residual error principle is introduced to calculate the optimal further time step and reduce the impacts of measuring error on inversion results [32][33][34].
The boundary heat flux value will be inversed in this paper.Therefore, the forward problem should be solved firstly by the predicted value of heat flux to obtain the calculated temperature value T x S of measuring point S at time x.In addition, when the measured temperature value T x mea of measuring point S has a measuring error, T x mea can be represented by the sum of the actual temperature and measuring error, namely, where ω is the standard normal random number in interval −2 576, 2 576 , ωσ is the measuring error, and σ is the standard deviation of the measured value.The form of σ is shown below:

Complexity
As the measured temperature and calculated temperature are known, the standard temperature error in the whole inversed time domain can be obtained.The standard temperature error formula can be defined below: In ideal conditions, there exits According to the residual error principles, when Formula ( 43) is satisfied, R B is the optimal future time step.The standard temperature error is equal to the standard deviation of the measured value, namely, The obtained future time step R B is substituted into the SFSM to calculate the heat flux to be inversed and thus reduce the impacts of the measuring error on the inversion results.

Inverse Problem Solving Process
Step 1. Select the initial predicted value q X p of heat flux at a certain time.
Step 2. Calculate to obtain the calculated temperature value of measuring point S at time R after such time through q X p value and Formula (33).
Step 3. Calculate the optimal future time step R B by Formula (44).
Step 4. Obtain the corresponding sensitivity ∂J q X /∂q X by Formula (37).
Step 5. Update the heat flux value q X to be inversed by Formula (39) to obtain the final inversion results.
Step 6. Push backward in time orientation, repeat Steps 1-5, and obtain the heat flux inversion value at different times.

Numerical Experiment and Analysis
The effectiveness of the above methods is verified by numerical experiment.The impacts of different future time steps, measuring point positions, and measuring error on the inversion results are analyzed.At the same time, the inversion results of the optimal future time step obtained through residual error principle and the inherent future time step are analyzed to verify the accuracy of proposed methods.
The two-dimensional plate heat transfer model for forward problems as shown in Figure 1 is used.In such simulation example, the length Lx and width Ly of the plate is 0.05 m, the heat conductivity coefficient is λ = 50 W/ m ⋅ K , the heat diffusion rate is a = 1 × 10 −5 m 2 /s, and the initial temperature is T 0 = 20 °C.The actual heat flux of boundary D 4 is distributed in step wave as shown in The value in Formula (45) is the exact solution of heat flux (exact).The temperature of measuring point S is obtained for actual heat flux distribution by the methods used in forward problems.The measured temperature of measuring point is obtained by Formula (40) for the above inverse problems.

Impacts of Future Time
Step on Results.In the case of standard deviation σ = 0 00 and the position of the measuring point L = 0 01 m away from the boundary of D 4 , the measured temperature history and the calculated temperature history when the future time step R = 2, 5, 8 are shown in Figure 4.The inversion results when the future time step R = 2, 5, 8 are shown in Figure 5. Figure 4 shows the history of measured and calculated temperatures for different future time steps.Figure 5 shows the impacts of different future time steps on the inversion results when σ = 0 00 and L = 0 01 m.
In the case of standard deviation σ = 0 001 and the position of the measuring point L = 0 01 m away from the boundary of D 4 , the inversion results when the future time step R = 2, 5, 8 are shown in Figure 6. Figure 6 shows the impacts of different future time steps on the inversion results when σ = 0 001 and L = 0 01 m.
Table 1 shows the relative average error of inversion results of different future time steps when σ = 0 00 and L = 0 01 m and σ = 0 001 and L = 0 01 m.Based on the data in Table 1 and by the comparison of the inversion results in Figures 5 and 6, when the future time step is increased, the relative average error of results is increased but the relative average error can be controlled within 9% by the proposed methods.In addition, when R value is smaller, the inversion results are sensitive to the measuring error.In case of different errors, the increase of R value can obtain the relatively smooth inverse value curve, control measuring error, and improve the inversion accuracy.

Impacts of Measuring Point Position on Results.
In the case of standard deviation σ = 0 00 and the future time step R = 2, the measured temperature history and the calculated temperature history when the position of the measuring point L = 0 01 m, 0 015 m, 0 02 m, 0 025 m are shown in Figure 7.The inversion results when the position of the measuring point L = 0 01 m, 0 015 m, 0 02 m, 0 025 m are shown in Figure 8. Figure 7 shows the history of measured and 7 Complexity calculated temperatures for different measuring point positions.Figure 8 shows the impacts of different measuring point positions on the inversion results when σ = 0 00 and R = 2.
Table 2 shows the relative average error of inversion results of different measuring point positions when σ = 0 00 and R = 2. Based on the data in Table 2 and the inversion results in Figure 8, the different measuring point positions have little impacts on the inversion results.Although the relative average error of result is increased, the proposed methods can track the exact solution of heat flux and obtain more accurate inversion results.In Figure 8, when the heat fluxes change, the inversion results will fluctuate.Furthermore, the farther the measuring point position, the larger the fluctuation.It is caused by damping and delay of unsteady heat conduction.8 Complexity measured temperature history and the calculated temperature history when the standard deviation σ = 0 001, 0 005, 0 01, 0 02 are shown in Figure 9.The inversion results when the standard deviation σ = 0 001, 0 005, 0 01, 0 02 are shown in Figure 10. Figure 9 shows the history of measured and calculated temperatures for different measuring errors.
Figure 10 shows the impacts of different measuring errors on the inversion results when R = 3 and L = 0 01 m.Table 3 shows the relative average error of inversion results of different measuring errors when L = 0 01 m and R = 3.Based on the data in Table 3 and the inversion results in Figure 10, when the measuring error is relatively smaller, the better inversion results can be obtained.When increasing the measuring error, the inversion results will be deteriorated and vibration will be more severe.
When there exists a higher measuring error, the residual error principle is introduced.For different measuring errors, the methods for the above inverse problem are used to calculate the optimal future time step.Taking the measuring point position L = 0 01 m away from the boundary of D 4 , when the standard deviation σ = 0 001, 0 005, 0 01, 0 02, the corresponding optimal future time step obtained from Formula (44) is R B = 3, 5, 6, 8.The inversion results are shown in Figure 11. Figure 11 shows the impacts of different measuring errors on the inversion results when L = 0 01 m and R = R B .
Table 4 shows the relative average error of inversion results of different measuring errors when L = 0 01 m and R = R B .By comparing the data in Table 4 and inversion results in Figure 11 with the data in Table 3 and inversion results in Figure 10, it can be known that the proposed methods can effectively restrain the impacts of the measuring error on inversion results and control the relative average error of inversion results within 9% when there exists the measuring error.

Conclusion
The boundary heat flux of the two-dimensional unsteady heat conduction system is inversed by the BEM and SFSM based on residual error principles.By solving and analyzing the algorithm example, it demonstrates that the proposed methods have higher accuracy in the inversion process.At the same time, by discussing the impacts of different future time steps, measuring point positions, and measuring errors on the results, it demonstrates that the obtained inversion results can better represent the variation trend of the exact solution on time orientation and have better stability when there exists the measuring error or the measuring point position is changed.

Figure 3 :
Figure 3: The model of a rectangular cross-section of an infinitely long beam.

25 LFigure 7 :Figure 4 :
Figure 7: The history of measured and calculated temperatures for different measuring point positions when σ = 0 00 and R = 2.

Figure 11 :
Figure 11: Impacts of different measuring errors on inversion results when L = 0 01 m and R = R B .

Table 1 :
Relative average error of inversion results of different future time steps when σ = 0 00 and L = 0 01 m.

Table 3 :
Relative average error of inversion results of different measuring errors when L = 0 01 m and R = 3.

Table 4 :
Relative average error of inversion results of different measuring errors when L = 0 01 m and R = R B .