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

Authors: Shoubin Wang; Yuanzheng Deng; Xiaogang Sun
Journal: Complexity
 (2018)

---

## Abstract

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.


---

## Body

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


## 2. 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:(1)∂2T∂x2+∂2T∂y2=1a∂T∂t,∈Ω,t>t0,T=T¯,∈Γ1,t>t0,q=−λ∂T∂n=q¯,∈Γ2,t>t0,∂T∂n+hλT=q¯′,∈Γ3,t>t0,T=T0,∈Ω,t=t0,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/λTf, where Tf is the ambient environment temperature and h is the surface heat transfer coefficient.The weight functionT∗ is introduced by the weighted residual method to get [28]
(2)∫t0tx∫Ω∇2T−1a∂T∂tT∗dΩdτ=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3∂T∂n+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯∂T∗∂ndΓdτ.The left side of Formula (2) is decomposed into
(3)∫t0tx∫Ω∇2T−1a∂T∂tT∗dΩdτ=∫t0tx∫ΩT∗∇2TdΩdτ−∫t0tx∫Ω1a∂T∂tT∗dΩdτ.The first item of Formula (3) is converted into
(4)∫t0tx∫ΩT∗∇2TdΩdτ=∫t0tx∫ΩT∗∇2T−T∇2T∗dΩdτ+∫t0tx∫ΩT∇2T∗dΩdτ,Use Green’s second formula:∫Dv∇2u−u∇2vdD=∫sv∂u/∂n−u∂v/∂nds, where s is the boundary curve of region D and ds is the arc differential.Based on Green’s formula, Formula (4) is converted as
(5)∫t0tx∫ΩT∗∇2TdΩdτ=∫t0tx∫ΓT∗∂T∂n−T∂T∗∂ndΓdτ+∫t0tx∫ΩT∇2T∗dΩdτ.Formula (5) is substituted into Formula (3), and Formula (3) is substituted into Formula (2) to get
(6)∫t0tx∫ΓT∗∂T∂n−T∂T∗∂ndΓdτ+∫t0tx∫ΩT∇2T∗dΩdτ−∫t0tx∫Ω1a∂T∂tT∗dΩdτ=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3∂T∂n+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯∂T∗∂ndΓdτ.The single time integration by parts is conducted for∫t0tx∫Ω1/a∂T/∂tT∗dΩdτ∫udv=uv−∫vdu to get
(7)∫t0tx∫Ω1a∂T∂tT∗dΩdτ=∫Ω1aT∗TdΩt=t0t=tx−∫t0tx∫Ω1a∂T∗∂tTdΩdτ.Formula (7) is substituted into Formula (6) to get
(8)∫t0tx∫ΩT∇2T∗dΩdτ+∫t0tx∫Ω1a∂T∗∂tTdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3∂T∂n+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯∂T∗∂ndΓdτ−∫t0tx∫ΓT∗∂T∂n−T∂T∗∂ndΓdτ.BecauseΓ=Γ1+Γ2+Γ3, ∂T/∂n=q, and ∂T∗/∂n=q∗, Formula (8) can be converted as
(9)∫t0tx∫ΩT∇2T∗+1a∂T∗∂tdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3q+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯q∗dΓdτ−∫t0tx∫ΓT∗q−Tq∗dΓdτ.∫t0tx∫ΓT∗q−Tq∗dΓdτ is decomposed to get
(10)∫t0tx∫ΓT∗q−Tq∗dΓdτ=∫t0tx∫Γ1T∗q−Tq∗dΓdτ+∫t0tx∫Γ2T∗q−Tq∗dΓdτ+∫t0tx∫Γ3T∗q−Tq∗dΓdτ.Formula (10) is substituted into Formula (9) to get
(11)∫t0tx∫ΩT∇2T∗+1a∂T∗∂tdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ3hλTT∗dΓdτ+∫t0tx∫Γ1T¯q∗dΓdτ+∫t0tx∫Γ2Tq∗dΓdτ+∫t0tx∫Γ3Tq∗dΓdτ−∫t0tx∫Γ1qT∗dΓdτ+∫t0tx∫Γ2q¯T∗dΓdτ+∫t0tx∫Γ3q¯′T∗dΓdτ.Further simplify and combine the second item and the third item at the right side of Formula (11) to get
(12)∫t0tx∫ΩT∇2T∗+1a∂T∗∂tdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ3hλTT∗dΓdτ+∫t0tx∫Γq∗TdΓdτ−∫t0tx∫ΓT∗qdΓdτ.The basic solution of such equation is shown below:(13)T∗=14πatx−td/2exp−r24atx−t,where d is the space dimension. The two-dimensional problem will be discussed in this paper; d=2 and r=x−xi2+y−yi2. The basic solution has the following characteristics:
(14)∇2T∗+1a∂T∗∂t=0,∈Ω,t<tx,∫ΩT∗TdΩ=Ti,t=tx.Formula (13) is derived:
(15)q∗=∂T∗∂n=−r8πa2tx−t2exp−r24atx−t∗∂r∂n=−D8πa2tx−t2exp−r24atx−t,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
(16)1αCiTi+∫t0tx∫Γq∗TdΓdτ+∫t0tx∫Γ3hλTT∗dΓdτ=∫t0tx∫ΓT∗qdΓdτ+∫Ω1aT∗TdΩt=t0,where
(17)Ci=1,∈Ω,12,∈Γ,Smoothboundary.


### 2.2. Dispersion of Boundary Integral Equation

When the time domain is divided, functionsT and q will change with 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).
(18)CiTi+a∫ΓT∫t0txq∗dτdΓ+∫Γ3hλT∫t0txT∗dτdΓ=a∫Γq∫t0txT∗dτdΓ+∫ΩT∗TdΩt=t0.The internal time level is integrated:(19)∫t0txq∗dτ=−D2πar2exp−r24atx−t0,∫t0txT∗dτ=14πaEib,where
(20)b=r24atx−t0,where Eib is the exponential integral function which can be calculated by a series:
(21)Eib=−C−lnb+∑k=1∞−1k−1bkk⋅k!,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
(22)CiTitx+a∫ΓTtxqt∗dΓ+a∫Γ3hλTtxTt∗dΓ=a∫ΓqtxTt∗dΓ+∫ΩT∗Tt0dΩ.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
(23)CiTitx+a∑j=1N1+N2∫ΓjTjtxqt∗dΓ=a∑j=1N1+N2∫ΓjqjtxTt∗dΓ+∑m=1M∫ΩMT∗Tt0dΩ.The linear interpolation is used in this paper. Because the interpolation function of the linear unit is(24)φ1ξ=1−ξ2,φ2ξ=1+ξ2,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
(25)CiTitx+∑j=1N1+N2hij1Tjtx+hij2Tj+1tx=∑j=1N1+N2gij1qjtx+gij2qj+1tx+∑m=1M∫ΩMT∗Tt0dΩ,where
(26)hijy=a∫Γjφyqt∗dΓ,gijy=a∫ΓjφyTt∗dΓ,y=1,2.Given(27)Hijy=hij−12+hij1,Gijy=gij−12+gij1.Formula (25) is converted as
(28)CiTitx+∑j=1N1+N2HijTjtx=∑j=1N1+N2Gijqjtx+∑m=1M∫ΩMT∗Tt0dΩ.Similarly, whenj is within the boundary, Formula (22) can be written as
(29)CiTitx+∑j=1N3hij1+hλgij1Tjtx+hij2+hλgij2Tj+1tx=∑j=1N3gij1qjtx+gij2qj+1tx+∑m=1M∫ΩMT∗Tt0dΩ.Given(30)Hijy=hij−12+hλgij−12+hij1+hλgij1,Gijy=gij−12+gij1.Formula (29) is converted as
(31)CiTitx+∑j=1N3HijTjtx=∑j=1N3Gijqjtx+∑m=1M∫ΩMT∗Tt0dΩ.Formulas (28) and (31) are written in a matrix form:
(32)HTtx=GQtx+P,P=∑m=1M∫ΩMT∗Tt0dΩ.The value ofT and q at the boundary node can be obtained through Formula (32). Taking Ci=1, the temperature at any internal point can be obtained through Formulas (16), (25), and (29).


### 2.3. Mathematical Model of Rectangle Plate Heat Transfer Process

Figure1 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 D1,D2,D3 is heat insulated. The boundary D4 has the time-dependent heat flux qt. By using the relevant mathematical model, Formula (1) can be converted as
(33)∂2T∂x2+∂2T∂y2=1a∂T∂t,∈Ω,t>t0,−λ∂T∂x=0,∈D1,t>t0,−λ∂T∂x=0,∈D2,t>t0,−λ∂T∂y=0,∈D3,t>t0,−λ∂T∂y=qt,∈D4,t>t0,T=T0,∈Ω,t=t0.Figure 1
Two-dimensional unsteady plate heat conduction system.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–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 Figures2 and 3 and the inversion method in the literature. The determination of the heating process of a simple two-dimensional 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].Figure 2
The model of square plate heat conduction.Figure 3
The model of a rectangular cross-section of an infinitely long beam.


## 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:(1)∂2T∂x2+∂2T∂y2=1a∂T∂t,∈Ω,t>t0,T=T¯,∈Γ1,t>t0,q=−λ∂T∂n=q¯,∈Γ2,t>t0,∂T∂n+hλT=q¯′,∈Γ3,t>t0,T=T0,∈Ω,t=t0,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/λTf, where Tf is the ambient environment temperature and h is the surface heat transfer coefficient.The weight functionT∗ is introduced by the weighted residual method to get [28]
(2)∫t0tx∫Ω∇2T−1a∂T∂tT∗dΩdτ=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3∂T∂n+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯∂T∗∂ndΓdτ.The left side of Formula (2) is decomposed into
(3)∫t0tx∫Ω∇2T−1a∂T∂tT∗dΩdτ=∫t0tx∫ΩT∗∇2TdΩdτ−∫t0tx∫Ω1a∂T∂tT∗dΩdτ.The first item of Formula (3) is converted into
(4)∫t0tx∫ΩT∗∇2TdΩdτ=∫t0tx∫ΩT∗∇2T−T∇2T∗dΩdτ+∫t0tx∫ΩT∇2T∗dΩdτ,Use Green’s second formula:∫Dv∇2u−u∇2vdD=∫sv∂u/∂n−u∂v/∂nds, where s is the boundary curve of region D and ds is the arc differential.Based on Green’s formula, Formula (4) is converted as
(5)∫t0tx∫ΩT∗∇2TdΩdτ=∫t0tx∫ΓT∗∂T∂n−T∂T∗∂ndΓdτ+∫t0tx∫ΩT∇2T∗dΩdτ.Formula (5) is substituted into Formula (3), and Formula (3) is substituted into Formula (2) to get
(6)∫t0tx∫ΓT∗∂T∂n−T∂T∗∂ndΓdτ+∫t0tx∫ΩT∇2T∗dΩdτ−∫t0tx∫Ω1a∂T∂tT∗dΩdτ=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3∂T∂n+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯∂T∗∂ndΓdτ.The single time integration by parts is conducted for∫t0tx∫Ω1/a∂T/∂tT∗dΩdτ∫udv=uv−∫vdu to get
(7)∫t0tx∫Ω1a∂T∂tT∗dΩdτ=∫Ω1aT∗TdΩt=t0t=tx−∫t0tx∫Ω1a∂T∗∂tTdΩdτ.Formula (7) is substituted into Formula (6) to get
(8)∫t0tx∫ΩT∇2T∗dΩdτ+∫t0tx∫Ω1a∂T∗∂tTdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3∂T∂n+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯∂T∗∂ndΓdτ−∫t0tx∫ΓT∗∂T∂n−T∂T∗∂ndΓdτ.BecauseΓ=Γ1+Γ2+Γ3, ∂T/∂n=q, and ∂T∗/∂n=q∗, Formula (8) can be converted as
(9)∫t0tx∫ΩT∇2T∗+1a∂T∗∂tdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ2q−q¯T∗dΓdτ+∫t0tx∫Γ3q+hλT−q¯′T∗dΓdτ−∫t0tx∫Γ1T−T¯q∗dΓdτ−∫t0tx∫ΓT∗q−Tq∗dΓdτ.∫t0tx∫ΓT∗q−Tq∗dΓdτ is decomposed to get
(10)∫t0tx∫ΓT∗q−Tq∗dΓdτ=∫t0tx∫Γ1T∗q−Tq∗dΓdτ+∫t0tx∫Γ2T∗q−Tq∗dΓdτ+∫t0tx∫Γ3T∗q−Tq∗dΓdτ.Formula (10) is substituted into Formula (9) to get
(11)∫t0tx∫ΩT∇2T∗+1a∂T∗∂tdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ3hλTT∗dΓdτ+∫t0tx∫Γ1T¯q∗dΓdτ+∫t0tx∫Γ2Tq∗dΓdτ+∫t0tx∫Γ3Tq∗dΓdτ−∫t0tx∫Γ1qT∗dΓdτ+∫t0tx∫Γ2q¯T∗dΓdτ+∫t0tx∫Γ3q¯′T∗dΓdτ.Further simplify and combine the second item and the third item at the right side of Formula (11) to get
(12)∫t0tx∫ΩT∇2T∗+1a∂T∗∂tdΩdτ−∫Ω1aT∗TdΩt=t0t=tx=∫t0tx∫Γ3hλTT∗dΓdτ+∫t0tx∫Γq∗TdΓdτ−∫t0tx∫ΓT∗qdΓdτ.The basic solution of such equation is shown below:(13)T∗=14πatx−td/2exp−r24atx−t,where d is the space dimension. The two-dimensional problem will be discussed in this paper; d=2 and r=x−xi2+y−yi2. The basic solution has the following characteristics:
(14)∇2T∗+1a∂T∗∂t=0,∈Ω,t<tx,∫ΩT∗TdΩ=Ti,t=tx.Formula (13) is derived:
(15)q∗=∂T∗∂n=−r8πa2tx−t2exp−r24atx−t∗∂r∂n=−D8πa2tx−t2exp−r24atx−t,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
(16)1αCiTi+∫t0tx∫Γq∗TdΓdτ+∫t0tx∫Γ3hλTT∗dΓdτ=∫t0tx∫ΓT∗qdΓdτ+∫Ω1aT∗TdΩt=t0,where
(17)Ci=1,∈Ω,12,∈Γ,Smoothboundary.


## 2.2. Dispersion of Boundary Integral Equation

When the time domain is divided, functionsT and q will change with 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).
(18)CiTi+a∫ΓT∫t0txq∗dτdΓ+∫Γ3hλT∫t0txT∗dτdΓ=a∫Γq∫t0txT∗dτdΓ+∫ΩT∗TdΩt=t0.The internal time level is integrated:(19)∫t0txq∗dτ=−D2πar2exp−r24atx−t0,∫t0txT∗dτ=14πaEib,where
(20)b=r24atx−t0,where Eib is the exponential integral function which can be calculated by a series:
(21)Eib=−C−lnb+∑k=1∞−1k−1bkk⋅k!,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
(22)CiTitx+a∫ΓTtxqt∗dΓ+a∫Γ3hλTtxTt∗dΓ=a∫ΓqtxTt∗dΓ+∫ΩT∗Tt0dΩ.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
(23)CiTitx+a∑j=1N1+N2∫ΓjTjtxqt∗dΓ=a∑j=1N1+N2∫ΓjqjtxTt∗dΓ+∑m=1M∫ΩMT∗Tt0dΩ.The linear interpolation is used in this paper. Because the interpolation function of the linear unit is(24)φ1ξ=1−ξ2,φ2ξ=1+ξ2,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
(25)CiTitx+∑j=1N1+N2hij1Tjtx+hij2Tj+1tx=∑j=1N1+N2gij1qjtx+gij2qj+1tx+∑m=1M∫ΩMT∗Tt0dΩ,where
(26)hijy=a∫Γjφyqt∗dΓ,gijy=a∫ΓjφyTt∗dΓ,y=1,2.Given(27)Hijy=hij−12+hij1,Gijy=gij−12+gij1.Formula (25) is converted as
(28)CiTitx+∑j=1N1+N2HijTjtx=∑j=1N1+N2Gijqjtx+∑m=1M∫ΩMT∗Tt0dΩ.Similarly, whenj is within the boundary, Formula (22) can be written as
(29)CiTitx+∑j=1N3hij1+hλgij1Tjtx+hij2+hλgij2Tj+1tx=∑j=1N3gij1qjtx+gij2qj+1tx+∑m=1M∫ΩMT∗Tt0dΩ.Given(30)Hijy=hij−12+hλgij−12+hij1+hλgij1,Gijy=gij−12+gij1.Formula (29) is converted as
(31)CiTitx+∑j=1N3HijTjtx=∑j=1N3Gijqjtx+∑m=1M∫ΩMT∗Tt0dΩ.Formulas (28) and (31) are written in a matrix form:
(32)HTtx=GQtx+P,P=∑m=1M∫ΩMT∗Tt0dΩ.The value ofT and q at the boundary node can be obtained through Formula (32). Taking Ci=1, the temperature at any internal point can be obtained through Formulas (16), (25), and (29).


## 2.3. Mathematical Model of Rectangle Plate Heat Transfer Process

Figure1 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 D1,D2,D3 is heat insulated. The boundary D4 has the time-dependent heat flux qt. By using the relevant mathematical model, Formula (1) can be converted as
(33)∂2T∂x2+∂2T∂y2=1a∂T∂t,∈Ω,t>t0,−λ∂T∂x=0,∈D1,t>t0,−λ∂T∂x=0,∈D2,t>t0,−λ∂T∂y=0,∈D3,t>t0,−λ∂T∂y=qt,∈D4,t>t0,T=T0,∈Ω,t=t0.Figure 1
Two-dimensional unsteady plate heat conduction system.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–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 Figures2 and 3 and the inversion method in the literature. The determination of the heating process of a simple two-dimensional 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].Figure 2
The model of square plate heat conduction.Figure 3
The model of a rectangular cross-section of an infinitely long beam.


## 3. 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 fluxqX 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 tX, the heat flux value q1,q2,q3,…,qX−1 at time X−1 and the measured value TX,TX+1,…,TX+R−1 of measuring point temperature at later time R are known.Its corresponding objective function can be defined as below:(34)JqX=∑n=1RTX+n−1qX−TmeaX+n−12,where qX is the parameter to be inversed. TX+n−1qX is the calculated temperature value of measuring point S at time tX+n−1 and obtained by solving the forward problem through the predicted value of qX. TmeaX+n−1 is the measured temperature value of measuring point S at time tX+n−1.


### 3.2. Sequential Function Specification Method for Inverse Problem

Suppose that the boundary heat flux satisfies the certain function relationship with the time domaintX,tX+R−1:
(35)qX=qX+1=qX+2=⋯=qX+R−1.Given the initial predicted value of heat fluxqX to be inversed as qpX, considering that the model described in Formula (33) has linear characteristics, the calculated value TX+n−1 of measuring point is
(36)TX+n−1qX=TX+n−1qpX+ZX+n−1qX−qpX,where ZX+n−1=∂TX+n−1/∂qX is the sensitivity coefficient of heat flux to be inversed. qX is derived by the equation in Formula (33) and the corresponding sensitivity equation can be obtained:
(37)∂2Z∂x2+∂2Z∂y2=1a∂Z∂t,∈Ω,tX≤t≤tX+R−1,−λ∂Z∂x=0,∈D1,tX≤t≤tX+R−1,−λ∂Z∂x=0,∈D2,tX≤t≤tX+R−1,−λ∂Z∂y=0,∈D3,tX≤t≤tX+R−1,−λ∂Z∂y=1,∈D4,tX≤t≤tX+R−1,Z=Z0=0,∈Ω,t=tX.For the sensitivity equation solving, asZ is unrelated to qX, the BEM for the forward problems is used to obtain the sensitivity coefficient.Substitute Formula (36) into the objective function (34) and given ∂JqX/∂qX=0 to get
(38)∑n=1RZX+n−12qX−qpX=∑n=1RTX+n−1qX−TmeaX+n−1ZX+n−1.Further simplify to get(39)qX=qpX+∑n=1RTX+n−1qX−TmeaX+n−1ZX+n−1∑n=1RZX+n−12.


### 3.3. Residual Error Principle

The residual error principle is introduced to calculate the optimal further time step and reduce the impacts of measuring error on inversion results [32–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 valueTSx of measuring point S at time x. In addition, when the measured temperature value Tmeax of measuring point S has a measuring error, Tmeax can be represented by the sum of the actual temperature and measuring error, namely,
(40)Tmeax=Tactx+ωσ,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:
(41)σ=1X−1∑x=1XTmeax−Tactx2.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:(42)RT=1X−1∑x=1XTmeax−TSxR2.In ideal conditions, there exits(43)Tactx=TSxR.According to the residual error principles, when Formula (43) is satisfied, RB is the optimal future time step. The standard temperature error is equal to the standard deviation of the measured value, namely,
(44)RT=σ.The obtained future time stepRB 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.


### 3.4. Inverse Problem Solving Process

Step 1.
Select the initial predicted valueqpX of heat flux at a certain time.Step 2.
Calculate to obtain the calculated temperature value of measuring pointS at time R after such time through qpX value and Formula (33).Step 3.
Calculate the optimal future time stepRB by Formula (44).Step 4.
Obtain the corresponding sensitivity∂JqX/∂qX by Formula (37).Step 5.
Update the heat flux valueqX to be inversed by Formula (39) to obtain the final inversion results.Step 6.
Push backward in time orientation, repeat Steps1–5, and obtain the heat flux inversion value at different times.


## 3.1. Objective Function of Inverse Problem

The inverse problem corresponding to the forward problem is to inverse the unknown boundary heat fluxqX 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 tX, the heat flux value q1,q2,q3,…,qX−1 at time X−1 and the measured value TX,TX+1,…,TX+R−1 of measuring point temperature at later time R are known.Its corresponding objective function can be defined as below:(34)JqX=∑n=1RTX+n−1qX−TmeaX+n−12,where qX is the parameter to be inversed. TX+n−1qX is the calculated temperature value of measuring point S at time tX+n−1 and obtained by solving the forward problem through the predicted value of qX. TmeaX+n−1 is the measured temperature value of measuring point S at time tX+n−1.


## 3.2. Sequential Function Specification Method for Inverse Problem

Suppose that the boundary heat flux satisfies the certain function relationship with the time domaintX,tX+R−1:
(35)qX=qX+1=qX+2=⋯=qX+R−1.Given the initial predicted value of heat fluxqX to be inversed as qpX, considering that the model described in Formula (33) has linear characteristics, the calculated value TX+n−1 of measuring point is
(36)TX+n−1qX=TX+n−1qpX+ZX+n−1qX−qpX,where ZX+n−1=∂TX+n−1/∂qX is the sensitivity coefficient of heat flux to be inversed. qX is derived by the equation in Formula (33) and the corresponding sensitivity equation can be obtained:
(37)∂2Z∂x2+∂2Z∂y2=1a∂Z∂t,∈Ω,tX≤t≤tX+R−1,−λ∂Z∂x=0,∈D1,tX≤t≤tX+R−1,−λ∂Z∂x=0,∈D2,tX≤t≤tX+R−1,−λ∂Z∂y=0,∈D3,tX≤t≤tX+R−1,−λ∂Z∂y=1,∈D4,tX≤t≤tX+R−1,Z=Z0=0,∈Ω,t=tX.For the sensitivity equation solving, asZ is unrelated to qX, the BEM for the forward problems is used to obtain the sensitivity coefficient.Substitute Formula (36) into the objective function (34) and given ∂JqX/∂qX=0 to get
(38)∑n=1RZX+n−12qX−qpX=∑n=1RTX+n−1qX−TmeaX+n−1ZX+n−1.Further simplify to get(39)qX=qpX+∑n=1RTX+n−1qX−TmeaX+n−1ZX+n−1∑n=1RZX+n−12.


## 3.3. Residual Error Principle

The residual error principle is introduced to calculate the optimal further time step and reduce the impacts of measuring error on inversion results [32–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 valueTSx of measuring point S at time x. In addition, when the measured temperature value Tmeax of measuring point S has a measuring error, Tmeax can be represented by the sum of the actual temperature and measuring error, namely,
(40)Tmeax=Tactx+ωσ,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:
(41)σ=1X−1∑x=1XTmeax−Tactx2.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:(42)RT=1X−1∑x=1XTmeax−TSxR2.In ideal conditions, there exits(43)Tactx=TSxR.According to the residual error principles, when Formula (43) is satisfied, RB is the optimal future time step. The standard temperature error is equal to the standard deviation of the measured value, namely,
(44)RT=σ.The obtained future time stepRB 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.


## 3.4. Inverse Problem Solving Process

Step 1.
Select the initial predicted valueqpX of heat flux at a certain time.Step 2.
Calculate to obtain the calculated temperature value of measuring pointS at time R after such time through qpX value and Formula (33).Step 3.
Calculate the optimal future time stepRB by Formula (44).Step 4.
Obtain the corresponding sensitivity∂JqX/∂qX by Formula (37).Step 5.
Update the heat flux valueqX to be inversed by Formula (39) to obtain the final inversion results.Step 6.
Push backward in time orientation, repeat Steps1–5, and obtain the heat flux inversion value at different times.


## 4. 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 Figure1 is used. In such simulation example, the length Lx and width Ly of the plate is 0.05 m, the heat conductivity coefficient is λ=50W/m⋅K, the heat diffusion rate is a=1×10−5m2/s, and the initial temperature is T0=20°C. The actual heat flux of boundary D4 is distributed in step wave as shown in
(45)qt=0Wm,0<t<60,240<t<300,4000Wm2,60≤t<120,180<t≤240,8000Wm2,120≤t≤180.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.


### 4.1. Impacts of Future Time Step on Results

In the case of standard deviationσ=0.00 and the position of the measuring point L=0.01m away from the boundary of D4, 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.01m.Figure 4
The history of measured and calculated temperatures for different future time steps whenσ=0.00 and L=0.01m.Figure 5
Impacts of different future time steps on inversion results whenσ=0.00 and L=0.01m.In the case of standard deviationσ=0.001 and the position of the measuring point L=0.01m away from the boundary of D4, 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.01m.Figure 6
Impacts of different future time steps on inversion results whenσ=0.00 and L=0.01m.Table1 shows the relative average error of inversion results of different future time steps when σ=0.00 and L=0.01m and σ=0.001 and L=0.01m. 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.Table 1
Relative average error of inversion results of different future time steps whenσ=0.00 and L=0.01m.
Future time stepR
2
5
8
σ=0.00 relative average error η (%)
3.52
6.19
8.15
σ=0.001 relative average error η (%)
5.10
6.24
8.17


### 4.2. 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.01m,0.015m,0.02m,0.025m are shown in Figure 7. The inversion results when the position of the measuring point L=0.01m,0.015m,0.02m,0.025m are shown in Figure 8. Figure 7 shows the history of measured and 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.Figure 7
The history of measured and calculated temperatures for different measuring point positions whenσ=0.00 and R=2.Figure 8
Impacts of different measuring point positions on inversion results whenσ=0.00 and R=2.Table2 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.Table 2
Relative average error of inversion results of different measuring point positions whenσ=0.00 and R=2.
Distance from the boundaryL (m)
0.010
0.015
0.020
0.025
Relative average errorη (%)
3.52
3.53
3.83
5.09


### 4.3. Impacts of Measuring Error on Results

In the case of the future time stepR=3 and the position of the measuring point L=0.01m away from the boundary of D4, the 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.01m.Figure 9
The history of measured and calculated temperatures for different measuring errors whenL=0.01m and R=3.Figure 10
Impacts of different measuring errors on inversion results whenL=0.01m and R=3.Table3 shows the relative average error of inversion results of different measuring errors when L=0.01m 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.Table 3
Relative average error of inversion results of different measuring errors whenL=0.01m and R=3.
Standard deviationσ
0.001
0.005
0.01
0.02
Relative average errorη (%)
4.81
9.41
16.21
31.12When 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 positionL=0.01m away from the boundary of D4, when the standard deviation σ=0.001,0.005,0.01,0.02, the corresponding optimal future time step obtained from Formula (44) is RB=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.01m and R=RB.Figure 11
Impacts of different measuring errors on inversion results whenL=0.01m and R=RB.Table4 shows the relative average error of inversion results of different measuring errors when L=0.01m and R=RB. 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.Table 4
Relative average error of inversion results of different measuring errors whenL=0.01m and R=RB.
Standard deviationσ
0.001
0.005
0.01
0.02
Relative average errorη (%)
4.81
6.55
7.56
8.97


## 4.1. Impacts of Future Time Step on Results

In the case of standard deviationσ=0.00 and the position of the measuring point L=0.01m away from the boundary of D4, 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.01m.Figure 4
The history of measured and calculated temperatures for different future time steps whenσ=0.00 and L=0.01m.Figure 5
Impacts of different future time steps on inversion results whenσ=0.00 and L=0.01m.In the case of standard deviationσ=0.001 and the position of the measuring point L=0.01m away from the boundary of D4, 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.01m.Figure 6
Impacts of different future time steps on inversion results whenσ=0.00 and L=0.01m.Table1 shows the relative average error of inversion results of different future time steps when σ=0.00 and L=0.01m and σ=0.001 and L=0.01m. 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.Table 1
Relative average error of inversion results of different future time steps whenσ=0.00 and L=0.01m.
Future time stepR
2
5
8
σ=0.00 relative average error η (%)
3.52
6.19
8.15
σ=0.001 relative average error η (%)
5.10
6.24
8.17


## 4.2. 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.01m,0.015m,0.02m,0.025m are shown in Figure 7. The inversion results when the position of the measuring point L=0.01m,0.015m,0.02m,0.025m are shown in Figure 8. Figure 7 shows the history of measured and 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.Figure 7
The history of measured and calculated temperatures for different measuring point positions whenσ=0.00 and R=2.Figure 8
Impacts of different measuring point positions on inversion results whenσ=0.00 and R=2.Table2 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.Table 2
Relative average error of inversion results of different measuring point positions whenσ=0.00 and R=2.
Distance from the boundaryL (m)
0.010
0.015
0.020
0.025
Relative average errorη (%)
3.52
3.53
3.83
5.09


## 4.3. Impacts of Measuring Error on Results

In the case of the future time stepR=3 and the position of the measuring point L=0.01m away from the boundary of D4, the 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.01m.Figure 9
The history of measured and calculated temperatures for different measuring errors whenL=0.01m and R=3.Figure 10
Impacts of different measuring errors on inversion results whenL=0.01m and R=3.Table3 shows the relative average error of inversion results of different measuring errors when L=0.01m 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.Table 3
Relative average error of inversion results of different measuring errors whenL=0.01m and R=3.
Standard deviationσ
0.001
0.005
0.01
0.02
Relative average errorη (%)
4.81
9.41
16.21
31.12When 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 positionL=0.01m away from the boundary of D4, when the standard deviation σ=0.001,0.005,0.01,0.02, the corresponding optimal future time step obtained from Formula (44) is RB=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.01m and R=RB.Figure 11
Impacts of different measuring errors on inversion results whenL=0.01m and R=RB.Table4 shows the relative average error of inversion results of different measuring errors when L=0.01m and R=RB. 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.Table 4
Relative average error of inversion results of different measuring errors whenL=0.01m and R=RB.
Standard deviationσ
0.001
0.005
0.01
0.02
Relative average errorη (%)
4.81
6.55
7.56
8.97


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



---

*Source: 6741632-2018-08-29.xml*