515
Views
4
CrossRef citations to date
0
Altmetric
Articles

Identification of the timewise thermal conductivity in a 2D heat equation from local heat flux conditions

ORCID Icon
Pages 903-919 | Received 17 Dec 2019, Accepted 17 Aug 2020, Published online: 03 Sep 2020

ABSTRACT

The aim of this paper is to identify numerically the timewise thermal conductivity coefficients in the two-dimensional heat equation in a rectangular domain using initial and Dirichlet boundary conditions and the local heat flux as over-specification conditions. The measurement data represented by the local heat flux is shown to ensure the unique solvability of the inverse problem solution. The two-dimensional inverse problem is discretized using an alternating direction explicit method. The resulting constrained optimization problem is minimized iteratively by employing the MATLAB subroutine. Exact and noisy input data are inverted numerically. The  root mean square error values for various noise levels p for both smooth and non-smooth continuous timewise thermal conductivity coefficients Examples are compared. Numerical results are presented and discussed in order to illustrate the performance of the inversion for thermal conductivity components identification. This study will be significant to researchers working on computational and mathematical methods for solving inverse coefficient identification problems with applications in heat transfer and porous media.

2010 Mathematics Subject Classifications:

1. Introduction

The identification problems of the thermal conductivity/diffusivity coefficients in the parabolic partial differential equation from the local heat flux and/or non-local boundary measurements have received important attention from many researchers in the past, see [Citation1–6] to mention only a few. The reconstruction of thermal conductivity coefficient along with the temperature in the heat equation from the non-local linear combination of heat flux and with local heat flux measurements has been investigated in [Citation7] and [Citation8], respectively. The determination of the orthotropic thermal conductivity coefficients and the temperature in the two-dimensional parabolic heat equation from the non-local heat flux over-specification conditions has been investigated in [Citation9].

In recent years, the inverse heat conduction problems have received widespread attention in engineering applications [Citation10,Citation11]. Mera et al. [Citation12] estimated the unknown conductivity coefficients and the unknown boundary data for the two-dimensional heat conduction problem in an anisotropic medium. Reddy and Dulikravich [Citation13] determined simultaneously the spatially varying thermal conductivity and specific heat using boundary temperature measurements. Reddy et al. [Citation14] estimated a single spatially varying material property (thermal conductivity) based on the solutions of a steady-state version of the heat conduction equation.

In a recent paper, [Citation15], the authors have investigated the recovery of the time-dependent thermal conductivity coefficients a1(t) and a2(t) of an orthotropic rectangular conductor along with the temperature u(x,y,t) in a two-dimensional problem given by the parabolic heat equation ut=a1(t)uxx+a2(t)uyy+f(x,y,t), from non-local over-specification conditions a1(t)(ν11(t)ux(0,y0,t)+ν12(t)ux(h,y0,t))=ϰ1(t), and a2(t)(ν21(t)uy(x0,0,t)+ν22(t)uy(x0,l,t))=ϰ2(t), where x0, y0 are fixed values from the intervals (0,h) and (0,l), and h, l are given positive quantities. When ν11=ν21=1 and ν21=ν22=0, the above non-local over-specification conditions become (Equation5) and (Equation6), respectively. In this paper, we generalize the above non-local over-specifications to the local heat flux measurements (Equation5) and (Equation6) for a more general two-dimensional parabolic heat equation, as described in the next section. The inverse problem investigated in this paper has already been proved to be locally uniquely solvable, but no numerical identification has been attempted so far, and it is the aim of the current research to undertake the numerical solution of this problem. In this work, the novelty consists in the development of a numerical optimization method for solving this nonlinear inverse coefficient problem for the two-dimensional heat equation. Numerically, the implementation is realized using the MATLAB subroutine lsqnonlin.

The structure of the paper is as follows. The mathematical formulation of the two-dimensional heat equation is formulated in Section 2. The numerical solution of the direct problem based on an alternating direction explicit method is presented in Section 3. The numerical solution of the inverse problem based on the minimization of the nonlinear least-squares objective function is introduced in Section 4. Numerical results are presented and discussed in Section 5. Finally, conclusions are given in Section 6.

2. Statement of the 2D heat equation problem

We consider an inverse problem of identifing the timewise thermal conductivity coefficients a1(t)>0 and a2(t)>0 in the two-dimensional heat equation [Citation16], (1) ut(x,y,t)=a1(t)b1(x)2ux2(x,y,t)+a2(t)b2(y)2uy2(x,y,t)+c1(x,y,t)ux(x,y,t)+c2(x,y,t)uy(x,y,t)+c3(x,y,t)u(x,y,t)+f(x,y,t),(x,y,t)QT,(1) where b1>0, b2>0 are known space variations, c1, c2 are known velocity representing convective flow, c3 is a known absorption coefficient, f is a known heat source, u is a unknown temperature in the rectangular domain QT={(x,y,t):0<x<h,0<y<l,0<t<T}, subject to the initial condition (2) u(x,y,0)=φ(x,y),(x,y)[0,h]×[0,l],(2) the Dirichlet boundary conditions (3) u(0,y,t)=μ1(y,t),u(h,y,t)=μ2(y,t),(y,t)[0,l]×[0,T],(3) (4) u(x,0,t)=μ3(x,t),u(x,l,t)=μ4(x,t),(x,t)[0,h]×[0,T],(4) and the local heat flux over-specification conditions (5) a1(t)ux(0,Y0,t)=ν1(t),t[0,T],(5) (6) a2(t)uy(X0,0,t)=ν2(t),t[0,T],(6) where X0, Y0 are fixed values from the intervals (0,h) and (0,l), where h, l are given positive numbers, and φ(x,y), μ1(y,t),μ2(y,t),μ3(x,t),μ4(x,t) are given functions satisfying compatibility conditions. The geometry of the inverse problem under investigation is shown in Figure .

Figure 1. Geometry of the inverse problem under investigation.

Figure 1. Geometry of the inverse problem under investigation.

The local existence and uniqueness of solution of the inverse problem (Equation1)–(Equation6) were established in [Citation16] and read as stated in the following two theorems.

Theorem 2.1

Suppose that the following assumptions are satisfied:

(A1)

φ(x,y)C2(D¯),whereD:={(x,y):0<x<h,0<y<l},μi(y,t)C2,1([0,l]×(0,T]),i=1,2,μi(x,t)C2,1([0,h]×(0,T]),i=3,4,νi(t)C[0,T],i=1,2,f(x,y,t)C1,1,0(QT¯),b1(x)C1[0,h],b2(y)C1[0,l],ci(x,y,t)C1,1,0(QT¯),i=1,2,3;

(A2)

φ(x,y)0,φx(x,y)>0,φy(x,y)>0,(x,y)D¯,μ1(y,t)0,μ2(y,t)0,μ1y(y,t)>0,μ2y(y,t)>0,μ1yy(y,t)0,μ2yy(y,t)0,(y,t)[0,l]×[0,T],μ3(t)0,μ4(t)0,μ3x(t)>0,μ4x(t)>0,μ3xx(x,t)0,μ4xx(x,t)0,(x,t)[0,h]×[0,T],f(x,y,t)0,fx(x,y,t)>0,fy(x,y,t)>0,c1(0,y,t)<0,c1(h,y,t)>0,c2(x,0,t)<0,c2(x,l,t)>0,c1y(x,y,t)0,c2x(x,y,t)0,c3x(x,y,t)0,c3y(x,y,t)0,(x,y,t)(QT¯),b1(x)>0,x[0,h],b2(y)>0,y[0,l],νi(t)0,i=1,2,t[0,T];

(A3)

μ1t(y,t)c2(0,y,t)μ1y(y,t)c3(0,y,t)μ1(y,t)f(0,y,t)<0,μ2t(y,t)c2(h,y,t)μ2y(y,t)c3(h,y,t)μ2(y,t)f(h,y,t)>0,μ3t(x,t)c1(x,0,t)μ3x(x,t)c3(x,0,t)μ3(x,t)f(x,0,t)<0,μ4t(x,t)c1(x,l,t)μ4x(x,t)c3(x,l,t)μ4(x,t)f(x,l,t)>0;

(A4)

μ1(y,0)=φ(0,y),μ2(y,0)=φ(h,y),μ3(x,0)=φ(x,0),μ4(x,0)=φ(x,l),μ1(0,t)=μ3(0,t),μ1(l,t)=μ4(0,t),μ2(0,t)=μ3(h,t),μ2(l,t)=μ4(h,t).

Then, there exists T0(0,T], which is determined by the input data, such that the problem (Equation1)–(Equation6) has a solution (a1(t),a2(t),u(x,y,t))C1[0,T0]×C1[0,T0]×C2,1(QT0¯), with ai(t)>0 for t[0,T0], i = 1, 2.

Theorem 2.2

Suppose that the condition νi(t)>0, i = 1,2 for t[0,T] is satisfied. Then, the inverse problem (Equation1)–(Equation6) cannot have more than one solution in the class (a1(t),a2(t),u(x,y,t))C1[0,T]×C1[0,T]×C2,1(QT¯), with ai(t)>0 for t[0,T], i = 1, 2.

3. Numerical solution of direct problem

In this section, we consider the direct initial boundary value problem (Equation1)–(Equation4), where a1(t), a2(t), b1(x), b2(y), ci(x,y,t), i=1,3¯, f(x,y,t), φ(x,y), μi(y,t), i = 1, 2, and μi(x,t), i = 3, 4 are known and the solution u(x,y,t) is to be determined together with the quantities of interest νi(t), i=1,2. We subdivide the solution domain QT into M1, M2 and N subintervals of equal step lengths Δx, Δy and Δt, where Δx=h/M1, Δy=l/M2, and Δt=T/N, respectively. At the node (i,j,n), we denote ui,jn:=u(xi,yj,tn), where xi=iΔx, yj=jΔy, tn=nΔt, a1n:=a1(tn), a2n:=a2(tn), c1i,jn:=c1(xi,yj,tn), c2i,jn:=c2(xi,yj,tn), c3i,jn:=c3(xi,yj,tn), b1i:=b1(xi), b2j:=b2(yj) and fi,jn:=f(xi,yj,tn) for i=0,M1¯, j=0,M2¯, n=0,N¯.

Alternating direction explicit (ADE) method makes use of two approximations that are implemented for computations proceeding in alternating directions, e.g. from left to right and from right to left, with each approximation being explicit in its respective direction of computation [Citation17–19]. The most important advantage of numerical simulation is the low cost and high computational speed. Furthermore, the numerical scheme possesses the advantage of the implicit methods, i.e. no severe limitation in the size of the time increment. Also, the ADE has another advantage, namely, that is unconditionally stable. For solving the direct problem (Equation1)–(Equation4) the ADE is as follows.

Let u~i,jn and v~i,jn be the solutions of the following equations which are multilevel finite difference discretization of Equation (Equation1): (7) u~i,jn+1u~i,jnΔt=a1nb1i(u~i+1,jnu~i,jnu~i,jn+1+u~i1,jn+1(Δx)2)+a2nb2j(u~i,j+1nu~i,jnu~i,jn+1+u~i,j1n+1(Δy)2)+c1i,jn(u~i+1,jnu~i1,jn+12(Δx))+c2i,jn(u~i,j+1nu~i,j1n+12(Δy))+c3i,jn(u~i,jn+1+u~i,jk2)+12(fi,jn+1+fi,jn),i=1,M11¯,j=1,M21¯,n=0,N1¯,(7) (8) v~i,jn+1v~i,jnΔt=a1nb1i(v~i+1,jn+1v~i,jn+1v~i,jn+v~i1,jn(Δx)2)+a2nb2j(v~i,j+1n+1v~i,jn+1v~i,jn+v~i,j1n(Δy)2)+c1i,jn(v~i+1,jn+1v~i1,jn2(Δx))+c2i,jn(v~i,j+1n+1v~i,j1n2(Δy))+c3i,jn(v~i,jn+v~i,jn+12)+12(fi,jn+1+fi,jn),i=M11,1¯,j=M21,1¯,n=0,N1¯.(8)

Furthermore, let the u~i,jn and v~i,jn also satisfy the initial and boundary conditions (Equation2)–(Equation4), namely, (9) u~i,j0=v~i,j0=φ(xi,yj),i=0,M1¯,j=0,M2¯,(9) (10) u~0,jn=v~0,jn=μ1(yj,tn),u~M1,jn=v~M1,jn=μ2(yj,tn),j=0,M2¯,n=1,N¯,(10) (11) u~i,0n=v~i,0n=μ3(xi,tn),u~i,M2n=v~i,M2n=μ4(xi,tn),i=0,M1¯,n=1,N¯.(11) Rearranging the terms in (Equation7) and (Equation8), we obtain the explicit calculations of u~i,jn+1 and v~i,jn+1 as follows: (12) u~i,jn+1=Ai,jnu~i,jn+Bi,jn(u~i+1,jn+u~i1,jn+1)+Ci,jk(u~i,j+1n+u~i,j1n+1)+Di,jn(u~i+1,jnu~i1,jn+1)+Ei,jn(u~i,j+1nu~i,j1n+1)+Gi,j,i=1,M11¯,j=1,M21¯,n=0,N1¯,(12) (13) v~i,jn+1=Ai,jnv~i,jn+Bi,jn(v~i+1,jn+1+v~i1,jn)+Ci,jn(v~i,j+1n+1+v~i,j1n)+Di,jn(v~i+1,jn+1v~i1,jn)+Ei,jn(v~i,j+1n+1v~i,j1n)+Gi,j,i=M11,1¯,j=M21,1¯,n=0,N1¯,(13) where (14) Ai,jn=1λi,jn1+λi,jn,Bi,jn=(Δt)a1nb1i(Δx)2(1+λi,jn),Ci,jn=(Δt)a2nb2j(Δy)2(1+λi,jn),Di,jn=(Δt)c1i,jn2(Δx)(1+λi,jn),Ei,jn=(Δt)c2i,jn2(Δy)(1+λi,jn),Gi,j=Δt2(1+λi,jn)(fi,jn+1+fi,jn),λi,jn=Δt(a1nb1i(Δx)2+a2nb2j(Δy)2c3i,jn2).(14) From (Equation12) and (Equation9)–(Equation11) for u~, u~i,jn+1 can be computed explicitly. In this case, calculations proceed from the grid point close to the boundaries x = 0 and y = 0, as i, j are increasing. The needed values such as u~i1,jn+1, u~i,jn and u~i+1,jn will be known from initial and boundary conditions (Equation9)–(Equation11). Similarly, v~i,jn+1 can be calculated explicitly from (Equation13) and (Equation9)–(Equation11) for v~, beginning at the boundaries x = 1 and y = 1 and marching in a sequence of decreasing i and j, i.e. i=M11,M12,,1, j=M21,M22,,1. These values are then substituted into the simple arithmetic mean approximation (15) ui,jn+1=u~i,jn+1+v~i,jn+12(15) to obtain the solution ui,jn+1.

The local heat flux over-specification conditions (Equation5) and (Equation6) can be calculated using finite difference method approximation as follows: (16) ν1(tn)=a1n(4u(1,Y0,tn)u(2,Y0,tn)3u(0,Y0,tn)2Δx),n=1,N¯,(16) (17) ν2(tn)=a2n(4u(X0,1,tn)u(X0,2,tn)3u(X0,0,tn)2Δy),n=1,N¯.(17)

4. Numerical solution of inverse problem

In this section, our goal is to obtain simultaneously stable identifications for the timewise thermal conductivity coefficients a1(t) and a2(t), satisfying Equations (Equation1)–(Equation6). We reduce the inverse problem to a nonlinear minimization of the least-squares objective function (18) F(a1,a2)=a1(t)ux(0,Y0,t)ν1(t)L2[0,T]2+a2(t)uy(X0,0,t)ν2(t)L2[0,T]2,(18) or, in discretizations form (19) F(a1,a2)=n=1N[a1nux(0,Y0,tn)ν1(tn)]2+n=1N[a2nuy(X0,0,tn)ν2(tn)]2,(19) where u(x,y,t) solves numerically using the ADE method the direct problem (Equation1)–(Equation4) for given a1(t) and a2(t), respectively. It is worth mentioning that in (Equation19) at the first time step, i.e. n = 0, the derivative ux(0,Y0) and uy(X0,0) are obtained from the initial condition (Equation2), via (Equation16) and (Equation17), as (20) ux(0,Y0)=4φ1,Y0φ2,Y03φ0,Y02(Δx),(20) (21) uy(X0,0)=4φX0,1φX0,23φX0,02(Δy),(21) where φi,j=φ(xi,yj) for i=0,M1¯, j=0,M2¯. Also, one can remark that at initial time t=0 the values a1(0) and a2(0) can be obtained from the local heat flux over-specification conditions (Equation5) and (Equation6) as (22) a1(0)=ν1(0)φx(0,Y0),a2(0)=ν2(0)φy(X0,0).(22) The minimization of the objective function (Equation19) is performed using the MATLAB subroutine lsqnonlin, [Citation20]. This routine attempts to find the minimum of a sum of squares by starting from the initial guesses. Simple bounds on the variables are allowed and the explicit calculation (analytical or numerical) of the gradient is not required to be supplied by the user. Furthermore, within lsqnonlin, we use the Trust-Region-Reflective (TRR) algorithm [Citation21], which is based on the interior-reflective Newton method. Each iteration involves a large linear system of equations whose solution, based on a preconditioned conjugate gradient method, allows a regular and sufficiently smooth decrease of the objective functional (19).

The inverse problem (Equation1)–(Equation6) is solved subject to both exact and noisy measurements (Equation5) and (Equation6). The noisy data is numerically simulated as (23) ν1ϵ1(tn)=ν1(tn)+ϵ1n,ν2ϵ2(tn)=ν2(tn)+ϵ2n,n=0,N¯,(23) where ϵ1n and ϵ2n are random variables generated from a Gaussian normal distribution with mean zero and standard deviations σ1 and σ2, respectively, given by (24) σ1=p×maxt[0,T]|ν1(t)|,σ2=p×maxt[0,T]|ν2(t)|,(24) where p represents the percentage of noise. We use the MATLAB function normrnd to generate the random variables ϵ1_=(ϵ1n)n=1,N¯ and ϵ2_=(ϵ2n)n=1,N¯ as follows: (25) ϵ1_=normrnd(0,σ1,N),ϵ2_=normrnd(0,σ2,N).(25) In the case of noisy data (Equation23), we replace ν1(tn) and ν2(tn) by ν1ϵ1(tn) and ν2ϵ1(tn), respectively, in (Equation19).

5. Numerical results and discussion

In this section, we present numerical results for the thermal conductivity a1(t) and a2(t) together with the temperature u(x,y,t), in the case of exact and noisy data (Equation23)–(Equation25). To quantify the accuracy of the approximate solution, we employ the root mean square errors (rmse) defined by (26) rmse(a1)=[TNn=1N(a1Numerical(tn)a1Exact(tn))2]1/2,(26) (27) rmse(a2)=[TNn=1N(a2Numerical(tn)a2Exact(tn))2]1/2.(27) For simplicity, we take h = l = T = 1. The lower bounds and upper bounds for the coefficients a1(t) and a2(t) are 1010 for (LB) and 102 for (UB). These bounds allow a wide search range for the unknown positive physical components a1(t) and a2(t).

From the physical point of view, the phenomenon of thermal conductivity/diffusivity is a transfer of kinetic energy. For instance, in metallic crystals the heat energy is transferred by two types of carries, conductivity electronics and oscillations of the crystalline lattice. Furthermore, physical situations in which the thermal conductivity a1(t) and a2(t) depends on time (and are unknown) occur in damage and radioactive decay applications [Citation22].

5.1. Example 1

Consider the inverse problem (Equation1)–(Equation6) with unknown coefficients a1(t) and a2(t), with the input data: (28) b1(x)=1+x,b2(y)=1+y,c1(x,y,t)=3x+y100t1,X0=0.5,Y0=0.5,c2(x,y,t)=x100+3yt1,c3(x,y,t)=x+y+t10,φ(x,y)=1+x+3y+3xy,μ1(y,t)=1+3y+ty310,μ2(y,t)=2+(6+t)y+ty310,μ3(x,t)=1+x,μ4(x,t)=4(1+x)+t(110+x2),f(x,y,t)=2t(1+t)(1+x)y35t(1+2t)y(1+y)(1t+3x+y100)(1+(3+2tx)y)+y(x2+y210)(1t+x100+3y)(3+3x+tx2+3ty210)(10+t+x+y)(1+x+3y+3xy+tx2y+ty310),(28) (29) ν1(t)=52(1+t),ν2(t)=14(18+t)(1+2t).(29) One can notice that the conditions of Theorems 2.1 and 2.2 are satisfied, and therefore, the uniqueness of the solution is guaranteed. In fact, it can easily be checked by direct substitution that the analytical solution is given by (30) u(x,y,t)=1+x+3y+3xy+tx2y+110ty3,(x,y,t)Q¯T,(30) (31) a1(t)=1+t,a2(t)=1+2t,t[0,1].(31) First, we assess the accuracy of the direct problem (Equation1)–(Equation4) with the input data (Equation28) when a1(t) and a2(t) are known and given by (Equation31). The numerical results for the interior temperature u(x,y,t) have been obtained in excellent agreement with the analytical solution (Equation30) and therefore they are not presented. Apart from the interior temperature, other output of interest is the data (Equation5) and (Equation6), which analytically is given by (Equation29). Figure  shows that the analytical and numerical solutions for this quantity obtained with M1=M2=10 and with various numbers of time steps N{20,40,80} are in very good agreement. Also, the rmse defined by (32) rmse(ν1)=[1Nn=1N(ν1Numerical(tn)ν1Exact(tn))2]1/2,(32) (33) rmse(ν2)=[1Nn=1N(ν2Numerical(tn)ν2Exact(tn))2]1/2,(33) indicated in Table , show more clearly the convergence of the numerical ADE solution to the analytical solution (Equation29).

Figure 2. The exact (Equation29) and numerical solutions for (a) ν1(t) and (b) ν2(t), with M1=M2=10 and with various numbers of time steps N{20,40,80}, for direct problem.

Figure 2. The exact (Equation29(29) ν1(t)=52(1+t),ν2(t)=14(18+t)(1+2t).(29) ) and numerical solutions for (a) ν1(t) and (b) ν2(t), with M1=M2=10 and with various numbers of time steps N∈{20,40,80}, for direct problem.

Table 1. The rmse values ((Equation32) and (Equation33)) for ν1(t) and ν2(t), with M1=M2=10 and with various N{20,40,80}, for direct problem.

Next we investigate the inverse problem. We fix M1=M2=10 and N = 40 and we start the investigation for reconstructing the unknown timewise thermal conductivity coefficients a1(t) and a2(t) and the temperature u(x,y,t) for p{0,1%,3%,5%} noise in the measured data (Equation23). We take the initial guesses for the vectors a1 and a2 as follows: (34) a10(tn)=a1(0)=1,a20(tn)=a2(0)=1,n=1,N¯.(34) Note that the values of a1(0) and a2(0) are available from (Equation22). The objective function (Equation19), as a function of the number of iterations, is plotted in Figure . From this figure, it can be seen that a fast monotonic decreasing convergence is achieved in about 5–9 iterations to reach a very low prescribed tolerance of O(1026) in about 6 to 8 minutes CPU time. The corresponding numerical results for a1(t) and a2(t) are shown in Figure  and summarized in Table . From this figure and table, it can be seen that as the percentage of noise p decreases from 5% to 3% to 1% and then to zero the numerically obtained results becomes more stable and accurate. For more, information about the number of iterations, the values of the objective function (19) at final iteration, the rmse values ((26) and (27)) and the computational time, see Table . Figure  compares the analytical (Equation30) and numerical solutions for u(x,y,t) showing good agreement and stability. No regularization was found necessary to penalize the nonlinear least-squares objective functional (19) for amounts of noise up to p=5%, hinting towards the conclusion that the inverse problem under investigation is not severely ill-posed. Nevertheless, for higher amounts of noise the instability in retrieving the coefficients a1(t) and a2(t) will become apparent and regularization would need to be employed.

Figure 3. The objective function (Equation19), as a function of the number of iterations, for Example 1 with p{0,1%,3%,5%} noise.

Figure 3. The objective function (Equation19(19) F(a1,a2)=∑n=1N[a1nux(0,Y0,tn)−ν1(tn)]2+∑n=1N[a2nuy(X0,0,tn)−ν2(tn)]2,(19) ), as a function of the number of iterations, for Example 1 with p∈{0,1%,3%,5%} noise.

Figure 4. The exact (Equation31) and numerical solutions for: (a) a1(t) and (b) a2(t), for Example 1 with p{0,1%,3%,5%} noise.

Figure 4. The exact (Equation31(31) a1(t)=1+t,a2(t)=1+2t,t∈[0,1].(31) ) and numerical solutions for: (a) a1(t) and (b) a2(t), for Example 1 with p∈{0,1%,3%,5%} noise.

Figure 5. The analytical (Equation30) and approximate solutions for the temperature u(x,y,1), for Example 1 with (a) no noise, (b) p=1% noise, (c) p=3% noise, and (d) p=5% noise. The absolute error between them is also included.

Figure 5. The analytical (Equation30(30) u(x,y,t)=1+x+3y+3xy+tx2y+110ty3,(x,y,t)∈Q¯T,(30) ) and approximate solutions for the temperature u(x,y,1), for Example 1 with (a) no noise, (b) p=1% noise, (c) p=3% noise, and (d) p=5% noise. The absolute error between them is also included.

Table 2. The number of iterations, the values of the objective function (19) at final iteration, the rmse values ((26) and (27)) and the computational time, for p{0,1,3,5}%, for Examples 1 and 2.

5.2. Example 2

The previous example has recovered the smooth timewise thermal conductivity coefficients a1(t) and a2(t) as given in (Equation31). In this example, we assess the numerical ADE scheme for reconstructing a non-smooth test: (35) a1(t)=1+cos2(2πt),a2(t)=1+cos2(3πt),t[0,1].(35) We also have the same analytical solution for the temperature u(x,y,t) given by (Equation30). The rest of the input data are the same as in Example 1 and (36) ν1(t)=52(1+cos2(2πt)),ν2(t)=18(18+t)(3+cos(6πt)),f(x,y,t)=(1t+3x+y100)(1+(3+2tx)y)+y(x2+y210)(1t+x100+3y)(3+3x+tx2+3ty210)(10+t+x+y)(1+x+3y+3xy+tx2y+ty310)t(1+x)y(3+cos(4πt))310ty(1+y)(3+cos(6πt)).(36) Also, with this data, the conditions of Theorems 2.1 and 2.2 are fulfilled, hence the uniqueness of the solution is also guaranteed. The initial guesses for the vectors a1 and a2 for this example has been taken as (37) a10(tn)=a1(0)=2,a20(tn)=a2(0)=2,n=1,N¯.(37) Note that the values of a1(0) and a2(0) are also available from (Equation22). We investigate the inverse problem as we did in Example 1. We take M1=M2=10 and N = 40 and reconstructing the unknown thermal conductivity coefficients a1(t) and a2(t) and the temperature u(x,y,t) for exact and noisy measured input data (Equation23), for various amounts of noise p = 0 (no noise) and p{1%,3%,5%} noisy data (Equation24). The objective function (Equation19), as a function of the number of iterations, is plotted in Figure . From this figure one can notice that a rapid monotonically decreasing convergence is achieved in about 6–10 iterations, with the objective function reaching a very small value of O(1026). The numerical results for the timewise thermal conductivity coefficients a1(t) and a2(t) are depicted in Figure  and summarized in Table . From this figure and table, it can be seen that as the percentage of noise p decreases from 5% to 3% to 1% and then to zero the numerically obtained results becomes more stable and accurate.

Figure 6. The objective function (Equation19), as a function of the number of iterations, for Example 2 with p{0,1%,3%,5%} noise.

Figure 6. The objective function (Equation19(19) F(a1,a2)=∑n=1N[a1nux(0,Y0,tn)−ν1(tn)]2+∑n=1N[a2nuy(X0,0,tn)−ν2(tn)]2,(19) ), as a function of the number of iterations, for Example 2 with p∈{0,1%,3%,5%} noise.

Figure 7. The exact (Equation35) and numerical solutions for: (a) a1(t) and (b) a2(t), for Example 2 with p{0,1%,3%,5%} noise.

Figure 7. The exact (Equation35(35) a1(t)=1+cos2⁡(2πt),a2(t)=1+cos2⁡(3πt),t∈[0,1].(35) ) and numerical solutions for: (a) a1(t) and (b) a2(t), for Example 2 with p∈{0,1%,3%,5%} noise.

Finally, Figure  shows the analytical (Equation35) and numerical solutions for the temperature u(x,y,t) for various amounts of noise p{0,1,3,5}%. From this figure it can be seen that the numerical solution is stable and furthermore, its accuracy improves as the noise level p decreases. The same conclusions can be drawn about the stable reconstructions for the unknown thermal conductivity coefficients.

Figure 8. The analytical (Equation30) and approximate solutions for the temperature u(x,y,1), for Example 2 with (a) no noise, (b) p=1% noise, (c) p=3% noise, and (d) p=5% noise. The absolute error between them is also included.

Figure 8. The analytical (Equation30(30) u(x,y,t)=1+x+3y+3xy+tx2y+110ty3,(x,y,t)∈Q¯T,(30) ) and approximate solutions for the temperature u(x,y,1), for Example 2 with (a) no noise, (b) p=1% noise, (c) p=3% noise, and (d) p=5% noise. The absolute error between them is also included.

6. Conclusions

In this paper, the inverse problem involving simultaneous identification of the timewise thermal conductivity coefficients a1(t) and a2(t) and the temperature u(x,y,t) in the two-dimensional parabolic heat equation (1) in a rectangular domain from the local heat flux over-specification conditions (5) and (6) has been numerically investigated. The direct solver based on the ADE scheme is employed. The inverse problem approach based on a nonlinear least-squares minimization problem using the MATLAB optimization toolbox subroutine lsqnonlin is developed. The rmse values for various noise levels p for both smooth and non-smooth continuous timewise thermal diffusivity coefficients Examples were compared. Numerical results presented and discussed for both exact and noisy data show that accurate and stable solutions have been obtained. No regularization has been found necessary indicating that the inverse problem is not severely ill-posed. In principle, the numerical analysis of this paper can be extended to three-dimensional problems; however, this non-trivial investigation is deferred to future work.

Nomenclature

a1=

thermal conductivity (W m1 K1 )

a2=

thermal conductivity (W m1 K1 )

c1=

velocity (m s1 )

c2=

velocity (m s1 )

c3=

absorbance (AU)

f=

heat source/force (J)

u=

temperature (K or C)

ν1=

heat flux at x = 0 (W m2 )

ν2=

heat flux at y = 0 (W m2 )

t=

time variable (s)

x, y=

space variables (m)

xi,yj=

space nodes (–)

tn=

time node (–)

ui,jn=

values of u at the node (i,j,n) (–)

l=

end of space interval (–)

h=

end of space interval (–)

T=

end of time interval (–)

M1=

number of finite difference in x-coordinate (–)

M2=

number of finite difference in y-coordinate (–)

N=

number of finite difference in t-coordinate (–)

F=

nonlinear objective least-squares function (19) (–)

p=

percentage of noise (–)

σ1,σ2=

standard deviations (–)

QT=

fixed domain (0,h)×(0,l)×(0,T) (–)

Q¯T=

closure of the solution domain QT (–)

Acknowledgements

The comments and suggestions made by the editor and the referees are gratefully acknowledged.

Disclosure statement

No potential conflict of interest was reported by the author(s).

References

  • Cannon JR, Yin HM. A class of nonlinear nonclassical parabolic equations. J Differ Equ. 1989;79:266–288. doi: 10.1016/0022-0396(89)90103-4
  • Cannon JR, Matheson AL. A numerical procedure for diffusion subject to the specification of mass. Int J Eng Sci. 1993;31:347–355. doi: 10.1016/0020-7225(93)90010-R
  • Cannon JR, Lin Y, Wang S. An implicit finite difference scheme for the diffusion equation subject to mass specification. Int J Eng Sci. 1990;28:573–578. doi: 10.1016/0020-7225(90)90086-X
  • Ivanchov MI. Some inverse problems for the heat equation with nonlocal boundary condition. Ukrainian Math J. 1993;45:1066–1071. doi: 10.1007/BF01070965
  • Ivanchov MI. On an inverse problem of heat conduction with nonlocal overdetermination condition. Visnyk L'vivskogo Universytetu Ser Mech Mat. 1994;40:12–15.
  • Ivanchov MI. Inverse problems for equations of parabolic type. Liviv: VNTL; 2003.
  • Huntul MJ, Lesnic D. Reconstruction of the timewise conductivity using a linear combination of heat flux measurements. J King Saud Univ Sci. 2020;32:928–933. doi: 10.1016/j.jksus.2019.05.006
  • Huntul MJ, Lesnic D. An inverse problem of finding the time-dependent thermal conductivity from boundary data. Int Commun Heat Mass Transfer. 2017;85:147–154. doi: 10.1016/j.icheatmasstransfer.2017.05.009
  • Huntul MJ, Hussein MS, Lesnic D, et al. Reconstruction of an orthotropic thermal conductivity. Int J Math Model Numer Optim. 2020;10:102–122.
  • Alifanov OM, Budnik SA, Nenarokomov AV, et al. Estimation of thermal properties of materials with application for inflatable spacecraft structure testing. Inverse Probl Sci Eng. 2012;20:677–690. doi: 10.1080/17415977.2012.665909
  • Huang CH, Huang CY. An inverse problem in estimating simultaneously the effective thermal conductivity and volumetric heat capacity of biological tissue. Appl Math Model. 2007;31:1785–1797. doi: 10.1016/j.apm.2006.06.002
  • Mera NS, Elliott L, Ingham DB, et al. An iterative BEM for the Cauchy steady state heat conduction problem in an anisotropic medium with unknown thermal conductivity tensor. Inverse Prob Eng. 2000;8:579–607. doi: 10.1080/174159700088027748
  • Reddy SR, Dulikravich GS. Simultaneous determination of spatially varying thermal conductivity and specific heat using boundary temperature measurements. Inverse Probl Sci Eng. 2019;27:1635–1649. doi: 10.1080/17415977.2019.1578352
  • Reddy SR, Dulikravich GS, Zeidi SMJ. Non-destructive estimation of spatially varying thermal conductivity in 3D objects using boundary thermal measurements. Int J Thermal Sci. 2017;118:488–496. doi: 10.1016/j.ijthermalsci.2017.05.011
  • Hussein MS, NSEP Kinash, Lesnic D, et al. Retrieving the time-dependent thermal conductivity of an orthotropic rectangular conductor. Appl Anal. 2017;96:1–15. doi: 10.1080/00036811.2016.1232401
  • Sagaydak R. Large on existence and uniqueness of solution for the inverse problem of determination major coefficients in two-dimensional parabloic equation. Visnyk L'vivskogo Universytetu Ser Mech Mat. 2005;64:236–244.
  • Barakat HZ, Clark AJ. On the solution of the diffusion equations by numerical methods. J Heat Transfer. 1966;88:421–427. doi: 10.1115/1.3691590
  • Campbell LJ, Yin B. On the stability of alternating-direction explicit methods for advection-diffusion equations. Numer Meth Partial Differ Equ. 2007;23:1429–1444. doi: 10.1002/num.20233
  • Ozisik MN. Finite difference methods in heat transfer. Boca Raton (FL): CRC Press; 1994.
  • Mathworks. Documentation optimization toolbox-least squares (model fitting) algorithms; 2016. Available from: www.mathworks.com
  • Coleman TF, Li Y. An interior trust region approach for nonlinear minimization subject to bounds. SIAM J Optim. 1996;6:418–445. doi: 10.1137/0806023
  • Cannon JR. The one-dimensional heat equation. Menlo Park (CA): Addison-Wesley; 1984.

Reprints and Corporate Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

To request a reprint or corporate permissions for this article, please click on the relevant link below:

Academic Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

Obtain permissions instantly via Rightslink by clicking on the button below:

If you are unable to obtain permissions via Rightslink, please complete and submit this Permissions form. For more information, please visit our Permissions help page.