Abstract
The Inverse Heat Conduction Problem (IHCP) refers to the inversion of the internal characteristics or thermal boundary conditions of a heat transfer system by using other known conditions of the system and according to some information that the system can observe. It has been extensively applied in the fields of engineering related to heat-transfer measurement, such as the aerospace, atomic energy technology, mechanical engineering, and metallurgy. The paper adopts Finite Difference Method (FDM) and Model Predictive Control Method (MPCM) to study the inverse problem in the third-type boundary heat-transfer coefficient involved in the two-dimensional unsteady heat conduction system. The residual principle is introduced to estimate the optimized regularization parameter in the model prediction control method, thereby obtaining a more precise inversion result. Finite difference method (FDM) is adopted for direct problem to calculate the temperature value in various time quanta of needed discrete point as well as the temperature field verification by time quantum, while inverse problem discusses the impact of different measurement errors and measurement point positions on the inverse result. As demonstrated by empirical analysis, the proposed method remains highly precise despite the presence of measurement errors or the close distance of measurement point position from the boundary angular point angle.
1. Introduction
The Inverse Heat Conduction Problem is able to retrieve the unknown parameters such as boundary conditions, material thermophysical parameters [1–3], internal heat sources, and boundary geometry by measuring the temperature information at the boundary or at some point in the heat-transfer system [4, 5]. The research of inverse heat conduction problem has a very wide application background. It has been applied in almost all fields of scientific engineering, including power engineering, aerospace engineering, metallurgical engineering, biomedical engineering, mechanical manufacturing, chemical engineering, nuclear physics, material processing, geometry optimization of equipment, and nondestructive testing. In view of the inverse problem of heat-transfer, experts and scholars have done quite a lot of research [6–11]. 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 [12]. Beck raised the concept of sensitivity coefficient and introduced it into the inverse problem, thereby successfully obtaining the steady-state and unsteady-state heat conduction application [13–15]. Luo et al. proposed the decentralized fuzzy inference method applicable to unsteady IHCP by dispersion and coordination of measurement information on the time domain on the basis of researching steady IHCP by using the decentralized fuzzy inference method [16]. Qian et al. solved the unsteady IHCP by using the SFSM and conjugate gradient method, which sufficiently demonstrated effectiveness of these two methods, analyzed, and compared the advantages and disadvantages of these two methods [17]. Jian Su et al. solved the heat conduction inverse problem of one transient heat-transfer coefficient by employing Alifanov’s iterative regularization algorithm [18]. Blanc G et al. investigated the one-dimensional transient inverse problem, finding that residual principle can optimize the key parameter in the heat conduction problem [19]. Wu Zhaochun studied the measurement point arrangement pattern during the solving process of the two-dimensional heat conduction inverse problem with DSC method, accordingly making some rational suggestions regarding the measurement points [20]. Lesnic et al. identified the thermophysical parameters in one-dimensional transient heat conduction problems by using the BEM [21]. Ershova et al. used the residual error principles in the Tikhonov regularization method and completed the crystal phonon inspection identification tasks [22]. Zhao Luyao combined the particle swarm optimization (PSO) and conjugate gradient method and applied the combined method to the inversion of the heat-transfer coefficient of one-dimensional unsteady-state system. It has been reported that the method exhibits high preciseness [23]. Li et al. 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]. Li Yanhao resolved the heat-flow problem found in the two-dimensional transient inverse problem by using the model prediction control algorithm and inversion result was relatively precise [27]. Fan Jianxue also adopted the model prediction control algorithm to solve the heat-transfer coefficient in the inner wall of two-dimensional transient steam drum and achieved good calculation results [28].
Regarding the boundary heat transfer in the heat conduction system, in the paper, FDM is adopted to solve the direct problem of the two-dimensional unsteady-state heat conduction without internal heat source and model prediction control method is used to solve the inverse problem. Besides, residual principle is introduced to optimize the regularization parameter during the inversion process, thereby improving the efficiency of inversion in terms of speed and time.
2. Unsteady-State Direct Problem
The Inverse Heat Conduction Problem usually involves the multiple deduction of the forward problem, and its inversion accuracy is directly affected by the calculation accuracy of the forward problem. Positive problem refers to the solution of historical temperature field through given boundary conditions, initial temperature, and thermal conductivity differential equation. Common solutions are Lattice Boltzmann Method, Finite Volume Method, Adomain Decomposition Method, Boundary Element Method, and Finite Difference Method.
LBM (Lattice Boltzmann Method) [29] is a mesoscopic research method based on molecular kinetics, which can well describe the complex and small interfaces in porous media. It is widely used in small-scale numerical simulation of porous media and other objects with complex interface structures.
The basic idea of FVM (Finite Volume Method) [30] is to divide the computational region into a series of nonrepeated control volumes and make each grid point have a control volume around it. A set of discrete equations is obtained by integrating the differential equations to be resolved with each control volume. It is commonly used in the case of discrete and complex grids.
ADM (Adomain Decomposition Method) [31] is a research method that decomposes the true solution of an equation into the sum of the components of several solutions and then tries to find the components of the solutions and make the sum of the components of the solutions approximate the true solution with any desired high precision. But it can only obtain accurate results under the condition that the energy conservation equation is not nonlinear.
BEM (Boundary Element Method) [32] is a research method which divides elements on the boundary of the defined domain and approximates the boundary conditions by functions satisfying the governing equations. The basic advantage of BEM is that it can reduce the dimensionality of the problem, but when it comes to solving the basic solution of the problem, the process of solving the basic solution is generally complicated.
In this paper, a square rectangular plate is selected as an experimental physical model, which is a very common physical model. The finite-difference method [33] can reduce the amount of calculation required by other research methods for the positive problem, and it is also convenient to query the temperature change curves of the required measuring points in the negative problem.
2.1. Mathematical Model
The mathematical model of two-dimensional unsteady-state heat conduction without internal heat source is expressed as follows.
where is the Dirichlet (first-type) boundary condition, is the Neumann (second-type) boundary condition, is the Robin (third-type) boundary condition, and is the boundary of the whole region . denotes the thermal diffusivity, and . , , and denote the specific heat, the density, and the heat conduction coefficient of the object, respectively. represents the temperature, is the temperature given by the Dirichlet boundary condition, is the environment temperature, and is the initial temperature. refers to the heat flux, refers to the boundary outer normal vector, and refers to the surface heat-transfer coefficient.
2.2. Discretization and Difference Scheme
The discrete rules of two-dimensional unsteady-state heat conduction problem without internal heat source in geometry and time domain are as follows.
Assuming that after the domain of uniform discretization, the step length of x-axis is and that of y-axis is , obviously, , , and .
n() is used to uniformly discretize the time domain t≧0 and the step length between two time moments , where and the temperature at node in the time moment is .
The explicit difference array of the two-dimensional unsteady-state heat conduction without internal heat source is expressed as follows.
Applying the first heat conduction equation in (1) to node at the time moment of , the equation can be rewritten as
The partial differential in the two sides of (2) can be approximated by difference quotient. The temperature item in the right of the equal sign can be approximated by first-order time forward difference quotient.
The second-order partial differential in the left of the equal sign can be approximated by the central difference quotient.
Substituting (3), (4), and (5) into (2), we can get the difference equation of (2).
Equation (6) is the difference equation of heat conduction equation, and the truncation error [34] is .
Assuming that and substituting it into (6), we can obtain
where refers to the Fourier coefficient and .
The stability condition of explicit finite difference equation of two-dimensional unsteady-state heat conduction without internal heat source is in interior node, ; in boundary node, ; in boundary angular point, .
2.3. Boundary Conditions
First-type boundary condition, i.e., the temperature, is given. In general, when the FDM is used for calculation, it shall be processed as in the moment of the initial, the boundary node temperature is ; then the boundary node temperature remains at .
In the second- and third-type boundary condition, it is necessary to set virtual node outside the boundary to make the boundary node into interior node. The node numbering is shown in Figure 2.
Second-type boundary condition, i.e., the heat flow boundary, is given,. Setting boundary as the given heat flow boundary condition and keeping it stable, the second-type boundary condition can be expressed as . Using central difference quotient to replace the first-order partial differential equation
(8) is rewritten as
Node 1 is changed into interior node and (9) is substituted into the interior node explicit difference equation of (7) to get
Substituting (9) into (10), we can obtain
Third-type boundary condition, i.e., the heat transfer boundary, is given, . Setting boundary as the given heat transfer boundary condition and keeping it stable, the third-type boundary condition can be expressed as . Using central difference quotient to replace the first-order partial differential equation
(12) is rewritten as
Node 1 is changed into interior node and (13) is substituted into the interior node explicit difference equation of (7) to get
Substituting (13) into (14), we can obtain
where is the Biot number, ; the truncation error of the second- and third-type boundary condition is .
Adiabatic boundary condition is . Similarly, the second- and third-type boundary condition is
Boundary angular point is 0 node. Virtual nodes 1’ and 3’ are set in the symmetric position of node 1 and node 3, respectively, and the central different quotient is applied in the - and -direction, respectively.
Equation (17) is rewritten as
Node 1 is changed into interior node and (17) and (18) are substituted into (7) to get
By (7), (11), (15), (16), and (19), the temperature value in any point of the model can be obtained.
2.4. Mathematical Model about the Heat Transfer Process of Rectangular Plate
Figure 3 shows the model of two-dimensional unsteady-state heat conduction system without internal heat source. The rectangle plate in Figure 3 is adopted; boundary is for heat insulation and is the third-type boundary condition. is the heat-transfer coefficient. Then, (1) can be changed by the corresponding mathematical model as follows:
2.5. Direct Problem Verification
Figures 4, 5, 6, and 7 display the temperature field distribution when . Figure 8 is the curve of measuring points with time. The simulation result of direct problem solving can demonstrate the rationality of explicit finite difference, which is convenient for performing the inversion algorithm of inverse problem.
The length Lx and width Ly of the plate is . The heat conductivity coefficient , thermal diffusivity , initial temperature , environment temperature , and heat-transfer coefficient .
3. Unsteady-State Problem
Predictive control is a model-based control algorithm, which focuses on the function of the model rather than the form of the model. Compared with other control methods, its characteristics are reflected in the use of rolling optimization and rolling implementation of the control mode to achieve the control effect, but also did not give up the traditional control feedback. Therefore, the predictive control algorithm is based on the future dynamic behavior of the process model prediction system under a certain control effect, uses the rolling optimization to obtain the control effect under the corresponding constraint conditions and performance requirements, and corrects the prediction of future dynamic behavior in the rolling optimization process by detecting real-time information.
3.1. Prediction Model of Inverse Problem
The step response of heat-transfer coefficient in the boundary of direct problem model is taken as the prediction model of inverse problem. The increment of heat-transfer coefficient in the boundary in future time, i.e., , is used to predict the temperature at in the boundary in the moment of , i.e., , where and is the future time step.
According to the principle of linear superposition [35], after loading group of increment on system since the time moment of , the temperature at , i.e., , is obtained.
Equation (21) is changed into
Equation (22) can be reduced to
Where,
The step response system function from the time moment to is defined as the impact of heat-transfer coefficient on . After the derivation of by (20), the step response equation can be obtained as follows.
From (25), it can be seen that is related to . Therefore, while solving inverse problem, it is necessary to use the explicit difference algorithm used in the direct problem solving and keep updating the calculation of .
The corresponding discrete value is obtained by (25); hence, the dynamic step response coefficient is further determined as
3.2. Rolling Optimization of Inverse Problem
Measurement value and predictive value of temperature can be seen in the time range from to . According to the finite optimization, parameters to be inverted are obtained, so the predictive value can be as close as possible to the measurement value in the future time domain; at this moment, the quadratic performance index of system can be launched.
where is the regularization parameter matrix and is the regularization parameter.
After the derivation of based on (27) and making , the optimal control rate can be obtained as
The optimal heat-transfer coefficient at moment can be obtained following (29).
where . There is no need to preset the function form for the inverse problem calculation by the above optimization algorithm.
3.3. Regularization Parameter
The residual principle [36–38] is introduced to calculate the optimal regularization parameter, aiming to reduce the impact of measurement errors on the inversion results.
To invert the boundary heat-transfer coefficient, it is necessary to firstly solve the direct problem using the predictive value of heat-transfer coefficient, to get the temperature calculation at in the moment, . Besides, in the case that the temperature measurement value at , sees measurement error, the temperature measurement value can be expressed by the actual temperature plus the measurement error.
where is the random number of standard normal in the range and is the standard deviation of measurement value, which is expressed as
where the constant K is the number of iterations.
The residual of heat-transfer coefficient in the whole inversion time domain is defined as
In (33), and are the actual value and inversion value of heat-transfer coefficient, respectively.
Since is unknown, it is available to calculate the temperature measurement value at () using with direct problem algorithm; thus, the temperature residual in the inversion time domain can be obtained.
In ideal condition,
Similarily,
From the residual principle, the regularization parameter is the optimal when both (35) and (36) are satisfied.
3.4. Solving Procedure of Inverse Problem
Select the initial predictive value of the heat-transfer coefficient at a time moment to perform inversion.
Obtain the temperature calculation values in measurement point S at R time moments after that moment based on and (20).
Calculate the optimal regularization parameter based on (37).
Assume the heat-transfer coefficient in the initial stage of inversion , and obtain the step response matrix by Eq. (25);
Confirm the heat-transfer coefficient at the time moment according to (30) and then use the direct problem algorithm to reconstruct the temperature field when the heat-transfer coefficient is .
Following the time direction backward, change the value in the initial stage of inversion and repeat steps and ; then get the inverse value of heat-transfer coefficient at different time moments.
4. Numerical Experiment and Analysis
Numerical experiments are performed to validate whether the proposed method is effective, with the focus on analyzing the impact of different measurement errors and measurement point positions on the inversion result. Also, the inversion result obtained in the condition without measurement error is compared with the practical result, which verifies the precision of the proposed method.
The two-dimensional plate heat transfer model (Figure 1) used in the above-mentioned direct problem is adopted. In the simulation example, the length Lx and width Ly of the plate is . The heat conductivity coefficient , thermal diffusivity , initial temperature , environment temperature , and heat-transfer coefficient . The purpose is to obtain the actual heat-transfer coefficient of the boundary D4.
4.1. Impact When the Measurement Error Is Zero
Given the measurement error , when the measurement point is in the of D4 boundary and the future time step , the inversion result is shown in Figure 9.
Figure 9 displays that except the transitory vibration in the initial stage. The inversion value is basically identical to the practical value, demonstrating the effectiveness of the inversion algorithm.
4.2. Impact of Measurement Error
Given the future time step and the measurement point is in the of D4 boundary, the inversion results when the measurement error is , , and are displayed in Figures 10, 11, and 12, respectively.
According to Table 1 and Figures 10, 11, and 12, smaller relative measurement error contributes to better inversion results. And enlarging measurement error will worsen the inversion results and aggravate the fluctuation.
4.3. Impact of Measurement Point Position
Given measurement error and future time step , the inversion results when the measurement point is in the of D4 boundary and when the measurement point is from the D4 boundary are shown in Figures 13, 14, and 15, respectively.
Analyze the contents of Table 2 and Figures 13 and 14. The explicit FDM is used for direct problem, when the measurement point in boundary is closer to the boundary angular point, which, however, imposes a little impact on the inversion result. Despite the increased relative average error, the proposed method still exhibits a better ability to track the exact solution of heat-transfer coefficient and the inversion result is relatively precise. In Figure 15, considering that the position of the measuring point is 0.001m away from the boundary and the initial time temperature of the position is 20, the temperature cannot change for a period of time, so the inversion result fluctuates greatly in the initial stage and increases when the distance of measurement point position from the boundary angular point becomes farther.
5. Conclusion
The boundary heat-transfer coefficient of the two-dimensional unsteady heat conduction system is inversed by the FDM and model prediction control method. By solving and analyzing the algorithm example, it demonstrates that the proposed methods have higher accuracy in the inversion process. Model predictive control method focuses on the model function rather than the structural form, so that we only need to know the step response or impulse response of the object; we can directly get the prediction model and skip the derivation process. It absorbs the idea of optimization control and replaces global optimization by rolling time-domain optimization combined with feedback correction, which avoids a lot of calculation required by global optimization and constantly corrects the influence caused by uncertain factors in the system. At the same time, by discussing the impacts of error free, measuring point positions, and measuring errors on the results, it demonstrates that the obtained inversion results, except the early oscillation, can better represent the stability of the exact solution.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Authors’ Contributions
Shoubin Wang and Rui Ni contributed to developing the ideas of this research. All of the authors were involved in preparing this manuscript.
Acknowledgments
This work was financially supported by the National Key Foundation for Exploring Scientific Instrument of China (2013YQ470767), Tianjin Municipal Education Commission Project for Scientific Research Items (2017KJ059), and Tianjin Science and Technology Commissioner Project (18JCTPJC62200, 18JCTPJC64100).