Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 48, 1, pp. 111-134, Warsaw 2010 A SEQUENTIAL AND GLOBAL METHOD OF SOLVING AN INVERSE PROBLEM OF HEAT CONDUCTION EQUATION Michał Ciałkowski Poznan University of Technology, Poznań, Poland e-mail: michal.cialkowski@put.poznan.pl Krzysztof Grysa Kielce University of Technology, Kielce, Poland e-mail: grysa@tu.kielce.pl The paper presents a solution to an inverse problem based on the analy- tical form of the direct problem solution in the convolutional form. The analytical form T(r,t) is a surface that is fitted to temperature patterns measured (and charged with errors) at the internal points. In the case of quickly-varying patterns, the solution to the inverse problem is highly sen- sitive to measurement errors (short sampling times). In order to obtain reliable results, the method of sequential (step by step) and global solving of the inverse problemwas used together with smoothing themeasurement results with the help of hyperbolic spline functions. The numerical results confirm effectiveness of the methods presented in the paper. Key words: inverse problem, heat conduction, sequential method, global method 1. Introduction In most heat conduction problems, the governing partial differential equation has to be solved with appropriate initial and boundary conditions. These pro- blems are called direct problems. The problems with incomplete boundary conditions, when the solution is prescribed at some internal points in the do- main, are called inverse problems. Solving the inverse problems is not easy because of their ill-posedness in the Hadamard sense, i.e. the existence, uniqueness and stability of their so- lutions are not always guaranteed. There have appeared various approaches to the inverse heat conduction problem. They are widely reported in Xian- wu and Atluri (2006). Some of the approaches were the following: the use of 112 M. Ciałkowski, K. Grysa Laplace transform, Ciałkowski andGrysa (1980), Grysa et al. (1981), the use of Helmholtz equation, Grysa (1989), the finite difference method Hensel and Hills (1984), the finite element approaches, Hore et al. (1977), Bass (1980), Tikhonov regularization method, Tikhonov and Arsenin (1977), Alifanov ite- rative regularization, Alifanov (1994), approaches based on Trefftz functions, Jirousek (1978), Qing-Hua Qin (2000), Ciałkowski (2006, 2007), Grysa and Leśniewska (2009). In the case of transient heat conduction problems with unknown bounda- ry condition on a part of the boundary, the unknown boundary temperature can be determined after each time step (sequential method) based onmeasu- red temperatures at the internal points of the solution domain. However, the measured temperatures are chargedwith errorswhat usually results in inaccu- rate solution to the inverse problem, in particular for small time steps. In the paper byCiałkowski (2007), splines are used to smooth noisy results ofmeasu- rements. Smoothingwith the use of splines decreases essentially oscillations of the inverse problem solution. In the paper byGrysa and Leśniewska (2009), smoothing of an inaccurate temperature internal response with the use of T-functions appeared to be a very efficient method ensuring stable solution of the considered IHCP. Here we use the Laplace transform technique to find the best approximate solution for the case of inaccurate input data. The analytical form of the solu- tion, T(r,t), can be presented as a surface in the r,t,T system of coordinates. When all results of the measurement are known, one can fit it into all the measured results simultaneously (global method). In transient problems, discretization of the temperature time derivative with the back finite difference may cause instability of the solution for a very small time step. Therefore, a solution to an inverse problem is based here on an analytical formof the direct problem solution. Such an approach is possible only for the simple geometries and for linear heat conduction equation. Here, theproblemofdeterminationof theboundarytemperatureonthe inner surface of the cylindrical layer is considered. 2. The direct problem solution The linear heat conduction equation for a cylindrical layer can be written as ρc ∂T ∂t =λ (∂2T ∂x2 + 1 x ∂T ∂x ) x∈ (R1,R2), t > 0 (2.1) A sequential and global method of solving... 113 In technical applications (cylindrical sleeves, turbine tubes, combustion chambers andothers) the inner surface, x=R1, is a lotmore loaded termically than the outer one, and in the case of thermal shock on the inner surface the change of heat flux on the outer surface, x = R2, is negligible in the initial period of time. Hence, for the heat flux on the outer surface we assume the following condition −λ∂T ∂x ∣ ∣ ∣ x=R2 = q(t)≈ 0 t> 0 (2.2) Fig. 1. Cross-section of the cylindrical layer In the case of heat turbine, condition (2.2) is satisfied thanks to isolation of the turbine tube. However, in general q(t) 6=0. Assume the initial temperature is a constant, i.e. T(x,t=0)=T0 = const x∈ 〈R1,R2〉 (2.3) and the inner surface temperature is a function of time T(x=R1, t1)=Tf(t) t> 0 (2.4) If cooling is considered, then T0 ­Tf. We define dimensionless coordinates ξ= x R2 τ = λ ρc t R22 ϑ= T T0 ϑf = Tf T0 g= R2 λT0 q 114 M. Ciałkowski, K. Grysa Here τ (the dimensionless time) is the Fourier number. Equation (2.1) with conditions (2.2) to (2.4) take the form as follows — the heat conduction equation ∂ϑ ∂τ = ∂2ϑ ∂ξ2 + 1 ξ ∂ϑ ∂ξ ξ∈ (ξ1,ξ2), τ > 0 (2.5) — the initial condition ϑ(ξ,0)=ϑ0 = const (2.6) — the boundary condition on the inner surface, ξ= ξ1 ϑ(ξ= ξ1,τ)=ϑf(τ) τ > 0 (2.7) — the boundary condition on the outer surface, ξ= ξ2 −∂ϑ(ξ= ξ2,τ) ∂ξ = g(τ) (2.8) Tosolve equation (2.5)with conditions (2.6) to (2.8), theLaplace transform is applied (Doetsch, 1964) ϑ(ξ,s)= ∞ ∫ 0 ϑ(ξ,τ)e−sτ dτ (2.9) As a result, we arrive to the following convolutional form of the considered direct problem solution ϑ(ξ,τ)= [ϑ′f(τ)+(ϑf(0)−ϑ0)δ(τ)]∗F(ξ,τ)+ (2.10) −[g′(τ)+g(0)δ(τ)]∗H(ξ,τ)+η(τ)ϑ0 with δ(t) standing for the Dirac function, η(t) being the Heaviside function and F(ξ,τ) = ∞ ∑ n=1 ress=sn 1 s I1( √ s)K0( √ sξ)−K1( √ s)I0( √ sξ) I0( √ sξ1)K1( √ s)+K0( √ sξ1)I1( √ s) est H(ξ,τ) = ∞ ∑ n=1 ress=sn 1 s √ s K0( √ sξ1)I0( √ sξ)− I0( √ sξ1)K0( √ sξ) I0( √ sξ1)K1( √ s)+K0( √ sξ1)I1( √ s) est and sn, n=0,1,2, . . ., being roots of the transcendental equation s ( I0( √ sξ1)K1( √ s)+K0( √ sξ1)I1( √ s) ) =0 A sequential and global method of solving... 115 For √ s= iµ, we obtain µ0 =0 and J0(µnξ1)Y1(µn)=Y0(µnξ1)J1(µn) n=1,2, . . . (2.11) For µξ1 ≫ 1 we arrive at the asymptotic form of equation (2.11), McLachlan (1964) cos[µ(1− ξ1)] = 0 with the roots µn = π 2 +(n−1)π 1− ξ1 n=1,2, . . . (2.12) Roots of the equation cos[µ(1− ξ1)] = 0, µn stand for the first approxi- mation for the roots of equation (2.11) when the Newton method is used to calculate themwith the required accuracy. In particular, we find F(ξ,τ)= 1+2 ∞ ∑ n=1 fn(ξ) e−µ 2 n τ µn (2.13) H(ξ,τ)= ln ξ ξ1 +2 ∞ ∑ n=1 hn(ξ) e−µ 2 n τ µ2n with fn(ξ)= 1 B [J1(µn)Y0(µnξ)−Y1(µn)J0(µnξ)] (2.14) hn(ξ)= 1 B [J0(µnξ1)Y0(µnξ)−Y0(µnξ1)J0(µnξ)] and B= ξ1[J1(µnξ1)Y1(µn)−J1(µn)Y1(µnξ1)]−J0(µnξ1)Y0(µn)+Y0(µnξ1)J0(µn) With regard to the convolutional form of the solution to problems (2.1) to (2.4), the functions ϑf and g that appear in formula (2.10)maybe replacedby their approximation with linear-wise functions in the successive time intervals (comp. Fig.2, where Θ′ = ϑf or Θ ′ = g). Consider the solution ϑ(ξ,τ) for a certain τ = τM < 1, M being a natural number. Because in (2.13) the function exp(−µ2nτ) appears,wefirst calculate the convolution of the functions ϑf and g with exp(−µ2nτ). Using the average value theorem and dividing the time interval 〈0,τM〉 into small subintervals, 〈0,τM〉 = M ⋃ k=1 〈τk−1,τk〉, we obtain 116 M. Ciałkowski, K. Grysa I1 =Θ ′(τ)∗ e−µ2nτ ∣ ∣ ∣ τ=τM = τM ∫ 0 Θ′(η)e−µ 2 n (τM−η) dη=e−µ 2 n τM τM ∫ 0 Θ′(η)eµ 2 n η dη≈ ≈ e−µ2nτM M ∑ k=1 τk ∫ τk−1 Θk −Θk−1 τk− τk−1 eµ 2 n η dη=e−µ 2 n τM M ∑ k=1 Θk−Θk−1 τk− τk−1 τk ∫ τk−1 eµ 2 n η dη≈ (2.15) ≈ e−µ2nτM M ∑ k=1 (Θk −Θk−1)eµ 2 n [τk−1+γk(τk−τk−1)] = M ∑ k=0 ΘkΦkn(τM) where γk ∈ (τk−1,τk),Θk ≡Θ(τk) and Θ=ϑf or Θ= g. Moreover Φkn(τM)=          −e−µ2n(τM−γ1τ) for k=0 B for 0 0 and ξ∈ 〈ξ1,1〉 ϑ(ξ,τM)=ϑf(τM)+ M ∑ k=0 ϑfk ∞ ∑ n=1 fn(ξ) µn Φkn(τM)+ (2.18) − M ∑ k=0 gk ∞ ∑ n=1 hn(ξ) µ2n Φkn(τM)−g(τM)ln ξ ξ1 or, in amore compact form, ϑ(ξ,τM)= {χ(ξ,τM)}⊤{ϑf}−V{ξ,τM}⊤{g} (2.19) where χk(ξ,τM)= ∞ ∑ n=1 fn(ξ) µn Φkn(τM) k=1,2, . . . ,M−1 χM(ξ,τM)= ∞ ∑ n=1 fn(ξ) µn ΦMn(τM)+1 Vk(ξ,τM)= ∞ ∑ n=1 hn(ξ) µ2n Φkn(τM) k=1,2, . . . ,M−1 VM(ξ,τM)= ∞ ∑ n=1 hn(ξ) µ2n Φkn(τM)+ ln ξ ξ1 ϑfM ≡ϑf(τM) Solution (2.10) in form (2.19) is convenient to solve an inverse problem. When calculating numerical sum of an infinite series, the summation is confined to such a number of its elements that leads to inaccuracy of the result not exceeding the assumed value. Let us consider the functions fn(ξ) and gn(ξ), formulas (2.14). For suffi- ciently big values of the roots µn, the functions canbe expressed in asymptotic forms 118 M. Ciałkowski, K. Grysa fn(ξ)≈− √ ξ1 ξ cos[µn(1− ξ)] (1− ξ1)sin[µn(1− ξ1)] gn(ξ)≈− 1√ ξ sin[µn(1− ξ)] (1−ξ1)sin[µn(1− ξ1)] By virtue of (2.12) we obtain sin[µn(1−ξ1)] = (−1)n−1. Hence |fn(ξ)| ¬ √ ξ1 ξ 1 1− ξ1 |gn(ξ)| ¬ 1√ ξ 1 1− ξ1 For convergence of the series in the formula (2.18) the expression Φkn(τM)/µn is crucial. By virtue of (2.16), we find max k¬M 1 µn Φkn(τM)= 1 µn e−µ 2 n (1−γM)(τM−τM−1) 0<γM <τM < 1 Thus, for anypositive εand sufficientlybig µn, byvirtueof the formula (2.12), the following inequality holds 1 µn e−µ 2 n (1−γM)(τM−τM−1) <ε Hence, we can find such a value of the number n for the first series in formula (2.18) that leads to its sufficiently accurate sum. The second series in (2.18) converges as well thanks to the expression µ2n in the denominator. 3. The inverse problem Owing to the practical application, let us a consider an inverse problem for the cylindrical layer. The unknown temperature, ϑf(τk), on the inner surface, for ξ = ξ1, is to be found based on the initial temperature of the layer, the temperature on its outer surface, ξ = 1, and known values of temperature, ϑmeas(ξ ∗ i ,τk), at some inner points, ξ ∗ i , i = 1, . . . ,Mmeas (see Fig.3). The measured values, ϑmeas(ξ ∗ i ,τk) is charged with a noise. In order to find the unknown inner boundary temperature, ϑf(τk), a di- stance between the surface ϑ(ξ,τk) and the given values ϑmeas(ξ ∗ i ,τk) is mi- nimized. The functional that describes the distance (the mean-square error) reads IM = I(τM)= Mmeas ∑ i=1 ‖ϑ(ξ∗i ,τM)−ϑmeas(ξ∗i ,τM)‖2 M =1,2, . . . (3.1) A sequential and global method of solving... 119 Fig. 3. Graphic form of the solution ϑ(ξ,τ) By virtue of the formula (2.19) the functional IM takes the form IM = Mmeas ∑ i=1 ( M ∑ k=0 ϑfkχk(ξ ∗ i ,τM)− M ∑ k=0 gkVk(ξ ∗ i ,τM)−ϑmeas(ξ∗i ,τM) )2 (3.2) In the functional IM, the inner boundary temperature, ϑfM , is an unk- nown.Hence, from the necessary condition of the functional IM minimum, i.e. from the equation ∂IM/∂ϑfM =0, we find ϑfM = M−1 ∑ k=0 ϑfkak+ Mmeas ∑ i=1 bi ( M ∑ k=0 gkVk(ξ ∗ i ,τM)+ϑmeas(ξ ∗ i ,τM) ) (3.3) with ak =− ∑Mmeas i=1 χM(ξ ∗ i ,τM)χk(ξ ∗ i ,τM) ∑Mmeas i=1 χ 2 M(ξ ∗ i ,τM) bi = χM(ξ ∗ i ,τM) ∑Mmeas i=1 χ 2 M(ξ ∗ i ,τM) With ϑfM known from formula (3.3) one can find successive values of the inner boundary temperature for M = 1,2, . . . based on the measured inner point temperature values, ϑmeas(ξ ∗ i ,τk). In this case, the inner boundary temperature is calculated sequentially. Thesequentialmethod is sensitive tonoiseofmeasured temperaturevalues, in particular for initial moments of time, Grysa (1988). The noise present in 120 M. Ciałkowski, K. Grysa the measured data can cause instabilities in the estimated heat fluxes. These instability problems demand proper numerical treatment. The global method seems to be the proper one. In this case, in order to diminish the measurement noise impact on the searched values of the inner boundary temperature, ϑf1,ϑf2, . . . ,ϑfM , we minimize the functional descri- bing themean-value distancebetween the surface ϑ(ξ,τk) and the given values ϑmeas(ξ ∗ i ,τk) in the whole time interval 〈0,τM〉 I = M ∑ k=1 Ik = M ∑ k=1 Mmeas ∑ i=1 ‖ϑ(ξ∗i ,τk)−ϑmeas(ξ∗i ,τk)‖2 (3.4) Using form (3.1) of the function ϑ(ξ,τ) for τ = τM, we obtain I = {ϑf}⊤[χχ]{ϑf}−2{ϑf}⊤([Vχ]{g}+{ϑχ})+C (3.5) Here [χχ] = M ∑ k=1 Mmeas ∑ i=1 {χ(ξ∗i ,τk)}{χ(ξ∗i ,τk)}⊤ [Vχ] = M ∑ k=1 Mmeas ∑ i=1 {V (ξ∗i ,τk)}{χ(ξ∗i ,τk)}⊤ {ϑχ}= M ∑ k=1 Mmeas ∑ i=1 ϑmeas(ξ ∗ i ,τk){χ(ξ∗i ,τk)}⊤ C = M ∑ k=1 Mmeas ∑ i=1 [{V (ξ∗i ,τk)}⊤{g}+ϑmeas(ξ∗i ,τk)]2 The functional I minimized with respect to the vector {ϑf} leads to a system of algebraic equations [χχ]{ϑf}= [Vχ]{g}+{ϑχ} (3.6) The solution reads {ϑf}= [χχ]−1[Vχ]{g}+[χχ]−1{ϑχ} (3.7) In the global method, solving an inverse problem leads to a system of algebraic equations. A sequential and global method of solving... 121 4. Numerical example In order to show the difference between the sequential and the globalmethod, a solution of the inverse problem of cooling the inner surface of a cylindrical layer is presented. The inner radius of the layer, R1 = 0.9m, and outer one, R2 = 1.0m (its thickness is equal to 0.10m). Thermophysical properties of the layer are as follows: λ=30W/(mK), ρ=7800kg/m3, c=500J/(kgK). The following location of two thermal probes are considered: a) s1 =2.0mm and s2 =2.5mm from the inner surface, b) s1 =1.5mm and s2 =2.5mm from the inner surface. The dimensionless initial temperature is assumed as ϑ(ξ,0) = 1. Hence, for T0 =1000 ◦Cwe have T(ξ,0)= 1000◦C. The dimensionless inner surface temperature is assumed in the following form ϑ(ξ1,τ)=ϑf(τ)=ϑmax [ 1− ( 1− ϑmin ϑmax ) ze1−z ] z= τ τmin For τ = τmin we have ϑf(τmin) = ϑmin. For τ ∈ (0,τmin) cooling takes place; for τ > τmin the temperature ϑ(ξ1,τ) increases. Intensity of cooling may be controlled with the minimal temperature value, τmin, and with the moment of time τmin. As the measured data for the inverse problem, numerical data obtained by solving the direct problem can be used. For the direct problem, values of the inner surface temperature, ϑ(ξ1,τ) = ϑf(τ), and the heat flux q(τ) on the outer surface, for ξ2 = 1, are prescribed. We add some random noise on the values of ϑmeas(ξ ∗ i ,τk) whose level is at most 5 ◦C. For this purpose, the noise with normal distribution N(0,0.005) generated numerically in Fortran has been used. The inverse problem has been solved: a) for noise present in the measured data, b) for measured data smoothed with hyperbolic splines. The time step was ∆t=0.2s. In thecase of sequentialmethod, the solution sensitivity for initialmoments of time is a result of solution ϑ(ξ,τ) ”response” for the input measured data with noise. Changes in time of the noise present in the data measured in the first point (the distance from the boundary s1 = 2mm) and in the second one 122 M. Ciałkowski, K. Grysa Fig. 4. Distribution of the random error of measured temperature in the first point for ∆t=0.2s Fig. 5. Distribution of the random error of measured temperature in the second point for ∆t=0.2s (s2 = 2.5mm) are presented in Fig.4 and Fig.5. In Fig.6, the change of temperature in time at the two points is shown. The global solution to the inverse problem for disturbed and undisturbed temperature measurements (here the measurements are simulated numerical- ly) is presented in Fig.7 (the inner boundary temperature) and Fig.8 (heat flux on the inner boundary). The heat flux determination strongly depends on the temperature measurement error. For some moments of time, the heat flux error exceeds 50%. In order to reduce that disadvantageous situation, the measurements are smoothed with hyperbolic splines (Kosma, 1999). The idea of using the hyperbolic splines follows the character of parabolic equation so- A sequential and global method of solving... 123 Fig. 6. Change of temperature in time at two points of measurement for ∆t=0.2s Fig. 7. Temperature on the inner boundary for disturbed and undisturbed temperature measurement at the inner points distant by s1 =2mm and s2 =2.5mm from the boundary for ∆t=0.2s and global approach Fig. 8. Dimensionless heat flux on the inner boundary for disturbed and undisturbed temperature measurement at the inner points distant by s1 =2mm and s2 =2.5mm from the boundary for ∆t=0.2s and global approach 124 M. Ciałkowski, K. Grysa lution, in which a function of hyperbolic type, exp(At)f(x), appears. In the hyperbolic splines, functions exp(As)p(s) are used, with s being time or spa- tial variable and p(s) standing for a polynomial of the third order. For c=0, a hyperbolic spline becomes a classic spline. Figs. 9 to 12 present the temperature and heat flux on the inner boundary obtained as a result of the smoothing. Fig. 9. Distribution of the random error of temperature measurements in the point distant by s1 =2.0mm from the inner boundary Fig. 10. Distribution of the random error of temperature measurements in the point distant by s2 =2.5mm from the inner boundary For the thermoelements distant by s1 = 2.0mm and s2 = 2.5mm from the inner boundary and for the time step ∆t=0.2s, the sequential approach leads to an unstable solution. Displacing the first thermoelement to the point s1 =1.5mm ensured stability of the inverse problem solution. A sequential and global method of solving... 125 Fig. 11. Temperature on the inner boundary for undisturbed and disturbed smoothed temperaturemeasurement at the inner points distant by s1 =2.0mm and s2 =2.5mm from the boundary for ∆t=0.2s and global approach Fig. 12. Dimensionless heat flux on the inner boundary for undisturbed and disturbed smoothed temperature measurement at the inner points distant by s1 =2mm and s2 =2.5mm from the boundary for ∆t=0.2s and global approach The results for the sequential approach for s1 =1.5mm and s2 =2.5mm are presented in Figs. 13 to 15 and for the global approach – in Figs. 16 to 19. Both solutions are compared in Figs. 20 and 21. The calculations show that in the case of sequential approach, placing the thermoelements far from the inner boundarymay lead to an unstable solution for the initial moments of time. The disadvantageous property is not observed in the global approach. 126 M. Ciałkowski, K. Grysa Fig. 13. Change of temperature in time at two points of measurement for ∆t=0.2s Fig. 14. Temperature on the inner boundary for disturbed and undisturbed temperature measurement at the inner points distant by s1 =1.5mm and s2 =2.5mm from the boundary for ∆t=0.2s and sequential approach Fig. 15. Dimensionless heat flux on the inner boundary for undisturbed and disturbed temperature measurement at the inner points distant by s1 =1.5mm and s2 =2.5mm from the boundary for ∆t=0.2s and sequential approach A sequential and global method of solving... 127 Fig. 16. Temperature on the inner boundary for disturbed and undisturbed temperature measurement at the inner points distant by s1 =1.5mm and s2 =2.5mm from the boundary for ∆t=0.2s and global approach Fig. 17. Dimensionless heat flux on the inner boundary for undisturbed and disturbed temperature measurement at the inner points distant by s1 =1.5mm and s2 =2.5mm from the boundary for ∆t=0.2s and global approach Fig. 18. Distribution of the random error of temperature measurements in the point distant by s2 =2.5mm from the inner boundary 128 M. Ciałkowski, K. Grysa Fig. 19. Distribution of the random error of temperature measurements in the point distant by s2 =1.5mm from the inner boundary Fig. 20. Dimensionless heat flux on the inner boundary for undisturbed and disturbed smoothed temperature measurement at the inner points distant by s1 =1.5mm and s2 =2.5mm from the boundary for ∆t=0.2s and sequential approach 5. Location of the measurement points versus accuracy of temperature determination Consider nowhow themeasurement points location affects accuracy of tempe- rature identification. In order to simplify the consideration, the case of a flat slab will be analysed. Conclusions hold also for a thin cylindrical layer which A sequential and global method of solving... 129 Fig. 21. Dimensionless heat flux on the inner boundary for undisturbed and disturbed smoothed temperature measurement at the inner points distant by s1 =1.5mmand s2 =2.5mm from the boundary for ∆t=0.2s and global approach can be approximately described as a flat slab. Instead of equation (2.5) with conditions (2.6) to (2.8) we consider dimensionless equation ∂ϑ ∂τ = ∂2ϑ ∂ξ2 τ > 0, ξ∈ (0,1) (5.1) with conditions ϑ(ξ,0)=ϑ0(ξ)= 0 ∂ϑ ∂ξ ∣ ∣ ∣ ∣ ∣ ξ=0 =0 ϑ(ξ∗,τ)=ϑ∗(τ)> 0 0¬ ξ∗ ¬ 1 Here, for simplicity we put ϑ(ξ,0) = ϑ0 = 0 and g(τ) = 0 (comp. (2.6) and (2.8)), and comparing with (2.7), ξ∗ stands for ξ1 and ϑ ∗ stands for ϑf. In order to find temperature distribution at the beginning of the heating process, we apply the backward finite difference to approximate the time de- rivative in Eq. (2.1) ∂ϑ(ξ,τ) ∂τ ≈ ϑ(ξ,τ)−ϑ(ξ,τ −∆τ) ∆τ =α2[ϑ(ξ,τ)−ϑ(ξ,τ −∆τ)] (5.2) α2 = 1 ∆τ Then, the approximate solution to the problem reads ϑ(ξ,∆τ)=ϑ∗ cosh(αξ) cosh(αξ∗) (5.3) 130 M. Ciałkowski, K. Grysa Consider the influence of the measured temperature noise δϑ∗ and inac- curacy of the thermoelement location δξ∗ on the value of temperature. From (5.3), we find ϑ(ξ,∆τ)+ δϑ(ξ,∆τ)= (ϑ∗+ δϑ∗) cosh(αξ) cosh[α(ξ∗+ δξ∗)] (5.4) Then, subtracting (5.3) from (5.4), we arrive at the temperature error δϑ=ϑ∗ cosh(αξ) cosh(αξ∗) ( cosh(αξ∗) cosh[α(ξ∗+ δξ∗)] −1 ) + δϑ∗ cosh(αξ) cosh[α(ξ∗+ δξ∗)] (5.5) and then we easily find the relative temperature error δϑ ϑ = cosh(αξ∗) cosh[α(ξ∗+ δξ∗)] −1+ δϑ ∗ ϑ∗ cosh(αξ∗) cosh[α(ξ∗+ δξ∗)] (5.6) The first derivative of (5.4) with respect to ξ leads to the following de- scription of the heat flux error δ δϑ(ξ,∆τ) dξ =(ϑ∗+δϑ∗) αsinh(αξ) cosh[α(ξ∗+δξ∗)] −dϑ(ξ,∆τ) dξ =αtanh(αξ)δϑ (5.7) Hence, we easily obtain the relative heat flux error δdϑ dξ dϑ dξ = δϑ ϑ (5.8) The maximum value of the error δϑ appears for ξ = 1, i.e. on the bo- undary with an unknown thermal condition. The relative heat flux error is proportional to the relative temperature error. Notice that the error δϑ/ϑ does not depend on the variable ξ. Consider now a flat slab with thickness h = 100mm made of steel with thermal coefficients c=500J/(kgK), λ=30W/(mK), ρ=7800kg/m3. The time step is ∆t = 0.1s. Then the Fourier number ∆τ = 0.00007623 and α = √ 1/∆τ = 114.0175. Moreover, let the inaccuracy of the thermoelement location be δx∗ =0.2mm and x∗ =99mm. Hence, we obtain ξ∗ =0.99 and δξ∗ =0.002. For the thermal shock (ξ∗ = ξ=1 and ϑ(ξ=1,τ)= 1) we have ϑ(ξ,∆τ)= cosh(αξ)/coshα (comp. (5.3)) and for ξ∗ =0.99, we obtain ϑ∗ =ϑ(ξ∗,∆τ)= cosh(αξ∗) coshα (5.9) A sequential and global method of solving... 131 Hence, for the accurate ϑ∗, we arrive at the formula ϑ(ξ,∆τ)=ϑ∗ cosh(αξ) cosh(αξ∗) = cosh(αξ) coshα = e−α(1−ξ)+e−α(1+ξ) 1+e−2α (5.10) From the input data we find that αξ∗ ∼ 100 and, therefore, we can write with a negligible error cosh(αξ) cosh[α(ξ∗+ δξ∗)] = eα(ξ−ξ ∗−δξ∗) 1+e −2αξ 1+e−2α(ξ ∗+δξ∗) =eα(ξ−ξ ∗−δξ∗) (5.11) Therefore, the errors of the temperature determination read δϑ=eα(ξ−ξ ∗)[ϑ∗(e−αδξ ∗ −1)+ δϑ∗e−αδξ∗] (5.12) δϑ ϑ =e−αδξ ∗ −1+ δϑ ∗ ϑ∗ e−αδξ ∗ For heat flux errors, we obtain δ (dϑ dξ ) =αtanh(αξ)eα(ξ−ξ ∗)[ϑ∗(e−αδξ ∗ −1)+ δϑ∗e−αδξ∗] (5.13) Notice that the temperature distribution as well as heat flux error di- stribution is reinforced exponentially for ξ > ξ∗ and is damped exponen- tially for ξ < ξ∗. The conclusion is that the closer ξ∗ from the boundary ξ = 1, the smaller value achieves the multiplier exp[α(ξ− ξ∗)], ξ > ξ∗. For ξ = 1, α = 114.0175, ξ∗ = 0.99, exp[α(1− ξ∗)] = 3.1273 and for ξ∗ = 0.98, exp[α(1− ξ∗)] = 9.78. For the heat flux, the coefficient of reinforcement is equal to αexp[α(1− ξ∗)]. One more conclusion results from the above analysis, namely that the number of measurement points in the slab (or a cylindrical layer) along the line perpendicular to the boundary should be as small as possible, because – according to the last conclusion – all inaccuracies of themeasurements will be reinforced in points placed closer to the boundary with the unknown thermal condition than the measurement points. Finally, we can conclude that themost important seems to be proper smo- othing of the input data. If the smoothing is appropriate to the considered problem, the results will be encumberedwith an error not exceeding the inac- curacy of the input data. 132 M. Ciałkowski, K. Grysa 6. Conclusions In the paper two methods of solving the inverse transient heat conduction problems are presented and compared. Bothmethods, the sequential and glo- bal one, rely on analytic form of the direct problem solution. The problem is formulated for a flat slab, a cylinder layer and a spherical one, but the detailed consideration is presented for the cylindrical layer only. In the direct problem solution, the unknown temperature of the inner boundary of the cylindrical layer is an element of convolution. In the sequential method, the successive values of the inner boundary temperature are calculated one by one after each measurement of temperature at chosen inner points of the considered body. In the global one, all values of the inner boundary temperature are calculated simultaneously based on all measured temperatures at the inner points. In the case of sequential method, the inner boundary temperature is sensi- tive to temperature measurement errors. Also for the initial moments of time and when the thermoelement is placed too far from the inner boundary, the results are not satisfactory and instability of the solution appears. The global method is free of these disadvantages. For larger time steps, bothmethods lead to practically the same results. From the numerical calculation a rather obvious conclusion follows, name- ly, that the smoothing of the results of temperaturemeasurement is important for the solution quality. The mean-square smoothing of the temperature me- asurementwith the use of hyperbolic splines (Kosma, 1999), leads to a natural approximation of real temperatures because the heat conduction equation so- lution may be expressed with the use of hyperbolic functions (Ciałkowski, 2006). Acknowledgements This work was carried out in the framework of the research projects 3134/B/T02/2007/33 and N51300332/0541, which were financed by the resources for the development of science in the years 2007-2009. References 1. Alifanov O.M., 1994, Inverse Heat Transfer Problems, Springer-Verlag, New York 2. BassB., 1980,Applicationof thefinite elementmethod to thenonlinear inverse heat conduction problem using Beck’s second method, Transaction of ASME, 102, 168-176 A sequential and global method of solving... 133 3. Ciałkowski M., 2006, Uogólnione funkcje cieplne (Generalized thermal func- tions), Zeszyty Naukowe Politechniki Poznańskiej. Maszyny Robocze i Trans- port, 61, 25-37 [in Polish] 4. Ciałkowski M., 2007, Sekwencyjna i globalna metoda rozwiązania zagad- nienia odwrotnego dla równania przewodnictwa ciepła (Sequential and global metod of solving the inverse heat conduction problems),XIII Sympozjum Wy- miany Ciepła i Masy, Koszalin 5. Ciałkowski M., Frąckowiak A., 2000, Funkcje cieplne i ich zastosowanie do rozwiązywania zagadnień przewodzenia ciepła i mechaniki (Heat functions and their application to solving heat conduction andmechanical problems),Wy- dawnictwo Politechniki Poznańskiej 6. Ciałkowski M.J., Grysa K., 1980, On a certain inverse problem of tempe- rature and thermal stresses fields,Acta Mechanica, 36, 169-185 7. Doetsch G., 1964, Guide to the Applications of Laplace Transforms, PWN Warszawa [Polish translation] 8. Grysa K., 1988, Uwagi o stabilności rozwiązań pewnych jednowymiarowych zagadnień odwrotnych przewodnictwa cieplnego (Remarks on stability of one- dimensional inverse heat conduction problem solutions), Prace Naukowe Poli- techniki Lubelskiej, 167,Mechanika, 39, 5-28 [in Polish] 9. Grysa K., 1989, On the exact and approximate methods of solving inverse problemsof temperaturefields,Rozprawy, PolitechnikaPoznańska,204,Poznań [in Polish] 10. Grysa K., Ciałkowski M.J., Kamiński H., 1981, An inverse temperature field in the theory of thermal stresses,Nucl. Engng. Design, 64, 2, 169-184 11. GrysaK.,LeśniewskaR., 2009,Differentfinite elementapproaches for inver- se heat conduction problems, Inv. Problems Sci. Eng., accepted for publication 12. Hensel E.C., Hills R.G., 1984,A spacemarching finite difference algorithm for the onedimensional inverse conductionheat transfer problem,ASMEpaper, No. 84-HT-48 13. HoreP.S.,KruttzG.W., SchoenhalsR.J., 1977,Application of the finite elementmethod to the inverse heat conduction problem,ASME paper, No. 77- WA/TM-4 14. Jirousek J., 1978, Basis for development of large finite elements locally satis- fying all field equations,Comp. Meth. Appl. Eng., 14, 65-92 15. Jirousek J., Wróblewski A., 1996, T-elements: State of the art and future trends,Arch. Comp. Meth. Eng., 3, 323-434 16. KosmaZ., 1999,Metody numeryczne dla zastosowań inżynierskich (Numerical methods for engineer applications), Politechnika Radomska, Radom [in Polish] 134 M. Ciałkowski, K. Grysa 17. McLachlan N.W., 1964, Bessel Functions for Engineers, PWN Warszawa [Polish translation] 18. Qing-HuaQin, 2000,The Trefftz Finite and Boundary ElementMethod.WIT Press, Southampton 19. Tikhonov A.N., Arsenin V.Y., 1977, Solution of Ill-posed Problems, Wiley & Sons,Washington, DC 20. Xianwu Ling, Atluri S.N., 2006, Stability analysis for inverse heat conduc- tion problems,Comput. Modeling Eng. Sci. (CMES), 13, 3, 219-228 Sekwencyjna i globalna metoda rozwiązania zagadnienia odwrotnego dla równania przewodnictwa ciepła Streszczenie W pracy przedstawiono rozwiązanie zagadnienia odwrotnego w oparciu o anali- tyczną postać rozwiązania zagadnienia prostegow postaci splotowej. Ta postać anali- tyczna T(r,t) jest powierzchnią, która jest dopasowywanadoprzebiegówtemperatury pomierzonej (obciążonej błędem)w punktachwewnętrznych.Dla przebiegów szybko- zmiennych rozwiązanie zagadnienia odwrotnego jest bardzowrażliwena błędypomia- rów (małe czasy próbkowania). Dla otrzymania wiarygodnychwyników zastosowano metodę sekwencyjnego (z kroku na krok) i globalnego rozwiązywania zagadnienia od- wrotnego w połączeniu z wygładzaniem wyników pomiarowych za pomocą hiperbo- licznych funkcji sklejanych.Wyniki obliczeń numerycznych potwierdzają efektywność przedstawionychmetod. Manuscript received January 22, 2009; accepted for print June 17, 2009