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.
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 and of an orthotropic rectangular conductor along with the temperature in a two-dimensional problem given by the parabolic heat equation from non-local over-specification conditions and where , are fixed values from the intervals and , and h, l are given positive quantities. When and , the above non-local over-specification conditions become (Equation5(5) (5) ) and (Equation6(6) (6) ), respectively. In this paper, we generalize the above non-local over-specifications to the local heat flux measurements (Equation5(5) (5) ) and (Equation6(6) (6) ) 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 and in the two-dimensional heat equation [Citation16], (1) (1) where , are known space variations, , are known velocity representing convective flow, is a known absorption coefficient, f is a known heat source, u is a unknown temperature in the rectangular domain , subject to the initial condition (2) (2) the Dirichlet boundary conditions (3) (3) (4) (4) and the local heat flux over-specification conditions (5) (5) (6) (6) where , are fixed values from the intervals and , where h, l are given positive numbers, and , are given functions satisfying compatibility conditions. The geometry of the inverse problem under investigation is shown in Figure .
The local existence and uniqueness of solution of the inverse problem (Equation1(1) (1) )–(Equation6(6) (6) ) were established in [Citation16] and read as stated in the following two theorems.
Theorem 2.1
Suppose that the following assumptions are satisfied:
(A1) | |||||
(A2) | |||||
(A3) | |||||
(A4) |
Then, there exists which is determined by the input data, such that the problem (Equation1(1) (1) )–(Equation6(6) (6) ) has a solution with for i = 1, 2.
Theorem 2.2
Suppose that the condition i = 1,2 for is satisfied. Then, the inverse problem (Equation1(1) (1) )–(Equation6(6) (6) ) cannot have more than one solution in the class with for i = 1, 2.
3. Numerical solution of direct problem
In this section, we consider the direct initial boundary value problem (Equation1(1) (1) )–(Equation4(4) (4) ), where , , , , i = 1, 2, and i = 3, 4 are known and the solution is to be determined together with the quantities of interest . We subdivide the solution domain into , and N subintervals of equal step lengths and where and respectively. At the node we denote where , , , , , , , and for , .
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(1) (1) )–(Equation4(4) (4) ) the ADE is as follows.
Let and be the solutions of the following equations which are multilevel finite difference discretization of Equation (Equation1(1) (1) ): (7) (7) (8) (8)
Furthermore, let the and also satisfy the initial and boundary conditions (Equation2(2) (2) )–(Equation4(4) (4) ), namely, (9) (9) (10) (10) (11) (11) Rearranging the terms in (Equation7(7) (7) ) and (Equation8(8) (8) ), we obtain the explicit calculations of and as follows: (12) (12) (13) (13) where (14) (14) From (Equation12(12) (12) ) and (Equation9(9) (9) )–(Equation11(11) (11) ) for , 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 , and will be known from initial and boundary conditions (Equation9(9) (9) )–(Equation11(11) (11) ). Similarly, can be calculated explicitly from (Equation13(13) (13) ) and (Equation9(9) (9) )–(Equation11(11) (11) ) for , beginning at the boundaries x = 1 and y = 1 and marching in a sequence of decreasing i and j, i.e. , These values are then substituted into the simple arithmetic mean approximation (15) (15) to obtain the solution .
The local heat flux over-specification conditions (Equation5(5) (5) ) and (Equation6(6) (6) ) can be calculated using finite difference method approximation as follows: (16) (16) (17) (17)
4. Numerical solution of inverse problem
In this section, our goal is to obtain simultaneously stable identifications for the timewise thermal conductivity coefficients and , satisfying Equations (Equation1(1) (1) )–(Equation6(6) (6) ). We reduce the inverse problem to a nonlinear minimization of the least-squares objective function (18) (18) or, in discretizations form (19) (19) where solves numerically using the ADE method the direct problem (Equation1(1) (1) )–(Equation4(4) (4) ) for given and respectively. It is worth mentioning that in (Equation19(19) (19) ) at the first time step, i.e. n = 0, the derivative and are obtained from the initial condition (Equation2(2) (2) ), via (Equation16(16) (16) ) and (Equation17(17) (17) ), as (20) (20) (21) (21) where for , . Also, one can remark that at initial time the values and can be obtained from the local heat flux over-specification conditions (Equation5(5) (5) ) and (Equation6(6) (6) ) as (22) (22) The minimization of the objective function (Equation19(19) (19) ) 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(1) (1) )–(Equation6(6) (6) ) is solved subject to both exact and noisy measurements (Equation5(5) (5) ) and (Equation6(6) (6) ). The noisy data is numerically simulated as (23) (23) where and are random variables generated from a Gaussian normal distribution with mean zero and standard deviations and , respectively, given by (24) (24) where p represents the percentage of noise. We use the MATLAB function normrnd to generate the random variables and as follows: (25) (25) In the case of noisy data (Equation23(23) (23) ), we replace and by and , respectively, in (Equation19(19) (19) ).
5. Numerical results and discussion
In this section, we present numerical results for the thermal conductivity and together with the temperature in the case of exact and noisy data (Equation23(23) (23) )–(Equation25(25) (25) ). To quantify the accuracy of the approximate solution, we employ the root mean square errors (rmse) defined by (26) (26) (27) (27) For simplicity, we take h = l = T = 1. The lower bounds and upper bounds for the coefficients and are for (LB) and for (UB). These bounds allow a wide search range for the unknown positive physical components and .
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 and depends on time (and are unknown) occur in damage and radioactive decay applications [Citation22].
5.1. Example 1
Consider the inverse problem (Equation1(1) (1) )–(Equation6(6) (6) ) with unknown coefficients and , with the input data: (28) (28) (29) (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) (30) (31) (31) First, we assess the accuracy of the direct problem (Equation1(1) (1) )–(Equation4(4) (4) ) with the input data (Equation28(28) (28) ) when and are known and given by (Equation31(31) (31) ). The numerical results for the interior temperature have been obtained in excellent agreement with the analytical solution (Equation30(30) (30) ) and therefore they are not presented. Apart from the interior temperature, other output of interest is the data (Equation5(5) (5) ) and (Equation6(6) (6) ), which analytically is given by (Equation29(29) (29) ). Figure shows that the analytical and numerical solutions for this quantity obtained with and with various numbers of time steps are in very good agreement. Also, the rmse defined by (32) (32) (33) (33) indicated in Table , show more clearly the convergence of the numerical ADE solution to the analytical solution (Equation29(29) (29) ).
Table 1. The rmse values ((Equation32(32) (32) ) and (Equation33(33) (33) )) for and , with and with various , for direct problem.
Next we investigate the inverse problem. We fix and N = 40 and we start the investigation for reconstructing the unknown timewise thermal conductivity coefficients and and the temperature for noise in the measured data (Equation23(23) (23) ). We take the initial guesses for the vectors and as follows: (34) (34) Note that the values of and are available from (Equation22(22) (22) ). The objective function (Equation19(19) (19) ), 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 in about 6 to 8 minutes CPU time. The corresponding numerical results for and 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 to to 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(30) (30) ) and numerical solutions for 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 , 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 and will become apparent and regularization would need to be employed.
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 , for Examples 1 and 2.
5.2. Example 2
The previous example has recovered the smooth timewise thermal conductivity coefficients and as given in (Equation31(31) (31) ). In this example, we assess the numerical ADE scheme for reconstructing a non-smooth test: (35) (35) We also have the same analytical solution for the temperature given by (Equation30(30) (30) ). The rest of the input data are the same as in Example 1 and (36) (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 and for this example has been taken as (37) (37) Note that the values of and are also available from (Equation22(22) (22) ). We investigate the inverse problem as we did in Example 1. We take and N = 40 and reconstructing the unknown thermal conductivity coefficients and and the temperature for exact and noisy measured input data (Equation23(23) (23) ), for various amounts of noise p = 0 (no noise) and noisy data (Equation24(24) (24) ). The objective function (Equation19(19) (19) ), 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 . The numerical results for the timewise thermal conductivity coefficients and 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 to to and then to zero the numerically obtained results becomes more stable and accurate.
Finally, Figure shows the analytical (Equation35(35) (35) ) and numerical solutions for the temperature for various amounts of noise . 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.
6. Conclusions
In this paper, the inverse problem involving simultaneous identification of the timewise thermal conductivity coefficients and and the temperature 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
= | thermal conductivity (W m K ) | |
= | thermal conductivity (W m K ) | |
= | velocity (m s ) | |
= | velocity (m s ) | |
= | absorbance (AU) | |
f | = | heat source/force (J) |
u | = | temperature (K or C) |
= | heat flux at x = 0 (W ) | |
= | heat flux at y = 0 (W ) | |
t | = | time variable (s) |
x, y | = | space variables (m) |
= | space nodes (–) | |
= | time node (–) | |
= | values of u at the node (–) | |
l | = | end of space interval (–) |
h | = | end of space interval (–) |
T | = | end of time interval (–) |
= | number of finite difference in x-coordinate (–) | |
= | 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 (–) |
= | standard deviations (–) | |
= | fixed domain (–) | |
= | closure of the solution domain (–) |
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.