Estimation of a Line Heat Source Using an Adjoint Free Gradient Based Inverse Analysis

1  Introduction

Due to the ever increasing power and capacity of computers as well as the development of different numerical methods over the past decades, inverse heat transfer problems have been subject of intensive research. These problems are mathematically classified as ill-posed as their solutions are unstable and very sensitive to the errors in the measurements used in the analysis [1]. Among the numerical methods to overcome the instabilities present in these problems are iterative regularization techniques which are based on the improving the inverse problem solution sequentially through an iterative minimization method such as steepest descent and conjugate gradient method and then stopping the inverse analysis process when a stable solution is obtained [2]. The determination of the functional form of an unknown quantity (such as heat flux, heat transfer coefficient, and heat source) in the inverse heat transfer problems is a difficult task as the problem falls into the category of function estimation approach [3,4]. In the inverse heat transfer problems with no a priori information on the functional form, the function estimation approach is usually concerned with the solution of two additional problems, namely, sensitivity and adjoint problems. The solution of these two problems requires additional mathematical manipulation and imposes the computational cost comparable to direct problem solution one on the inverse analysis [1].

The estimation of the heat source using an inverse analysis has received considerable attention over past decades [5–7]. It has many practical applications in various engineering fields such as welding processes [8,9]. In [10], the unknown functional form of a timewise varying strength of a line heat source is estimated using the conjugate gradient method. The heat source is placed at the specific location in a regular 2D physical domain with insulated boundaries. In [11], the time-varying strength of a heat source in a two-dimensional heat conduction problem is estimated using a method based on the Karhunen-Loeve Galerkin procedure and the conjugate gradient method. In [12], based on the finite element method and the conjugate gradient, both the location and the time-varying strength of point heat sources within a body are estimated from temperature boundary measurements. In [13], a genetic algorithm is used to estimate plane heat source in an one-dimensional domain with insulated boundaries. In [8], a multiple model adaptive inverse method is used to estimate a moving heat source. In [14], a meshless numerical scheme is presented for recovering the time-dependent heat source in general three-dimensional heat conduction problems. The second-order Crank–Nicolson scheme is employed for discretization in time-domain. In [15], a meshless numerical scheme is presented for solving the inverse heat source problem using the fundamental solution of the heat equation as a basis function. A regularized solution is obtained by employing the Tikhonov regularization method. In [16], the inverse problem of determination of the time-dependent heat source under a non-classical dynamic boundary condition has been considered. The numerical solution of the inverse problem includes the boundary element method combined with the second-order Tikhonov regularization. In [17], the time-varying heat source in one-dimensional heat conduction is estimated using conjugate gradient method combined with the adjoint problem for three types of source variations including constant, linearly increasing, and linearly decreasing ones. In [18], unknown parameters of a time-varying plane heat source in a one-dimensional inverse heat conduction problem are estimated using Repulsive Particle Swarm Optimization (RPSO) method. Then the results are compared with those of the Levenberg-Marquardt Method, a widely used minimization method in the inverse heat conduction problems. In [19], unknown heat source strength in three-dimensional inverse heat conduction problems is studied. The direct problem is solved by applying the Green’s function method. The Tikhonov and truncated singular value decomposition regularization methods are developed to identify the unknown source strength.

In the literature, however, there are still some limitations on the methods used to identify the variables involved in the inverse heat transfer problems such as the heat source. For example, most of the heated body shapes considered are regular (rectangular or circular) or most of the boundary conditions considered include a constant temperature (Dirichlet boundary condition) or insulated case. Thus a methodology considering a general irregular 2D domain with a variety of boundary conditions is required.

The numerical method proposed in this study deals with 2D irregular physical domains and a variety of boundary conditions. It can also be extended to 3D problems and other heat source types such as point and plane heat sources. In this method, the irregular physical domain and all associated boundary conditions are transformed from the physical domain to the computational one so that the finite difference method can be used to solve the heat conduction equation for the irregular body shapes. A very easy to implement sensitivity analysis is proposed based on chain rule and direct differentiation which enables us to calculate all elements of the sensitivity matrix during the direct problem solution (at each iteration) without dealing with the sensitivity and adjoint problems to calculate the gradient of defined least squares objective function. After computation of the sensitivity matrix, the gradient of the objective function and search step size can be computed easily through explicit expressions. The steepest descent method, an iterative gradient based minimization method, is used to minimize the objective function. In this study, both the steady-state (for a constant heat source, a parameter estimation problem) and transient (for a time-varying heat source, a function estimation problem) heat conduction problems are considered. In the function estimation part, it is assumed that there is no prior information on the functional form of the heat source and the estimation process can be performed with a reasonable initial guess for the heat source.

2  Governing Equation

2.1 Transient Heat Conduction

Consider an irregular heat-conducting body of thermal conductivity kT, density ρ, and specific heat c which is initially at the temperature T0 (Fig. 1). At time t>0, a heat flux q˙ is applied at boundary surface Γ1, a convective heat transfer is applied at the boundaries Γi,i=2,3,4 with corresponding heat transfer coefficients hi,i=2,3,4 and surrounding temperatures T∞i,i=2,3,4, and heat is generated continuously over time at the point Lx†,y† with a line heat source of strength gL(t)(Wm). The transient heat transfer equation can be expressed as

kT(∂2T(x,y,t)∂x2+∂2T(x,y,t)∂y2)+g(x,y,t)=ρc∂T(x,y,t)∂t(1)

Figure 1: Irregular heat-conducting body (physical domain) subjected to a heat flux q˙ on surface Γ1, a line heat source at point Lx†,y†, and convective heat transfer on surfaces Γi,i=2,3,4 (a) and the corresponding computational domain (b)

The line heat source of time-varying strength gL(t)(Wm) is related to the volumetric source g(x,y,t)(Wm3) by the delta function notation as [20]

g(x,y,t)=gL(t)δ(x−x†)δ(y−y†)(2)

where t is the time and δ(−) represents the Dirac delta function and approximates the source in the domain and is defined by [11]

δn(x−x†)≈n2cosh2⁡(n(x−x†)),δn(y−y†)≈n2cosh2⁡(n(y−y†))(3)

where (x†,y†) is the location of time-varying line heat source. The function δn becomes the Dirac delta function as n approaches infinity. In this study, we take n=80. Each of the delta functions in Eq. (2) has the units of (1m). Thus the unit for the volumetric source g(x,y,t) is as (Wm)(1m)(1m)=(Wm3). Eq. (1) can be expressed as

kT(∂2T(x,y,t)∂x2+∂2T(x,y,t)∂y2)+gL(t)δn(x−x†)δn(y−y†)=ρc∂T(x,y,t)∂tin the physical domain Ω(x,y),t>0(4)

with the boundary and initial conditions

∂T(x,y,t)∂n1=q˙kTon boundary surface Γ1(x,y)(5)

∂T(x,y,t)∂ni=−hikT(TΓi(x,y,t)−T∞i)on boundary surface Γi(x,y),i=2,3,4(6)

T(x,y,0)=T0(x,y)in physical domain Ω(x,y)(7)

First the heat-conducting body is mapped from the irregular x,y physical domain (Fig. 1a) to the regular ξ,η computational domain (Fig. 1b). The elliptic grid generation method is used to mesh the irregular physical domain. Then by transforming the heat conduction equation and the boundary and initial conditions from the (x,y,t) to the (ξ,η,t) variables, we obtain [21]

kT(α∂2T(ξ,η,t)∂ξ2−2β∂2T(ξ,η,t)∂ξ∂η+γ∂2T(ξ,η,t)∂η2J2+(∇2ξ)∂T(ξ,η,t)∂ξ+(∇2η)∂T(ξ,η,t)∂η)+

gL(t)n24cosh2⁡(n(x−x†))cosh2⁡(n(y−y†))=ρc∂T(ξ,η,t)∂t(8)

where ∇2ξ=P(ξ,η) and ∇2η=Q(ξ,η) are grid control functions. By considering P(ξ,η)=Q(ξ,η)=0 to obtain a smooth grid over the physical domain, we obtain

kT(α∂2T(ξ,η,t)∂ξ2−2β∂2T(ξ,η,t)∂ξ∂η+γ∂2T(ξ,η,t)∂η2J2)+

n2gL(t)4cos⁡h2(n(x−x†))cosh2⁡(n(y−y†))=ρc∂T(ξ,η,t)∂t in 1<ξ<M,1<η<N, for t>0(9)

where

α=xη2+yη2

β=xξxη+yξyη

γ=xξ2+yξ2

J=xξyη−xηyξ(Jacobian of transformation)(10)

The boundary and initial conditions can be transformed as

(−1Jγ(γ∂T(ξ,η,t)∂η−β∂T(ξ,η,t)∂ξ))Γ1=q˙kT at 1<ξ<M,η=1, for t>0(11)

(1Jγ(γ∂T(ξ,η,t)∂η−β∂T(ξ,η,t)∂ξ))Γ2=−h2kT(T(ξ,η,t)−T∞2) at 1<ξ<M,η=N, for t>0(12)

(−1Jα(α∂T(ξ,η,t)∂ξ−β∂T(ξ,η,t)∂η))Γ3=−h3kT(T(ξ,η,t)−T∞3) at 1<η<N,ξ=1, for t>0(13)

(1Jα(α∂T(ξ,η,t)∂ξ−β∂T(ξ,η,t)∂η))Γ4=−h4kT(T(ξ,η,t)−T∞4) at 1<η<N,ξ=M, for t>0(14)

T(ξ,η,0)=T0∗(ξ,η) in 1<ξ<M,1<η<N, for t=0(15)

Using the finite difference method to discretize the derivatives in the resulting equations in the computational domain, we obtain (assuming Δξ=Δη=1)

fξ=12(fi+1,j−fi−1,j)

fη=12(fi,j+1−fi,j−1)

fξξ=fi+1,j−2fi,j+fi−1,j

fηη=fi,j+1−2fi,j+fi,j−1

fξη=14(fi+1,j+1−fi−1,j+1−fi+1,j−1+fi−1,j−1)(16)

where f≡x,y,T. To discretize the boundary condition equations, one may use one-sided forward and one-sided backward expressions. Using the forward-time-central-space (FTCS) scheme (an explicit method), we obtain [21]

kTJ2(α(Ti+1,jn−2Ti,jn+Ti−1,jn)−2β(Ti+1,j+1n−Ti−1,j+1n−Ti+1,j−1n+Ti−1,j−1n)4+γ(Ti,j+1n−2Ti,jn+Ti,j−1n))+

n2gL(t)4cos⁡h2(n(x−x†))cosh2⁡(n(y−y†))=ρcTi,jn+1−Ti,jnΔt,i=2,…,M−1,j=2,…,N−1, for t>0(17)

where Δt is the time step. The time-marching method is used to calculate Ti,jn+1, the nodal temperatures at the time level n+1, by knowing the nodal temperatures at the previous time level n, Ti,jn. Thus we obtain the following explicit expression for the transient heat conduction equation

Ti,jn+1=Ti,jn+Δtρcn2gL(t)4cosh2⁡(n(x−x†))cosh2⁡(n(y−y†))+

kTΔtρcJ2(α(Ti+1,jn−2Ti,jn+Ti−1,jn)−2β(Ti+1,j+1n−Ti−1,j+1n−Ti+1,j−1n+Ti−1,j−1n)4+γ(Ti,j+1n−2Ti,jn+Ti,j−1n))(18)

2.2 Steady-State Heat Conduction

With the same explanations for transforming the transient heat condition equation, we can obtain the transformed relations for the steady-state heat conduction equation. Then by substituting the finite difference expressions given in Eq. (16) into the obtained steady-state equation, and then solving for Ti,j we obtain (i=2,…,M−1;j=2,…,N−1)

Ti,j=12(α+γ)(αTi+1,j+αTi−1,j−2βTξη+γTi,j+1+γTi,j−1+J2n2gLkT4cos⁡h2(n(x−x†))cosh2⁡(n(y−y†)))(19)

where from Eq. (16), Tξη=1/4(Ti+1,j+1−Ti−1,j+1−Ti+1,j−1+Ti−1,j−1). An iterative method such as Gauss–Seidel or Successive Over-Relaxation (SOR) method can be used to solve the steady-state heat conduction equation, Eq. (19), and the boundary conditions.

3  Inverse Analysis

3.1 Objective Function for Transient Heat Conduction

As in this study, two separate inverse heat transfer problems are considered (the transient and steady-state cases), the objective function expressions are defined separately for each case. For the transient heat transfer problem, the inverse analysis is concerned with the estimation of the time-varying line heat source strength gL(t) applied at the time ti(i=1,…,r where r is the number of time steps) and the place (x†,y†) using the transient temperature measurements of a single sensor S placed at the point (Si,Sj) inside the body. In the transient inverse analysis, the aim is to minimize the difference between the estimated temperatures at the sensor place, Te(Si,Sj,ti), which is obtained from the solution of the direct transient heat conduction equation using the estimated line heat source strength gL(ti), and the measured temperatures Tm(Si,Sj,ti) over the time domain 0<t<tr as much as possible. Thus the objective function can be mathematically expressed as

min{𝒥gL at (x†,y†):=‖Te(Si,Sj,t)−Tm(Si,Sj,t)‖2:Eq.(4) in Ω,BCs and IC in Eqs.(5)−(7)}(20)

where gL=[gL(t1),gL(t2),gL(t3),…,gL(tr)]T. Therefore,

𝒥=∑i=1r[Te(Si,Sj,ti)−Tm(Si,Sj,ti)]2(21)

Sensitivity Analysis for Transient Heat Conduction

In this study, a gradient based minimization method (steepest descent method) is used to minimize the objective function value for both the transient and steady-state problems. Thus the gradient of the objective function with respect to the unknown variables is required. In the transient inverse problem, the gradient of the objective function 𝒥 defined by Eq. (21) with respect to gL(ti), i=1,…,r can be obtained by

∂𝒥∂gL(ti′)=2∑i=1r[Te(Si,Sj,ti)−Tm(Si,Sj,ti)]∂Te(Si,Sj,ti)∂gL(ti′)(22)

This is a function estimation problem and there is no prior information on the functional form of the unknown variable (here the line heat source strength gL(t)). The function estimation problems are usually solved using the adjoint problem to obtain the gradient of the objective function. However, the adjoint method comes with mathematical complexity. In addition, the adjoint equation solution itself has a computational cost comparable to the direct solution. Here, as it will be shown, using a parameter estimation approach, the unknown functional form of the line heat source strength can be calculated efficiently and accurately without solving the additional equations. In the proposed method, the elements of the sensitivity matrix are calculated during the heat conduction equation solution using some explicit expressions. Using the chain rule, we can write [22]

∂Te(Si,Sj,ti)∂gL(ti′)=(∂Te(Si,Sj,ti)∂kT|at sensor place S)(∂gL(ti′)∂kT|at source place L)(23)

This implies that two variables at two different times ti and ti′ can be related using a constant variable (here the constant thermal conductivity kT). The term ∂Te(Si,Sj,ti)∂kT can be obtained by taking derivative of Ti,jn+1 in Eq. (18) with respect to kT, as follows

∂Ten+1(Si,Sj)∂kT=ΔtρcJ2(α(TSi+1,Sjn−2TSi,Sjn+TSi−1,Sjn)

−2β14(TSi+1,Sj+1n−TSi−1,Sj+1n−TSi+1,Sj−1n+TSi−1,Sj−1n)+γ(TSi,Sj+1n−2TSi,Sjn+TSi,Sj−1n))(24)

In Eq. (23), ∂gL(ti′)∂kT can be obtained by initially solving Eq. (9) for gL(t) and then taking the derivative of the resulting expression with respect to kT, as follows

gL(ti′)=ρc∂T(ξ,η,t)∂t(4cosh2⁡(n(x−x†))cosh2⁡(n(y−y†)))n2−

kT(4cos⁡h2(n(x−x†))cosh2⁡(n(y−y†)))J2n2(α∂2T(Li,Lj,ti′)∂ξ2−2β∂2T(Li,Lj,ti′)∂ξ∂η+γ∂2T(Li,Lj,ti′)∂η2)(25)

thus

∂gL(ti′)∂kT=−(4cos⁡h2(n(x−x†))cosh2⁡(n(y−y†)))(α∂2T(Li,Lj,ti′)∂ξ2−2β∂2T(Li,Lj,ti′)∂ξ∂η+γ∂2T(Li,Lj,ti′)∂η2)J2n2(26)

∂gL(ti′)∂kT is calculated at the heat source place at L(x=x†,y=y†) which corresponds to L(Li,Lj) in the computational domain. Therefore,

cosh2⁡(n(x−x†))cosh2⁡(n(y−y†))|(x=x†,y=y†)=1(27)

And Eq. (26) becomes

∂gL(ti′)∂kT=−(4)(α∂2T(Li,Lj,ti′)∂ξ2−2β∂2T(Li,Lj,ti′)∂ξ∂η+γ∂2T(Li,Lj,ti′)∂η2)J2n2(28)

Now all the sensitivity coefficients ∂Te(Si,Sj,ti)∂gL(ti′) can be calculated during the solution of the transient heat conduction equation (direct solution). The sensitivity matrix Ja can be explicitly expressed as

JagL(t)=[∂Te(Si,Sj,t1)∂gL(t1)∂Te(Si,Sj,t1)∂gL(t2)∂Te(Si,Sj,t1)∂gL(t3)…∂Te(Si,Sj,t1)∂gL(tr)∂Te(Si,Sj,t2)∂gL(t1)∂Te(Si,Sj,t2)∂gL(t2)∂Te(Si,Sj,t2)∂gL(t3)…∂Te(Si,Sj,t2)∂gL(tr)∂Te(Si,Sj,t3)∂gL(t1)∂Te(Si,Sj,t3)∂gL(t2)∂Te(Si,Sj,t3)∂gL(t3)…∂Te(Si,Sj,t3)∂gL(tr)         ⋮∂Te(Si,Sj,tr)∂gL(t1)         ⋮∂Te(Si,Sj,tr)∂gL(t2)         ⋮∂Te(Si,Sj,tr)∂gL(t3)⋱…         ⋮∂Te(Si,Sj,tr)∂gL(tr)]r×r(29)

At any time t, the current temperature distribution is independent of a future heat source component which means that the sensitivity matrix is lower-triangular. In other words, for i′>i (the elements above the main diagonal of the sensitivity matrix), ∂Te(Si,Sj,ti)∂gL(ti′)=0. So, the sensitivity matrix will be as follows

JagL(t)=[∂Te(Si,Sj,t1)∂gL(t1)           0         0…            0∂Te(Si,Sj,t2)∂gL(t1)∂Te(Si,Sj,t2)∂gL(t2)         0…            0∂Te(Si,Sj,t3)∂gL(t1)∂Te(Si,Sj,t3)∂gL(t2)∂Te(Si,Sj,t3)∂gL(t3)…            0          ⋮∂Te(Si,Sj,tr)∂gL(t1)           ⋮∂Te(Si,Sj,tr)∂gL(t2)         ⋮∂Te(Si,Sj,tr)∂gL(t3)⋱…        ⋮∂Te(Si,Sj,tr)∂gL(tr)]r×r(30)

3.2 Objective Function for Steady-State Heat Conduction

For the steady-state case, the sensors are placed at the boundary Γ2. Therefore, the objective function is defined as

min{𝒥gL at (x†,y†):=‖TΓ2−Tm‖2:Eq.(19) in Ω,BCs in Eqs.(20) and (21)}(31)

where Tm is the vector of the measured temperature. Eq. (31) can be expressed as

𝒥=∑i=2M−1[Ti,N−T(i,N)m]2(32)

where the Ti,N and T(i,N)m, i=2,…,M−1, are the estimated and measured temperatures at the outer surface Γ2, respectively. The gradient of the objective function 𝒥 with respect to the unknown heat source strength gL is expressed as

∂𝒥∂gL=2∑i=2M−1[Ti,N−T(i,N)m]∂Ti,N∂gL(33)

Again, the sensitivity coefficients can be calculated using the chain rule as

∂Ti,N∂gL=(∂Ti,N∂kT)(∂gL∂kT|at source place L)(34)

The term ∂Ti,N∂kT can be obtained from the boundary condition at the surface Γ2, Eq. (12). By solving the boundary condition equation for Ti,N and then taking derivative of Ti,N with respect to kT, we obtain

∂Ti,N∂kT=−2h2Jγ(−4γTi,N−1+γTi,N−2−βTi+1,N+βTi−1,N+3γT∞2)(3kTγ+2h2Jγ)2(35)

The term ∂gL∂kT is the same as the term for the transient case. Thus from Eq. (28), we have

∂gL∂kT=−(4)(α∂2TLi,Lj∂ξ2−2β∂2TLi,Lj∂ξ∂η+γ∂2TLi,Lj∂η2)J2n2(36)

Thus the sensitivity coefficients for the steady-state heat conduction case will be as

∂Ti,N∂gL=−(J2n2)2h2Jγ(−4γTi,N−1+γTi,N−2−βTi+1,N+βTi−1,N+3γT∞2)(−4)(α∂2TLi,Lj∂ξ2−2β∂2TLi,Lj∂ξ∂η+γ∂2TLi,Lj∂η2)(3kTγ+2h2Jγ)2=

J3n2h2γ(−4γTi,N−1+γTi,N−2−βTi+1,N+βTi−1,N+3γT∞2)(2)(α∂2TLi,Lj∂ξ2−2β∂2TLi,Lj∂ξ∂η+γ∂2TLi,Lj∂η2)(3kTγ+2h2Jγ)2(37)

Sensitivity Analysis for Steady-State Heat Conduction

In this case, the sensitivity matrix Ja can be explicitly written as

JagL=[∂T2,N∂gL∂T3,N∂gL     ⋮∂TM−1,N∂gL](M−2)×1(38)

3.3 Minimization: Steepest Descent Method

The steepest descent minimization method used in this study is expressed as

gL(k+1)=gL(k)+β(k)d(k)(transient heat conduction)(39)

gL(k+1)=gL(k)+β(k)d(k)(steady-state heat conduction)(40)

where d(k) (and d(k)) is the direction of steepest descent and β(k) is the search step size. The negative of the gradient direction ∇𝒥(k) denotes the direction of steepest descent d(k). Thus,

d(k)=−∇𝒥(k)(transient heat conduction)(41)

d(k)=−∇𝒥(k)(steady-state heat conduction)(42)

where the gradient direction ∇𝒥 for the transient and steady-state heat conduction cases are defined in Eqs. (22) and (33), respectively. The search step size β(k)>0 is given by [1]

β(k)=−[Ja(k)d(k)]T[Te−Tm][Ja(k)d(k)]T[Ja(k)d(k)](transient heat conduction)(43)

β(k)=−[Ja(k)d(k)]T[Te−Tm][Ja(k)d(k)]T[Ja(k)d(k)](steady-state heat conduction)(44)

Minimization Algorithm

The inverse analysis steps proposed in this study to estimate the constant and time-varying line heat source strengths can be summarized as below:

1.   The temperature measurements are taken at the sensor place SSi,Sj and the time ti (i=1,…,r), Tm(Si,Sj,ti), for the transient case and the temperature measurements are taken at the outer surface Γ2, T(i,N)m, for the steady-state case.

2.   The direct problem is solved to obtain the temperature values at the sensor place and the time ti (i=1,…,r), Te(Si,Sj,ti), through solving Eq. (18) for the transient case and to obtain the temperature values at the outer surface Γ2, Ti,N, for the steady-state case. During the direct solution, the sensitivity coefficients (Eq. (23) for the transient case and Eq. (37) for the steady-state case) are computed. For the transient case, after the computations of the sensitivity coefficients (during the direct solution), the sensitivity matrix (Eq. (30)) can be computed as follows

doi=1,r     do j=1,r           if(i.LT.j)then           Ja(i,j)=0.0           else           Ja(i,j)=TkT(i,1)gLkT(j,1)           endif     enddoenddo

3. The objective function (𝒥(k)) is computed using Eqs. (21) and (32) for the transient and steady-state cases, respectively.

4. The specified stopping criterion is examined for stopping or continuation of the minimization process.

5. The gradient direction ∇𝒥(k), the direction of descent d(k) (and d(k)), and the search step size β(k) are computed.

6. gL (transient case) and gL (steady-state case) are updated from Eqs. (39) and (40), respectively. Set the next iteration (k=k+1) and return to the step 2.

3.4 Stopping Criterion

The minimization process to estimate the unknown variables can be stopped when the objective function is minimized sufficiently and stable and satisfactory results are obtained. In this study the minimization process is terminated when the objective function value 𝒥(k) becomes less than a very small number (ε=10−5 or ε=10−6 depending on the test case). For the case of the temperature measurements with errors, the discrepancy principle is used. In the discrepancy principle, the minimization process can be stopped when the difference between the estimated and measured temperatures is of the order of magnitude of the measurement errors. This can be mathematically expressed by

|Te−Tm|≈σ(45)

where σ is the standard deviation of the measurement errors and is assumed constant. The stopping criteria can then be obtained by substituting Eq. (45) into the objective function expression for each of the transient (Eq. (21)) and steady-state (Eq. (32)) case, as follows

εtransient=rσ2(46)

εsteady-state=(M−2)σ2(47)

The measured temperatures containing random errors Tmeas are generated by adding an error term ωσ to the exact temperatures Texact to give

Tmeas=Texact+ωσ(48)

where ω is a random variable with normal (Gaussian) distribution, zero mean, and unitary standard deviation. Assuming 99% confidence for the measured temperature, ω lies in the range −2.576≤ω≤2.576 and it is randomly generated by using MATLAB.

4  Results

For each of the inverse transient and steady-state heat conduction problems, two test cases are presented. The data used in the test cases are so that they can reflect the accuracy, efficiency, and robustness of the proposed inverse analysis. At each test case, initially the direct heat conduction problem is solved with prescribed values of the line heat source strength to compute the temperature value at the sensor place. Then the computed temperatures are used as simulated measured ones to recover the prescribed values of the line heat source strength. For the function estimation problems (transient case), two different forms of timewise variation for the heat source are considered.

Before proceeding further, the steady-state and transient heat conduction equations are validated with the results from the commercial finite element software COMSOL. We assume that the heat conducting body is made of stainless steel (type 304) and the data used in the transient and steady-state test cases are listed in Tables 1 and 2, respectively. It should be added that all computations in this study are carried out using a Fortran compiler.

The temperature distribution in the irregular heat conducting body for the transient heat transfer problem using the explicit code used in this study and the software COMSOL is shown in Fig. 2. Moreover, the temperature distribution along a radial line passing over the heat source computed by our code and COMSOL (Fig. 3) are compared (Fig. 4). And the temperature distribution in the irregular heat conducting body for the steady-state heat transfer problem using the explicit code and COMSOL is shown in Fig. 5. It can be seen that the there is a good agreement between the results obtained from our code and the ones obtained from the software COMSOL confirming the accuracy of our code. Furthermore, the shape and the location of the function δn(x−x†)δn(y−y†)=n24cosh2⁡(n(x−x†))cosh2⁡(n(y−y†)) for n=80 and x†=0.13504165411 m and y†=0.303559482098 m corresponding to the node (Li,Lj)=(60,60) (Table 2) is shown in Fig. 6. It can be seen that the maximum value of the function is at x=x† and y=y†. Noting that in this study n=80 and cosh2⁡(n(x−x†))cosh2⁡(n(y−y†))|(x=x†,y=y†)=1, we obtain

δ80(x†−x†)δ80(y†−y†)=8024=1600(49)

Figure 2: Validation of the transient heat conduction solver using the finite element software COMSOL. The result obtained by using our code (a) and the result obtained by using COMSOL (b)

Figure 3: Temperature along a radial line passing over heat source obtained by software COMSOL. (a) the place of line (b,c) the temperature distribution

Figure 4: Comparison of temperature distribution along a radial line passing over heat source obtained by our code and software COMSOL

Figure 5: Validation of the steady-state heat conduction solver using the finite element software COMSOL. The result obtained by using our code (a) and the result obtained by using COMSOL (b)

Figure 6: The shape and location of the function δn(x−x†)δn(y−y†) for n=80 and x†=0.13504165411 m and y†=0.303559482098 m

4.1 Test Cases 1 & 2 (Transient)

As mentioned previously, initially the direct transient heat conduction problem is solved using prescribed values of the line heat source strength (gL(t)desired) to compute the temperature value at the sensor place SSi,Sj and time ti. Then the computed temperatures at the sensor place and time ti are used as simulated measured ones to recover the prescribed values of the line heat source strength using an initial guess (gL(t)initial). The different forms of time wise variation for the heat source strengths for the transient test cases are given in Eqs. (50)–(53). And the data used in Test cases 1 & 2 are the same as the ones used for the validation test case (Table 1) except for the parameters listed in Table 3. The physical domain and the grid size (60×60) used for all four test cases are shown in Fig. 7. The estimation results for Test cases 1 & 2 including a comparison of the initial, final, and desired line heat source strengths (functions) as well as the objective function history are shown in Figs. 8 and 9, respectively. To investigate the effect of the measurement error on the inverse problem solution, two different measurement error levels are considered in the inverse analysis for Test case 1, including σ=0.01 (Fig. 8c,d), and σ=0.1 (Fig. 8e,f). It can be seen that the relatively good results are also obtained when there is large measurement error. In Test case 2, a more complicated initial guess is considered for the heat source strength to investigate the robustness of the proposed inverse analysis. As can be seen from the results for both test cases, there is an excellent agreement between the desired and final line heat source strengths, showing that the inverse analysis proposed in this study is very accurate and robust. As stated previously, all elements of the sensitivity matrix are computed during the heat equation solver without the need to solve the sensitivity and adjoint problems. So the computational cost for the inverse analysis at each iteration is approximately the same as the direct problem solution one. Despite the large number of the iterations in Test cases 1 & 2, the computation time is only a few hours on a PC with Intel Core i5 and 6G RAM.

Figure 7: The physical domain and the grid size of 60×60 used for all test cases

Figure 8: Test case 1 results. Estimation of the time-varying line heat source strength gL(t) using an initial guess gL(t)initial,Test case 1=22300(Wm) and objective function versus iteration number for cases of no measurement error (a,b) and measurement errors of σ=0.01 (c,d) and σ=0.1 (e,f). The final time and time step are tr=300(s) and Δt=0.2(s), respectively