# Solving of Two-Dimensional Unsteady-State Heat-Transfer Inverse Problem Using Finite Difference Method and Model Prediction Control Method

Authors: Shoubin Wang; Rui Ni
Journal: Complexity
 (2019)

---

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


---

## Body

## 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.(1)∂2T∂x2+∂2T∂y2=1α∂T∂t∈Ω,t>t0T=T¯∈Γ1,t>t0-λ∂T∂n=q∈Γ2,t>t0λ∂T∂n=hT-Tf∈Γ3,t>t0T=T0∈Ω,t=t0whereΓ1  is the Dirichlet (first-type) boundary condition, Γ2 is the Neumann (second-type) boundary condition, Γ3  is the Robin (third-type) boundary condition, and Γ=Γ1+Γ2+Γ3 is the boundary of the whole region Ω. α  denotes the thermal diffusivity, and α=λ/cρ. c, ρ, and λ denote the specific heat, the density, and the heat conduction coefficient of the object, respectively.  T represents the temperature, T- is the temperature given by the Dirichlet boundary condition, Tf is the environment temperature, and T0 is the initial temperature. q refers to the heat flux, n refers to the boundary outer normal vector, and h 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 ofx-axis is Δx=xi+1-xi and that ofy-axis is Δy=yj+1-yj, obviously, xi=iΔx, yj=jΔy, and i,j=0,1,2,3….n(n=0,1,2,3…) is used to uniformly discretize the time domain t≧0 and the step length between two time moments Δt=tn+1-tn,tn=nΔt, where i,j⋯xi,yj,n⋯tn;Ti,jn⋯n and the temperature at node i,j in the time moment n is Txi,yj,tn.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 i,j at the time moment of n, the equation can be rewritten as(2)∂2T∂x2+∂2T∂y2i,jn=1α∂T∂ti,jnThe 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.(3)∂T∂ti,jn=Ti,jn+1-Ti,jnΔt+OΔtThe second-order partial differential in the left of the equal sign can be approximated by the central difference quotient.(4)∂2T∂x2i,jn=Ti+1,jn-2Ti,jn+Ti-1,jnΔx2+OΔx2(5)∂2T∂y2i,jn=Ti,j+1n-2Ti,jn+Ti,j-1nΔy2+OΔy2Substituting (3), (4), and (5) into (2), we can get the difference equation of (2).(6)Ti+1,jn-2Ti,jn+Ti-1,jnΔx2+Ti,j+1n-2Ti,jn+Ti,j-1nΔy2=1αTi,jn+1-Ti,jnΔtEquation (6) is the difference equation of heat conduction equation, and the truncation error [34] is OΔt+Δx2Δy2.Assuming thatΔx=Δy=Δ and substituting it into (6), we can obtain(7)Ti+1,jn+Ti-1,jn+Ti,j-1n-4-1F0Ti,jn=Ti,jn+1F0whereF0 refers to the Fourier coefficient and F0=αΔt/Δx2=αΔt/Δy2=αΔt/Δ2.The stability condition of explicit finite difference equation of two-dimensional unsteady-state heat conduction without internal heat source is in interior node,F0≤1/4; in boundary node, F0≤1/22+Bi; in boundary angular point, F0≤1/41+Bi.


### 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 isT¯/2; then the boundary node temperature remains at T¯.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 Figure2.Second-type boundary condition, i.e., the heat flow boundary, is given,-λ∂T/∂n=q. Setting boundary x=0 as the given heat flow boundary condition and keeping it stable, the second-type boundary condition can be expressed as -λ∂T/∂x=q. Using central difference quotient to replace the first-order partial differential equation(8)-λT2n-T2′n2Δ=q(8) is rewritten as(9)T2′n=T2n-2ΔλqNode 1 is changed into interior node and (9) is substituted into the interior node explicit difference equation of (7) to get(10)T2n+T2′n+T4n+T0n-4-1F0T1n=1F0T1n+1Substituting (9) into (10), we can obtain(11)T1n+1=F0T0n+2T2n+T4n-2Δλq+1F0-4T1nThird-type boundary condition, i.e., the heat transfer boundary, is given,λ∂T/∂n=hT-Tf. Setting boundary x=0 as the given heat transfer boundary condition and keeping it stable, the third-type boundary condition can be expressed as λ∂T/∂x=hT-Tf. Using central difference quotient to replace the first-order partial differential equation(12)λT2n-T2′n2Δ=hT1n-Tf(12) is rewritten as(13)T2′n=T2n-2BiT1n-TfNode 1 is changed into interior node and (13) is substituted into the interior node explicit difference equation of (7) to get(14)T2n+T2′n+T4n+T0n-4-1F0T1n=1F0T1n+1Substituting (13) into (14), we can obtain(15)T1n+1=F0T0n+2T2n+T4n+2BiTf+1F0-4-2BiT1nwhereBi is the Biot number, Bi=hΔ/λ; the truncation error of the second- and third-type boundary condition is OΔ2.Adiabatic boundary condition isT=0. Similarly, the second- and third-type boundary condition is(16)T1n+1=F0T0n+2T2n+T4n+1F0-4T1nBoundary 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 thex- and y-direction, respectively.(17)λT3n-T3′n2Δ=hT0n-TfλT1n-T1′n2Δ=hT0n-TfEquation (17) is rewritten as(18)T3′n=T3n-2BiT0n-TfT1′n=T1n-2BiT0n-TfNode 1 is changed into interior node and (17) and (18) are substituted into (7) to get(19)T3n+T3′n+T1n+T1′n-4-1F0T0n=1F0T0n+1T0n+1=2F0T1n+T3n+2BiTf+12F0-2-2BiT0nBy (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

Figure3 shows the model of two-dimensional unsteady-state heat conduction system without internal heat source. The rectangle plate in Figure 3 is adopted; boundary D1,D2,D3 is for heat insulation and D4 is the third-type boundary condition. h is the heat-transfer coefficient. Then, (1) can be changed by the corresponding mathematical model as follows:(20)∂2T∂x2+∂2T∂y2=1α∂T∂t∈Ω,t>t0-λ∂T∂x=0∈D1,t>t0-λ∂T∂x=0∈D2,t>t0-λ∂T∂y=0∈D3,t>t0-λ∂T∂y=hT-Tf∈D4,t>t0T=T0∈Ω,t=t0


### 2.5. Direct Problem Verification

Figures4, 5, 6, and 7 display the temperature field distribution when t=50,100,150,200s. 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 is0.2m. The heat conductivity coefficient λ=47W/m∗K, thermal diffusivity a=1.28∗10-6m2/s, initial temperature T0=20°C, environment temperature Tf=50°C, and heat-transfer coefficient h=2000W/m2∗K.


## 2.1. Mathematical Model

The mathematical model of two-dimensional unsteady-state heat conduction without internal heat source is expressed as follows.(1)∂2T∂x2+∂2T∂y2=1α∂T∂t∈Ω,t>t0T=T¯∈Γ1,t>t0-λ∂T∂n=q∈Γ2,t>t0λ∂T∂n=hT-Tf∈Γ3,t>t0T=T0∈Ω,t=t0whereΓ1  is the Dirichlet (first-type) boundary condition, Γ2 is the Neumann (second-type) boundary condition, Γ3  is the Robin (third-type) boundary condition, and Γ=Γ1+Γ2+Γ3 is the boundary of the whole region Ω. α  denotes the thermal diffusivity, and α=λ/cρ. c, ρ, and λ denote the specific heat, the density, and the heat conduction coefficient of the object, respectively.  T represents the temperature, T- is the temperature given by the Dirichlet boundary condition, Tf is the environment temperature, and T0 is the initial temperature. q refers to the heat flux, n refers to the boundary outer normal vector, and h 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 ofx-axis is Δx=xi+1-xi and that ofy-axis is Δy=yj+1-yj, obviously, xi=iΔx, yj=jΔy, and i,j=0,1,2,3….n(n=0,1,2,3…) is used to uniformly discretize the time domain t≧0 and the step length between two time moments Δt=tn+1-tn,tn=nΔt, where i,j⋯xi,yj,n⋯tn;Ti,jn⋯n and the temperature at node i,j in the time moment n is Txi,yj,tn.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 i,j at the time moment of n, the equation can be rewritten as(2)∂2T∂x2+∂2T∂y2i,jn=1α∂T∂ti,jnThe 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.(3)∂T∂ti,jn=Ti,jn+1-Ti,jnΔt+OΔtThe second-order partial differential in the left of the equal sign can be approximated by the central difference quotient.(4)∂2T∂x2i,jn=Ti+1,jn-2Ti,jn+Ti-1,jnΔx2+OΔx2(5)∂2T∂y2i,jn=Ti,j+1n-2Ti,jn+Ti,j-1nΔy2+OΔy2Substituting (3), (4), and (5) into (2), we can get the difference equation of (2).(6)Ti+1,jn-2Ti,jn+Ti-1,jnΔx2+Ti,j+1n-2Ti,jn+Ti,j-1nΔy2=1αTi,jn+1-Ti,jnΔtEquation (6) is the difference equation of heat conduction equation, and the truncation error [34] is OΔt+Δx2Δy2.Assuming thatΔx=Δy=Δ and substituting it into (6), we can obtain(7)Ti+1,jn+Ti-1,jn+Ti,j-1n-4-1F0Ti,jn=Ti,jn+1F0whereF0 refers to the Fourier coefficient and F0=αΔt/Δx2=αΔt/Δy2=αΔt/Δ2.The stability condition of explicit finite difference equation of two-dimensional unsteady-state heat conduction without internal heat source is in interior node,F0≤1/4; in boundary node, F0≤1/22+Bi; in boundary angular point, F0≤1/41+Bi.


## 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 isT¯/2; then the boundary node temperature remains at T¯.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 Figure2.Second-type boundary condition, i.e., the heat flow boundary, is given,-λ∂T/∂n=q. Setting boundary x=0 as the given heat flow boundary condition and keeping it stable, the second-type boundary condition can be expressed as -λ∂T/∂x=q. Using central difference quotient to replace the first-order partial differential equation(8)-λT2n-T2′n2Δ=q(8) is rewritten as(9)T2′n=T2n-2ΔλqNode 1 is changed into interior node and (9) is substituted into the interior node explicit difference equation of (7) to get(10)T2n+T2′n+T4n+T0n-4-1F0T1n=1F0T1n+1Substituting (9) into (10), we can obtain(11)T1n+1=F0T0n+2T2n+T4n-2Δλq+1F0-4T1nThird-type boundary condition, i.e., the heat transfer boundary, is given,λ∂T/∂n=hT-Tf. Setting boundary x=0 as the given heat transfer boundary condition and keeping it stable, the third-type boundary condition can be expressed as λ∂T/∂x=hT-Tf. Using central difference quotient to replace the first-order partial differential equation(12)λT2n-T2′n2Δ=hT1n-Tf(12) is rewritten as(13)T2′n=T2n-2BiT1n-TfNode 1 is changed into interior node and (13) is substituted into the interior node explicit difference equation of (7) to get(14)T2n+T2′n+T4n+T0n-4-1F0T1n=1F0T1n+1Substituting (13) into (14), we can obtain(15)T1n+1=F0T0n+2T2n+T4n+2BiTf+1F0-4-2BiT1nwhereBi is the Biot number, Bi=hΔ/λ; the truncation error of the second- and third-type boundary condition is OΔ2.Adiabatic boundary condition isT=0. Similarly, the second- and third-type boundary condition is(16)T1n+1=F0T0n+2T2n+T4n+1F0-4T1nBoundary 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 thex- and y-direction, respectively.(17)λT3n-T3′n2Δ=hT0n-TfλT1n-T1′n2Δ=hT0n-TfEquation (17) is rewritten as(18)T3′n=T3n-2BiT0n-TfT1′n=T1n-2BiT0n-TfNode 1 is changed into interior node and (17) and (18) are substituted into (7) to get(19)T3n+T3′n+T1n+T1′n-4-1F0T0n=1F0T0n+1T0n+1=2F0T1n+T3n+2BiTf+12F0-2-2BiT0nBy (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

Figure3 shows the model of two-dimensional unsteady-state heat conduction system without internal heat source. The rectangle plate in Figure 3 is adopted; boundary D1,D2,D3 is for heat insulation and D4 is the third-type boundary condition. h is the heat-transfer coefficient. Then, (1) can be changed by the corresponding mathematical model as follows:(20)∂2T∂x2+∂2T∂y2=1α∂T∂t∈Ω,t>t0-λ∂T∂x=0∈D1,t>t0-λ∂T∂x=0∈D2,t>t0-λ∂T∂y=0∈D3,t>t0-λ∂T∂y=hT-Tf∈D4,t>t0T=T0∈Ω,t=t0


## 2.5. Direct Problem Verification

Figures4, 5, 6, and 7 display the temperature field distribution when t=50,100,150,200s. 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 is0.2m. The heat conductivity coefficient λ=47W/m∗K, thermal diffusivity a=1.28∗10-6m2/s, initial temperature T0=20°C, environment temperature Tf=50°C, and heat-transfer coefficient h=2000W/m2∗K.


## 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 theD4 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., Δhk+P, is used to predict the temperature at S in the D4 boundary in the moment of R, i.e., Tk+P, where P=0,1,2…R-1 and R is the future time step.According to the principle of linear superposition [35], after loading R group of increment Δhk+P on system since the time moment of k, the temperature at S, i.e., Tk+P, is obtained.(21)Tk+P=T¯k+PΔhk=Δhk+1=⋯Δhk+r-1=0+∑ϕRΔhk+ϕR-1Δhk+1+⋯+ϕ1Δhk+R-1Equation (21) is changed into(22)Tk+P=T¯k+Phk=hk+1=⋯hk+r-1=0+∑ΔϕR-1hk+ΔϕR-2hk+1+⋯+Δϕ0hk+R-1Equation (22) can be reduced to(23)T=T¯h=0+AhWhere,(24)A=Δϕ000⋯0Δϕ1Δϕ00⋯0⋮⋮⋮⋮ΔϕR-1ΔϕR-2ΔϕR-3⋯Δϕ0,ΔϕY=ϕY+1-ϕYY=0,1,⋯R-1.The step response system function from the time momenttk to tk+R-1 is defined as the impact of heat-transfer coefficient on Tx,y,t. After the derivation of hk by (20), the step response equation can be obtained as follows.(25)∂2H∂x2+∂2H∂y2=1α∂H∂t∈Ω,tk≤t≤tk+R-1-λ∂H∂x=0∈D1,tk≤t≤tk+R-1-λ∂H∂x=0∈D2,tk≤t≤tk+R-1-λ∂H∂y=0∈D3,tk≤t≤tk+R-1-λ∂H∂x=Tk-Tf-hHx,y,t∈Γ3,tk≤t≤tk+R-1H=H0∈Ω,t=tkFrom (25), it can be seen that H is related to h. 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 h.The corresponding discrete valueHk+R-1 is obtained by (25); hence, the dynamic step response coefficient ϕR is further determined as(26)ϕR=Hk+R-1


### 3.2. Rolling Optimization of Inverse Problem

Measurement value and predictive value of temperature can be seen in the time range fromtk to tk+R-1. According to the finite optimization, R 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.(27)min⁡Jh=Tmea-TpreΤTmea-Tpre+hΤθhwhereθ is the regularization parameter matrix (28)θ=α0⋯00α⋯0⋮⋮⋱000⋯αand α is the regularization parameter.After the derivation ofh based on (27) and making dJh/dh=0, the optimal control rate can be obtained as(29)h=AΤA+θ-1AΤTmea-T¯preh=0The optimal heat-transfer coefficient at momentk can be obtained following (29).(30)h=RAΤA+θ-1AΤTmea-T¯preh=0whereR=1,0,…,0Τ. 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 atS in the k moment, TSk. Besides, in the case that the temperature measurement value at S, Tmeak sees measurement error, the temperature measurement value can be expressed by the actual temperature plus the measurement error.(31)Tmeak=Tactk+ωσwhereω is the random number of standard normal in the range and σ is the standard deviation of measurement value, which is expressed as(32)σ=1K-1∑x=1kTmeak-Tactk2where the constant K is the number of iterations.The residual of heat-transfer coefficient in the whole inversion time domain is defined as(33)Sα=1K-1∑k=1khactk-hSk2In (33), hactk and hSk are the actual value and inversion value of heat-transfer coefficient, respectively.SincehSk is unknown, it is available to calculate the temperature measurement value at S (TSk) using hSk with direct problem algorithm; thus, the temperature residual in the inversion time domain can be obtained.(34)Sα=1K-1∑k=1kTmeak-TS1k2In ideal condition,(35)hactk=hSkSimilarily,(36)TS1k=TSkFrom the residual principle, the regularization parameter is the optimal when both (35) and (36) are satisfied.(37)STα=σ


### 3.4. Solving Procedure of Inverse Problem

( 1 ) Select the initial predictive value of the heat-transfer coefficient hpk at a time moment to perform inversion.( 2 ) Obtain the temperature calculation values in measurement point S at R time moments after that moment based on hpk and (20).( 3 ) Calculate the optimal regularization parameter α based on (37).( 4 ) Assume the heat-transfer coefficient in the initial stage of inversion h=0, and obtain the step response matrix A by Eq. (25);( 5 ) Confirm the heat-transfer coefficient at the time moment hk according to (30) and then use the direct problem algorithm to reconstruct the temperature field when the heat-transfer coefficient is hk.( 6 ) Following the time direction backward, change the value in the initial stage of inversion and repeat steps (4) and (5); then get the inverse value of heat-transfer coefficient at different time moments.


## 3.1. Prediction Model of Inverse Problem

The step response of heat-transfer coefficient in theD4 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., Δhk+P, is used to predict the temperature at S in the D4 boundary in the moment of R, i.e., Tk+P, where P=0,1,2…R-1 and R is the future time step.According to the principle of linear superposition [35], after loading R group of increment Δhk+P on system since the time moment of k, the temperature at S, i.e., Tk+P, is obtained.(21)Tk+P=T¯k+PΔhk=Δhk+1=⋯Δhk+r-1=0+∑ϕRΔhk+ϕR-1Δhk+1+⋯+ϕ1Δhk+R-1Equation (21) is changed into(22)Tk+P=T¯k+Phk=hk+1=⋯hk+r-1=0+∑ΔϕR-1hk+ΔϕR-2hk+1+⋯+Δϕ0hk+R-1Equation (22) can be reduced to(23)T=T¯h=0+AhWhere,(24)A=Δϕ000⋯0Δϕ1Δϕ00⋯0⋮⋮⋮⋮ΔϕR-1ΔϕR-2ΔϕR-3⋯Δϕ0,ΔϕY=ϕY+1-ϕYY=0,1,⋯R-1.The step response system function from the time momenttk to tk+R-1 is defined as the impact of heat-transfer coefficient on Tx,y,t. After the derivation of hk by (20), the step response equation can be obtained as follows.(25)∂2H∂x2+∂2H∂y2=1α∂H∂t∈Ω,tk≤t≤tk+R-1-λ∂H∂x=0∈D1,tk≤t≤tk+R-1-λ∂H∂x=0∈D2,tk≤t≤tk+R-1-λ∂H∂y=0∈D3,tk≤t≤tk+R-1-λ∂H∂x=Tk-Tf-hHx,y,t∈Γ3,tk≤t≤tk+R-1H=H0∈Ω,t=tkFrom (25), it can be seen that H is related to h. 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 h.The corresponding discrete valueHk+R-1 is obtained by (25); hence, the dynamic step response coefficient ϕR is further determined as(26)ϕR=Hk+R-1


## 3.2. Rolling Optimization of Inverse Problem

Measurement value and predictive value of temperature can be seen in the time range fromtk to tk+R-1. According to the finite optimization, R 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.(27)min⁡Jh=Tmea-TpreΤTmea-Tpre+hΤθhwhereθ is the regularization parameter matrix (28)θ=α0⋯00α⋯0⋮⋮⋱000⋯αand α is the regularization parameter.After the derivation ofh based on (27) and making dJh/dh=0, the optimal control rate can be obtained as(29)h=AΤA+θ-1AΤTmea-T¯preh=0The optimal heat-transfer coefficient at momentk can be obtained following (29).(30)h=RAΤA+θ-1AΤTmea-T¯preh=0whereR=1,0,…,0Τ. 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 atS in the k moment, TSk. Besides, in the case that the temperature measurement value at S, Tmeak sees measurement error, the temperature measurement value can be expressed by the actual temperature plus the measurement error.(31)Tmeak=Tactk+ωσwhereω is the random number of standard normal in the range and σ is the standard deviation of measurement value, which is expressed as(32)σ=1K-1∑x=1kTmeak-Tactk2where the constant K is the number of iterations.The residual of heat-transfer coefficient in the whole inversion time domain is defined as(33)Sα=1K-1∑k=1khactk-hSk2In (33), hactk and hSk are the actual value and inversion value of heat-transfer coefficient, respectively.SincehSk is unknown, it is available to calculate the temperature measurement value at S (TSk) using hSk with direct problem algorithm; thus, the temperature residual in the inversion time domain can be obtained.(34)Sα=1K-1∑k=1kTmeak-TS1k2In ideal condition,(35)hactk=hSkSimilarily,(36)TS1k=TSkFrom the residual principle, the regularization parameter is the optimal when both (35) and (36) are satisfied.(37)STα=σ


## 3.4. Solving Procedure of Inverse Problem

( 1 ) Select the initial predictive value of the heat-transfer coefficient hpk at a time moment to perform inversion.( 2 ) Obtain the temperature calculation values in measurement point S at R time moments after that moment based on hpk and (20).( 3 ) Calculate the optimal regularization parameter α based on (37).( 4 ) Assume the heat-transfer coefficient in the initial stage of inversion h=0, and obtain the step response matrix A by Eq. (25);( 5 ) Confirm the heat-transfer coefficient at the time moment hk according to (30) and then use the direct problem algorithm to reconstruct the temperature field when the heat-transfer coefficient is hk.( 6 ) Following the time direction backward, change the value in the initial stage of inversion and repeat steps (4) and (5); 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 (Figure1) used in the above-mentioned direct problem is adopted. In the simulation example, the length Lx and width Ly of the plate is 0.2m. The heat conductivity coefficient λ=47W/m∗K, thermal diffusivity a=1.28∗10-6m2/s, initial temperature T0=20°C, environment temperature Tf=50°C, and heat-transfer coefficient h=2000W/m2∗K. The purpose is to obtain the actual heat-transfer coefficient of the boundary D4.Figure 1
Heat conduction model.Figure 2
Boundary node.Figure 3
The model of two-dimensional unsteady-state heat conduction system without internal heat source.Figure 4
The temperature field in the t=50s.Figure 5
The temperature field in the t=100s.Figure 6
The temperature field in the t=150s.Figure 7
The temperature field in the t=200s.Figure 8
The curve of the measuring point with time.


### 4.1. Impact When the Measurement Error Is Zero

Given the measurement errorσ=0.000, when the measurement point is in the L=0.1m of D4 boundary and the future time step R=5, the inversion result is shown in Figure 9.Figure 9
The heat transfer coefficient of the measuring point without error.Figure9 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 stepR=5 and the measurement point is in the L=0.1m of D4 boundary, the inversion results when the measurement error is σ=0.001, σ=0.005, and σ=0.01 are displayed in Figures 10, 11, and 12, respectively.Figure 10
The heat transfer coefficient of the measuring point withσ=0.001.Figure 11
The heat transfer coefficient of the measuring point withσ=0.005.Figure 12
The heat transfer coefficient of the measuring point withσ=0.01.According to Table1 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.Table 1
Relative average errors of inversion result under different measurement errors givenL=0.1m and R=5.
Measurement errorσ 0.001 0.005 0.01  Relative average errorη% 7.01% 11.29% 16.01%


### 4.3. Impact of Measurement Point Position

Given measurement errorσ=0.001 and future time step R=5, the inversion results when the measurement point is in the L=0.004m,0.008m of D4 boundary and when the measurement point is l=0.001m from the D4 boundary are shown in Figures 13, 14, and 15, respectively.Figure 13
The heat transfer coefficient of the measuring point at the boundaryL=0.008m.Figure 14
The heat transfer coefficient of the measuring point at the boundaryL=0.004m.Figure 15
The heat transfer coefficient of the measuring point at the distance from the boundaryl=0.001m.Analyze the contents of Table2 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.Table 2
Relative average errors of inversion result in different measurement point positions givenσ=0.001 and R=5.
Measurement point L = 0.004 m L = 0.008 m l = 0.001 m  Relative average errorη% 7.65% 7.06% 21.05%


## 4.1. Impact When the Measurement Error Is Zero

Given the measurement errorσ=0.000, when the measurement point is in the L=0.1m of D4 boundary and the future time step R=5, the inversion result is shown in Figure 9.Figure 9
The heat transfer coefficient of the measuring point without error.Figure9 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 stepR=5 and the measurement point is in the L=0.1m of D4 boundary, the inversion results when the measurement error is σ=0.001, σ=0.005, and σ=0.01 are displayed in Figures 10, 11, and 12, respectively.Figure 10
The heat transfer coefficient of the measuring point withσ=0.001.Figure 11
The heat transfer coefficient of the measuring point withσ=0.005.Figure 12
The heat transfer coefficient of the measuring point withσ=0.01.According to Table1 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.Table 1
Relative average errors of inversion result under different measurement errors givenL=0.1m and R=5.
Measurement errorσ 0.001 0.005 0.01  Relative average errorη% 7.01% 11.29% 16.01%


## 4.3. Impact of Measurement Point Position

Given measurement errorσ=0.001 and future time step R=5, the inversion results when the measurement point is in the L=0.004m,0.008m of D4 boundary and when the measurement point is l=0.001m from the D4 boundary are shown in Figures 13, 14, and 15, respectively.Figure 13
The heat transfer coefficient of the measuring point at the boundaryL=0.008m.Figure 14
The heat transfer coefficient of the measuring point at the boundaryL=0.004m.Figure 15
The heat transfer coefficient of the measuring point at the distance from the boundaryl=0.001m.Analyze the contents of Table2 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.Table 2
Relative average errors of inversion result in different measurement point positions givenσ=0.001 and R=5.
Measurement point L = 0.004 m L = 0.008 m l = 0.001 m  Relative average errorη% 7.65% 7.06% 21.05%


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



---

*Source: 7432138-2019-07-22.xml*