Back to Inverse Heat Conduction: Ill-Posed Problems Home
Table of Contents
SEC . 6.2 C HAPTER 6 D IFFERENCE M ETHODS FOR T HE S OLUTION OF T HE O NE-DIMENSIONAL I NVERSE H EAT C ONDUCTION P ROBLEM 6.1 - ~ I NTRODUCTION Chapter 3 indicated that difference methods offer considerable potential for solving a wide variety o f nonlinear direct heat conduction problems. This statement is also true for the inverse heat conduction problem. In particular, the nonlinearity associated with temperature dependent thermal properties can be easily accommodated with difference methods. The recommended procedure is to locally linearize the problem by evaluating all thermal properties at temperatures corresponding to the previous time step. Thermal properties are seldom known to sufficient accuracy to justify iteration. This approach is called quasi-linearization. Chapter 6 focuses entirely on difference methods for the solution of the IHCP. The techniques presented are equally applicable to lumped capacitance, finite control volume, and finite element techniques. Section 6.2 expands earlier discussions on sensitivity coefficients and how they can be efficiently calculated by difference methods. Section 6.3 considers exactly matching the calculated temperatures with measured values for a single sensor for each time step; Section 6.4 considers mUltiple sensors. Section 6.5 presents techniques for simultaneously estimating several heat flux components. Sections 6.6, 6.7 and 6.8 apply difference techniques to the function specification method. Section 6.9 considers the sequential regularization method. Section 6.10 introduces space marching techniques. Section 6.11 gives several examples and Section 6.12 lists some computer programs. 218 2 19 S ENSITIVITY COEFFICIENTS 6.2 S ENSITIVITY COEFFICIENTS A ND THEIR CALCULATION BY DIFFERENCE M ETHODS The concept of sensitivity coefficients was introduced in C hapter 1 and has been used extensively throughout this book. Since the key to understanding the material in this chapter is a thorough understanding of sensitivity coefficients, some concepts about them will be expanded. Two sensitivity coefficients occur repeatedly in this chapter: (1) that due to a step qM in heat flux for infinite time duration and (2) t hat due to a pulse qM in heat flux for time duration t M - tM - I ' B~th the heat flux time variation and the corresponding sensitivity coefficient variation with time are shown schematically in Figure 6.1. Actual step function sensitivity coefficients are shown in Figures 1.7 and 1.11 for planar slabs with an insulated surface and slabs that are semi-infinite, respectively. Since the duration of the step function is infinite, the step function sensitivity coefficient grows without bound. The curve labeled A in Figure 6.1 is for a step that begins at time t M-Io t hat labeled C is for the same step but beginning at time t M' The symbol Z is used to denote the step function sensitivity coefficient; this Z is the same as the 4> in Chapters 3 -5 i fthe thermal properties are constant. F or a sensor a t depth X l, the temperature response due to the step in heat flux can be written as T(Xko t , t M - I , q M-I, q M); t M - I denotes the time a t which the heat flux step o r pulse begins. The components of the heat flux vector qM - I , which are <1 I, <12' ... , <1M - I' contain al\ the information about the time variation of the heat flux prior to tM - 10 and hence contain all the (initial) temperature information necessary to continue the temperature calculations to t M , tM + I , • • • . The step function sensitivity coefficient is defined by OT(Xko t , tM - 10 11M - I, q M) (6.2.1) OqM aT aq q (l) 1 .,-1 I ., 1 .,+1 Time FIGURE 6.1 Schematic o f sensitivity coefficient for step in heat ftux a nd pulse in heat ftux. A, sensitivity coefficient for step in heat ftux at time 1 .,_ , . B, sensitivity coefficient for pulse in heat ftux, beginning at time 1., _,. C , sensitivity coefficient for step in heat ftux at time I.,. SEC. 6.2 220 C HAP.6 221 SENSITIVITY COEFFICIENTS O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM profile at tM - 1 • Using Eq. (6.2.2) and the notation of Eq. (6.2.4), one can show that for constant properties, is sensor location is arbitrary time is the time at the beginning of the heat flux step is the previous heat flux history is the magnitude o fthe heat flux step where: r Z l,r= 1, q M_.)=Z(Xl, t, t M-l' ~-I) - Z[Xlo t -(tM -tM- 1 ), t M-l, q M-t] which simply states that for the linear problem the response due to a step is the sum of the response due to a series o f pulses of constant magnitude and distributed over the same time interval. F or nonlinear problems, the sensitivity coefficients can be calculated from differences of the temperatures at the same time and location for two values of the heat flux. F or example, the temperature distribution is calculated for two different values of heat flux, q . and (1 + e)q·, where e is a small parameter o f ord~r 0.001. Then the sensitivity coefficient is given by X ( Xl, t, t M - A T.=.[x...::.:.,...:.t,-...:tM =---.:,:I,-...:q=::M:...-...:I.:.-'q..:...·....:.(_I_+-'e):,=.]_-_T_(:...X::..:1':...-t:..'.=M:...--...!.I.:.-'q..=:M:!..-_I!..:.'...!.q....:,·) t 1 eq· ) __ 1 , qM - 1 - (6.2.5) The only difference in the X and Z sensitivity coefficient calculations is the boundary condition: Sensitivity Coefficient Boundary Condition X(pulse) q(t) = { q., Z(step) q (t)=q·, (6.2.2) The two terms on the right of Eq. (6.2.2) are shown in Figure 6.1 as curves A and C ' for constant thermal properties, curve C is identical to curve A except it is di~placed in time by an amount tM -tM - 1 • Curve B is ( A-C) and is the pulse sensitivity coefficient X which is analogous to AljJ o f Chapters.3 to 5. When using difference methods, discrete values of X and Z are Important. At time t =tM' (6.2.3) because Z (Xl' t M- 1 , t M- 1 , qM- .)=0. F or the constant property case, a convenient subscript notation can be used. Z l.j-M+l = Z(Xl, t j -t M- X1,j j= 1 F or a linear problem, Z is independent of qM; that is, T is linear in qM' There are applications in which the thermal properties are held constant a few (typically 1- 10) future time steps while allowing some property ~anatlon over the entire problem time. F or this case, Z is independent of qM b ut It does depend on t - and qM-l through the temperature profile at t M-t· I f the thermal M1 properties are independent of temperature over the entire problem time, then Z is independent of qM a nd q M-l and the important time var~a~l.e is t - t M.- 1 , the time from initiation of the step. These three important sensltlVlty coefficient cases are summarized in Table 6.1. By definition, a sensitivity coefficient is identically zero prior to the initiation of the heat flux step o r pulse. The pulse sensitivity coefficient will be denoted by X with th~ s~me arguments as those for the step function sensitivity coefficient. F or quasl-hnear and constant property cases, superposition can be used to calculate X from Z. X (Xl' t, t M- L (6.2.4) 1) This same notation will also be used for the quasi-linear property case where it is implicitly understood that all properties are evaluated at the temperature 0, Case II I II None, general case Quasi-linear thermal properties Constant thermal properties - k a T! - =q(t) . <=0 az = ~ (k az) at = qdt) . <=L T (x, tM - d (6.2.6) - k a z! =1 ax . <=0 ax ax qM-I' qM- Z (Xh t , t M-I' q M) Z (Xh t , t M-1t d Z (Xh t -t M - d=4>(Xh t -t M- t~tM-l Sensitivity Coefficient, Z - k a T! Notation ~t~tM Temperature, T N otation f or S tep F unction S ensitivity C oefficients Restrictions 1 The evaluation of Eq. (6.2.5) requires two different temperature calculations (or direct solutions). I t is computationally more efficient to calculate X directly from its own partial differential equation, which can be derived from the transient heat conduction equation. F or example, consider the problem discussed in Section 3.3 for the quasi-linear case· for t > tM - 1 ; pc T ABLE 6.1 tM - all other t d = F (x) ax i- kJZ! aX ax (6.2.7) =0 . <=L Z (x, t M-I)=O · The quasi-linear case implies that k, p, a nd c are evaluated for the temperature profile at 1 ",-1' 222 CHAP. 6 ONE- DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM A comparison o fthe above two sets of equations for T a nd Z indicates that they both satisfy the same partial differential equation. The boundary and initial conditions are different and those for Z are considerably simpler. If X is to be calculated instead o f Z, then the foregoing x =o b oundary condition must be modified as follows - k OXI = {1 for pulse duration o x .. = 0 0 for all other times (6.2.8) The differential equation and the remaining boundary and initial conditions for X are identical to those for Z. D ue to the similarity o f the T a nd Z differential equations. considerable computational savings are possible. This point is explored in detail in a subsequent section. 6.3 S INGLE TEMPERATURE SENSOR, FUNCTION SPECIFICATION ( q=C), S INGLE F UTURE T IME STEP ( EXACT M ATCHING O F D ATA) A single temperature sensor is considered to be located at a depth X l below the active surface. If the heat flux qM is constant over the time interval tM - I ~ t~tM' the value of qM t hat forces a matching of the computed temperature at X l with 'the measured temperature can be calculated. F or realistic temperature d ata containing errors and small dimensionless times, this approach produces significant oscillations in the computed heat flux and would probably not be used in practice. The method is studied here because it aids in the understanding of inverse problems and also because an exact matching of temperature data is a building block that can be used in the development of other techniques. 6.3.1 M odification o f D ifference E quations o f t he D irect H eat C onduction P roblem f or t he S olution o f t he I HCP Difference methods for solving direct heat conduction problems were introduced in Section 3.3. By analogy with Eq. (3.3.45), the difference equations for a onedimensional direct heat conduction problem with a nonuniform grid can be written in symbolic form as SEC . 6.3 2 23 SINGLE T EMPERATURE SENSOR The coefficients ai' bi , Ci ' ai, bi, a nd ci depend on grid spacing and thermal properties a nd a re independent of the constant val~e o f ~eat flux qM for the linear and quasi-linear cases. F or the purpose o f diSCUSSion, only five nodes are used in what follows. Equation (6.3 .1) can be written more compactly as bl - C2 o o o -al b2 - c] 0 0 0 - a2 0 0 0 b] - a] - C4 0 b4 0 0 - a4 -C s bs Tf T 't d; +glqM d2 T~ d] T~ d4 ds T 't (6.3.2) The d's represent all of the right-hand side o f Eq. (6.3.1) except the. qM term. In the direct problem the T~,j=I, . .. , 5 values are found from a Simultaneous solution of Eq. (6.3.2). N ote the explicit dependence of the temperature o n the heat flux qM' . Instead o f solving the direct problem given by Eq. (6.3 .2), suppose qM IS a n unknown and a temperature sensor is located at node 3 ( x./L=O.5). T he direct problem unknown T~ is now a known measured temperature YM • F or t he IHCP, Eq. (6.3.2) becomes bl - C2 0 0 0 - al b2 - C] 0 0 gl 0 0 0 0 0 0 - a] b4 - Cs 0 0 0 - a4 bs Tf T 't qM T~ T 't d; d 2+a2 YM d ]-b]YM d4 +C4YM ds (6.3.3) Equation (6.3.3) represents a system of five linear algebraic equations for .the five unknowns ( Tf, T 't, qM, T~, T 't) a nd can be solved by any approprIate technique for linear algebraic equations. As indicated in previous, this formulation is linear in the unknown heat flux a nd s hould n ot reqUIre iteration. I t h as been demonstrated by Beck and Wolfl t hat the foregoing formulation is unstable for 0 = 1/2 a nd lumped thermal capacitance. Hence, when the d ata from a single temperature sensor are matched exactly, 0 = 1/2 should be avoided. Alifanov l6 used the same technique with 0 = 1. (See Sect. 3.3.3 for definition of 0.) 6.3.2 S ensitivity C oefficient A pproach f or E xactly M atching D ata F rom a S ingle S ensor (6.3.1) In Section 6.3.1 it was demonstrated that matching d ata exactly from a single temperature sensor produces a system of linear algebraic equations. ~~wever, the system of equations represented by Eq. (6 .3.3) is n o longer of the trIdiagonal form. (See Section 3.3 .4.) Hence, a direct solution of Eq. (6.3.3) is n ot the most computationally 'efficient approach. Because o f the linearity of the I HCP, it is possible to utilize the tridiagonal structure of the algebraic equations (for onedimensional problems) and develop a very efficient algorithm. 2 24 C HAP.6 O NE-DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM The temperature field T (x, t) depends in a continuous manner on the unknown heat flux qM (constant over the interval t M- I ~t~ tM)' This dependence is written as T (x, t, tM_ I , qM _ I , qM) where qM _ I is the vector of all previous heat flux values and tM - I indicates the time that" the heat flux step begins. Because the temperature field is continuous in qM, it can be expanded in a Taylor series about an arbitrary but known value o f heat flux q *, T (x, t, t M-I' qM-I , q M)= T (x, t, t M- I , ~-lo q *) + (q M-q S EC.6 .3 known, the complete temperature field can be calculated from Eq. (6.3.5). Equation (6.3.8) is equivalent to the Stolz equation given by Eq. (5.3.1) if q * =0. An efficient way o f calculating sensitivity coefficients is to solve the differential equation for X that is analogous to Eq. (6.2.7) . Since the T and X calculations are similar, considerable computational savings are possible. The development o f a n efficient algorithm requires a thorough understanding of the following tridiagonal matrix algorithm (TDMA): Linear eqUIJtion lorm : b l ul-a l u2=d l *) o T(x, t, t M-1! qM-I ' qM)\ OqM t M='. - Cjuj_l+bjuj-ajuj+l=dj , j =2,3, . .. , N-l + ( qM-q*)2 o2T(x, t, tM - I , q M-lo qM)\ + 2! o ql, t M=,. ... (6 .3.4) Forward elimination (starting a t node N ) : 1 WN= b ' eN=WNcN, N T (x, t M , t M-1! qM-l> q M)= T (x, t M , tM- I , ~-I' q *) W j=b (6.3.5) ' t M, M -I' M-I - !l uqM 'M='· (6.3.6) Note that X does not depend on qM for a linear problem. Equation (6.3.5) is a prescription for calculating the temperature field corresponding to qM provided the temperature field is known for a given q* a nd a means exists for calculating the sensitivity coefficient X. I f the problem is amenable to an analytical solution, X can be calculated as indicated in Section 1.6. Difference methods for obtaining X are discussed in Section 6.2 in connection with Eq. (6.2.5). F or the I HCP in which the sensor data are matched exactly, the left-hand side ofEq. (6.3.5) is replaced by the experimental temperature; hence, Eq. (6.3.5) becomes y M=tt' + (qM-q*)X t . 1 (6.3.7) Solving Eq. (6.3.7) for the estimate of the unknown heat flux, qM, gives: A q M=q IN=wNdN . J =N-l,N-2, . .. , 1 e j=wjcj , j =N-l,N-2, . .. , 2 (6.3.10) I I = wI(d l + ad2) ) _OT(X,tM,tM - I,qM-I,qM)\ q 1 ' j -ajej+1 J j=wj(dj+aJj+d, where the sensitivity coefficient is defined by X (x t (6.3.9) -CNUN-I +bNuN=dN F or linear problems, only the first derivative in Eq. (6.3.4) is nonzero, thus the following is a n exact result for location x a t time tM + (qM-q*)X(x, t M , t M- I , q M-d 226 SINGLE TEMPERATURE SENSOR * + y M-tt' X (6.3.8) t .1 where t t' is the temperature at tM for the sensor node k with q = q * over tM - I ~ t ~ t M • T he calculational procedure is to assume an arbitrary value of heat flux q *, calculate the temperature field T (x, tM, tM- lo q M-I' q*), and knowing the sensitivity coefficients X (x, tM, t M- lo qM-I), calculate the heat flux qM from Eq. (6.3.8) t hat exactly matches the temperature data YM • Once qM is Back substitution : u l=il u j=ejuj_I+Jj, j =2,3, . .. , N (6.3.11) In comparing the two problems for T and X, the following conclusions can be drawn: T aj' bj , Wj' ej Cj X +-same for b oth+-same for both because t hey_ depend only o n aj' bj, C j aj, bj , Wj ' ej Cj dl = gl dj =O,j=2, 3, . . . , N d l =d'i + glqM dj"/=O,j=2, 3, . .. , N Note that the difference equations for X can be derived by differentiating the difference equations for T with respect to qM ' G oing through the forward elimination portion of the algorithm for X, it can be shown that all the J.'s are identically zero except for I I a nd • I I = glwi (6.3.12) Consequently, the sensitivity coefficients are calculated as follows : X I.I = /I=glwl (6.3.13) 226 C HAP.6 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM SEC. 6.3 X j.l=ejXj _ I.I, (6.3.14) j =2,3, . .. , N The algorithm for exactly matching the temperature data from a single sensor can be summarized as follows: 1. Assume an arbitrary value o f q* and calculate the entire temperature field T =T(x , tM ' tM-1> q M-I, q*). q *=qM-1 is a common choice; however, q* = 0 is also perfectly acceptable for the linear problem and will reduce the total number of calculations. 2. Using the same matrix coefficients a n, bn , Cn , W n , and e n as were used for the T calculation in step 1, calculate the sensitivity coefficients from Eqs. (6.3.13) a nd (6.3.14) (6.3.13) (6.3.14) 3. Calculate the heat flux t hat exactly matches the experimental temperature d ata by using the Taylor series expansion OM l q M=q * + YM - T X ( 638) .. A 227 SINGLE T EMPERATURE SENSOR N a , N m , and Nd represent the cycle times required for a single addition, multiplication, and division, respectively. Using the information in Tables 6.2 and 6.3, the approximate fractional increase in computation time for the inverse problem relative to the direct problem is given by: ~T 2 (N + l )+ 2N(Tm/Ta) + (Td/Ta) T 3(N - 1)+5(N -1)(Tm/Ta)+N(T,,/To) "Computation index" is probably a more appropriate term for ~T/T because Eq . (6 .3.16) does not account for the fact that some cycle time is lost between the end of one operation and the beginning o f the next operation. F or a particular computer, ~T/T depends only on the total number o f nodes N. Using the computation times given in Table 6.3, ~T/T for various values of N are presented in Table 6.4. The I HCP causes less than a 40% increase in computations in comparison to the direct problem. Note that Eq. (6.3.16) does not depend on the node number of the temperature sensor; a temperature sensor a t node 1 requires the same number of computations as a backface sensor at node N. Some computational improvement is l .1 4. sion Calculate the complete temperature field from the Taylor series expan(6.3.15) Table 6.2 presents an arithmetic operation count to determine the relative increase in computations o f the inverse problem over the direct problem. The operational counts can be used to evaluate the computational efficiency for a particular algorithm. Table 6.3 presents the number of cycles on the C DC 7fIXJ computer required for the four arithmetic operations; one cycle is 27 ns. Let T ABLE 6 .2 O perational C ounts f or O ne T ime S tep o f I nverse A lgorithm t hat E xactlv M atches S ingle T emperature S ensor D ata ( N=total n umber o f n odes) E quation + o r- x (6.3.10) (6.3.11 ) S tep 3(N-l) S (N-I) C omments N D irect calculation T ABLE 6 .3 R elative S peeds o f V arious A rithmetic O perations o n t he C DC 7 600 C omputer 2 O peration N o . o f M achine Cycles Relative T ime ( f/r.) 4 5 20 1.0 1.25 5.0 + o rx T ABLE 6 .4 C omputation I ndex f or E xactlv M atching D ata f rom S ingle T emperature S ensor ( for C DC 7 600 c omputer. T able 6 .3). N o. o f e lements=N-1 N 2 - ~ ~ (6.3.13) (6.3.14) 3 (6.3.8) 2 4 S um(2+3+4) (6 .3.15) 2N 2 (N+ I ) 4 .5N+7 t n/f= - - - 14.25N - 9 .25 11 21 31 41 81 0 .383 0.350 0.339 0.333 0.324 N N 2N Inverse calculations only (6.3.16) 228 CHAP. 6 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM possible provided the algorithm is restructured. This restructuring is possible because it is not necessary t o k now the sensitivity coefficients X i , l a nd t~ for x > X x (the sensor depth). I f t he sensor is n ear node 1, t he active surface, a restructuring c an offer considerable savings. T he c omputational algorithm is r eordered in the following manner: 1. Assume a value for q* a nd perform the forward elimination portion o f t he tridiagonal matrix algorithm. (Note t hat t he only calculation in this step t hat d epends o n t he heat flux is f l ') 2. C alculate t~ from the back substitution portion o f t he tridiagonal matrix algorithm. t il = fl OM T J =eJ t M-I+fj, J j =2,3, . .. , K (6.3.17) X i,l 229 SINGLE TEMPERATURE SENSOR T ABLE 6.5 O perational C ount f or R estructured A lgorithm t hat E xactly M atches S ingle T emperature S ensor D ata ( K=sensor n ode n umber. N =number o f n odes. K :r;,N) Step E quation 1 5 1 +5 3 (6.3.13) (6.3.14) (6.3.8) C omments x + o r- (6 .3.11) (6.3.10) O ne equivalent direct problem solution 5 (N-1) 3 (N-I) N ;---------(6~;. ;;) -----~ ~-; -----~ ~-I- -----------t~ -c~~~~I~;i~~ - - - - - - - - - - -4 t (xx, t M)=t: 3. C alculate the sensitivity coefficient the previous algorithm SEC . 6.3 5 Sum (inverse) K X . ,I calculation 2 I I = W.(dl +a.J2) E xtra inverse calculations 2 K +2 2K + 1 in exactly the same way as for (6.3.13) (6.3.14) 4. Calculate the heat flux t hat exactly matches the experimental d ata by using the Taylor series expansion, (6.3.8) 5. With the heat flux qM calculated from step 4, recalculate f l = wl(d l + ad2) a nd repeat the back substitution portion o f t he tridiagonal matrix algorithm. T~=fl (6.3.18) Figure 6.2 presents the results of Eq. (6.3.19) for the C DC 7600 computer. I f the temperature sensor is located near the active surface (small numerical value o f K), t hen the restructured algorithm is very efficient a nd is faster t han the originally proposed algorithm. The foregoing type o f minimization o f calculations for the I HCP was first developed by Blackwell.3 However, the original approach did n ot m ake explicit use o f t he sensitivity coefficient concept. I n real-world problems, there will be errors in the temperature measurements. I f t he temperature d ata a re matched exactly, the temperature errors produce heat flux errors. Consequently, exact matching o f t he temperature d ata 0 . 40r----.....---~--~--~--, ~I. 0.30 '-- I !' ~ .M.-_ ~ .. N ote t hat it is n ot possible t o calculate the temperature field from the Taylor series result, Eq. (6.3.15), because the sensitivity coefficients for depths greater than the temperature sensor depth were n ot calculated in step 2. An operational c ount d emonstrates that the restructured algorithm is m ore efficient than the one first proposed. Table 6.5 summarizes the operational count for the new algorithm. N ote t hat when steps 1 a nd 5 a re combined, they are equivalent t o o ne direct problem solution plus one additional evaluation of f l = wl(d l + ati2)' / / ~T For ~f.W!i~5 / 00 20 40 60 80 K, Node number of temperature sensor T he c omputation index for the restructured algorithm is M t - 8 0.10 K + 2+(2K + 1)(r",/tQ)+(td/tQ) 3 (N - 1)+ 5 (N - 1)(t",/tQ)+ N(td/tQ) FIGURE 6 .2 C omputation index for algorithm that exactly matches single temperature sensor data ( K = sensor node number, N = number o f nodes, K ~ N ) . (6.3. 19) 230 C HAP. 6 ONE-DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM is n ot recommended unless the temperature data are smoothed prior to calculating the heat flux and/or relatively large time steps are taken (a.lit/xi. > 1). The function specification method using future temperatures (Chapter 4) can be used with difference methods and allows smaller computational time steps to be used without encountering stability problems; this method is discussed in Sections 6.6 and 6.7. 6 .4 M ULTIPLE T EMPERATURE SENSORS. FUNCTION SPECIFICATION ( q=C). SINGLE FUTURE T IME STEP Multiple temperature sensors are recommended in order t o o btain as much experimental information as possible. These include sensors at different depths as well as several sensors at the same depth. F or multiple sensors, only one sensor can be matched exactly. An alternative is t o determine the value of qM, c onstant over the interval tM - 1 ~ t~ tM • t hat minimizes the least squares error for all sensors considered. I f there is a total of J sensors, then the least squares error is defined as SEC . 6.4 M ULTIPLE TEMPERATURE SENSORS I f J = 1, then Eq. (6.4.5) reduces to Eq. (6.3 .8). Again, note that q* = 0 is a valid assumption. The computation index in Figure 6.2 is very nearly valid for the multiple sensor case provided K is the node number of the sensor located the greatest distance below the active surface. The only additional calculation of the multiplesensor problem compared to the single-sensor problem (with the sensor located at node K ) is the evaluation of Eq. (6.4 .5) instead of Eq. (6.3 .8). This additional calculation is generally trivial in comparison to the total number of calculations. Hence, Figure 6.2 is a n approximate indicator for the multiple temperature sensor case. An alternative to the use of Eq. (6.4.5) for the multiple-sensor problem for a single future time step is to calculate the heat flux Zit .u t hat exactly matches the individual (kth) measured temperature. > t,u, a nd then take a weighted average o f all Zit.M'S . This average can be derived by using a Taylor series slightly different from Eq. (6.4.4). Expanding the temperature field for q u in terms of ilt.M T(xto t M , t "'_I, q U-I' q u)= T (xt, t u, t U-lo q U-I, Zit.u) _ ) oT(Xtotu.tu-l,qU-I,qu)1 + (q u-qt.u 0 qu 9. . -4.... (6.4.6) J s= L t =l [>t . M-T(x to t u ,tU - lo qU-loqU)]2 (6.4.1) where k is the sensor index and M is the time index in >t . M'* I t is i mportant to remember that k is the sensor index, not the finite difference node index. Differentiating Eq. (6.4.1) with respect to qM, gives oS - 0 = 0=2 qu J L [ >t.u-T(Xt, tM.tU - t =l 1• q U-loqU)][-X(xtotU,tU-1,qU-I)] (6.4.2) But exactly matching the kth sensor at t u gives T (xt. t u, t U-I' q U-lo Zit.u) = >t.M T (xt, t u, t U-I' q U-I' q",)= f r + (qu-q*)Xt . 1 (6.4 .4) Substituting Eq. (6.4.4) into Eq. (6.4.2), replacing q u by q u and solving for q u, gives J L X t.I(>t.u-fr) t q u = q* + c..=-=..I- . J - - - - OT(Xl' t u , t"' _ I' q U-l, q u )1 o = X(xl , t u, t U-I, q u-d=Xt.1 qM 9. . - 4. ... The minimum least squares condition, Eq. (6.4.2), becomes L Xf.l t =1 ° The double subscript notation f (x, . t ",)= Yo.", corresponds t o the classical (x. t) notation. - ----- (6.4.8) J L X t.l(qu-Zit.u)Xt . 1= 0 t =1 (6.4.9) Solving for q u, J q u= L WtZit .u t= I (6.4.10) where the heat flux weighting factors are related t o the sensitivity coefficients through Xl. I (6.4.5) (6.4.7) and (6.4.3) As in Section 6.3, a n arbitrary value o f heat flux q* is assumed and the Taylor series expansion given by Eq. (6.3.5) is applied: 231 J (6.4.11) L x li t =1 Equation (6.4.10) provides an alternative interpretation to this solution of the multiple sensor, single future time step I HCP. T he procedure is to treat 232 C HAP.6 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM each sensor independently (over a single time step) and determine the heat flux qt.M that exactly matches the single data point l't.M ' This. heat flu~ is related t o the assumed arbitrary heat flux q* through the Taylor !.enes expansion ~ _ qt,M=q • + l 't,M-f: , k -1,2, . .. , J After evaluating Eq. (6.4.12), these qt.M heat flux values are weighted (averaged) according to Eq. (6.4.10). F rom the discussion of sensitivity coefficients in Section 1.6, the sensors closest to the active surface will have the largest sensitivity coefficients and hence will automatically be weighted more heavily. F or a problem with constant thermal properties, the sensitivity coefficients X t,l a nd hence heat flux weighting coefficients Wt need be calculated only once. I f the thermal properties are treated as being quasi-linear in temperature, both Xt • 1 and Wt must be recalculated a t each problem time step. 6.1 . Using the information in Table 1.1 for the response o f a p lanar slab with a n insulated inactive surface to a step in heat ftux, calculate the heat ftux weighting coefficients o f E q. (6.4.11) for two sensors a t dimensionless depths o f x t!L=O a nd X2/L=0.s a nd for two dimensionless time steps o f d t+ = alit/L2 = 0.05 a nd 0 .5. N ote t hat T able 1.1 is also shown in graphical form in Fig. 1.7. E XAMPLE Solution. Although Table 1.1 c ontains the dimensionless temperature response, the dimensionless sensitivity coefficients iJT/iJ(qL/k) a re identical t o T+. a lit 0 =0.5 X t,l = 0.252313 X i,I =0.458333 W I =0.7671 w2=0.2329 response for t > tM' One approach that makes use of all tempera.tur~ informat.ion simultaneously is the whole domain estimation procedure which IS the subject of the next section. 6.5 W HOLE D OMAIN E STIMATION W ITH DIFFERENCE M ETHODS The whole domain procedure has been discussed in Sections 4 .4.2. This procedure is most applicable when the physics o f the problem suggests that the heat flux variation with time is o f a particular analytical form, for example, forms such as given by Eqs. (4.4.1)- (4.4.5). This list is by no means exhaustive. A definite advantage o f the whole domain estimation procedure is t hat stability problems are minimized. However, more computations are required. . The whole domain procedure will be introduced first by means o f a speCific example and then will be generalized, Suppose a parabola adequately represents the time variation o f the heat flux, (6.5.1) The parameters P I' Pz, P3 are to be estimated so that the computed temperatures agree (within certain limits) with the experimentally measured temperatures. The extent of the fit o r agreement is determined by the ordinary least squares criteria. I f Yi is the experimentally measured temperature and 7; is the predicted temperature, both at time t ; a nd a t the same location, then the least squares error is given by X t,l =0.831876 X i,I =0.0153659 W I = 0.9963 (x/L=O) W 2 = 0 .0037 ( x/L=0.5) 233 W HOLE D OMAIN E STIMATION WITH DIFFERENCE METHODS (6412) .. X t,l + a lit d t = -2 = 0.05 L SEC. 6.5 r S= L (Yi _7;)2 (6.5.2) i =1 o The predicted temperatures depend implicitly on the parameters (P I , P2, P3) 7;=7;(PI,P2,P3) In comparing the W t'S for the foregoing example, several important conclusions can be drawn. F or at+ = 0.05 the sensor at x+ = 0.5 yields very little information because its sensitivity coefficient is small relative to that for the surface-mounted sensor. F or at+ = 0.5, the time step is sufficiently large that sensor 2 has time to respond; however, sensor 1 is still weighted more heavily because o fthe square of the sensitivity coefficients. Because temperature measurements invariably contain errors, the heat flux components determined by exactly matching a single sensor contain amplified errors. Even for the multiple sensors considered in Example 6.1, sensors located far apart tend to olTer very little smoothing o f the heat flux; consequently, other approaches should be considered. Smoothing o f the temperature data prior to the solution o f the I HCP certainly reduces the oscillations in the heat flux; however, this approach does not make use of the fact that a pulse of magnitude qM over the interval tM-I~t~tM has an influence on the sensor (6.5.3) The object is to choose (Ph Pl, P3) such that the ordinary least squares error S is minimized. This condition is determined by oS r 07; . (6.5.4) L (Yi-1;) U!>pj = 0, J = 1 ,2,3 oPj ; =1 Assuming that the computed temperature is continuous in the parameters (PI' P l, P3), the temperature can b e expanded in a Taylor series about trial values o fthe parameters. I f (PI. Pl, pj) represent the trial values of the parameters and (P~+ I, Pl+ I, pj+ I) represent new o r improved estimates o f the parameters, then - = -2 t ;+I=t ; +(Pi+I-PD07;1 OPI ,=,- +(P'2+ I _Pl)o7;1 OP2 ,=,- +( Pj+l_p;)o7;1,=,OP3 (6.5.5) _ _ _ __ 0 _. _ _ _ _ _ _ • _ _ _ _ _ SEC. 6.5 234 C HAP.6 W HOLE D OMAIN E STIMATION WITH DIFFERENCE METHODS 235 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM sensitivity-coefficient (11) differential equation is Note that higher-order terms are ignored in the Taylor series expansion; however, the expansion is exact for linear problems. Substituting Eq. (6.5.5) into each o f the three equations represented by (6.5.4) and simplifying, ±OoT; OoT; ±OoP2 OoT; ±aopJ OoT; • IW'-Pi T; T; PI PI PI PI ±OoT; OoT; ±OoP2 0oT; ±ooPJ OoP2 P't Pi T; T; T; PI P2 P2 ±OoT; ooPJ ±OoP2 0oT; T; T; Pt PJ j= I i= I ok OTI - 1I1"T" u u x x =o 011 / 1 - k" u x x =o = t j - I j-l,2,3 , ·_ (6.5.8) i= I I- i =' j =. i= I i =' _ 1I i= I r oT. r oT. i =1 OP2 r oT. L (~-1;)-' i~1 OPI L (~-1;)-' OkOT I _k01l/1 = 0 loT ox x =L ox x =L j =1,2, 3 v (6.5.10) (6.5.11) (6.5.6) L (~-1;)-' OP3 i =1 The superscripts v o n the coefficient matrix and the right-hand-side vector indicate that all terms should be evaluated a t P= p'. Equation (6.5.6) is a linear system of three algebraic equations with three unknowns; P't I - Pi is treated as an unknown correction factor. F or a linear problem, the sensitivity coefficients are independent o f the parameters, the coefficient matrix is in turn independent of parameters, and only one step (iteration) is required. Even though the heat flux in this example is linear in the parameters (PI, P2, PJ), temperaturedependent thermal properties make the problem nonlinear and require iteration. F or the nonlinear problem, one approach to the sensitivity coefficient calculation is to use differences. F or example, (6.5.7) where e is a small number, say 0.001. The evaluation of each sensitivity coefficient requires the solution o f two direct problems but one is the standard condition; for the foregoing example with three parameters, a total of four direct problems must be solved for each iteration in Eq. (6.5.6). Each iteration will also require the simultaneous solution o f three linear algebraic equations, Eq. (6.5.6). (See the comments below (4.5.19).) At first glance, it might appear that some computational savings are possible if the sensitivity coefficients are calculated directly from a differential equation similar to Eq. (6.2.7). However, this is not the case for the nonlinear temperaturedependent thermal property case. F or the whole domain problem, it may not be possible to locally linearize the problem by evaluating the thermal properties at the previous time step as it was for the sequential procedure o f Section 6.3. Considering Eq. (6.2.7) as the basic model for the temperature field, the This sensitivity-coefficient differential equation is linear in 111 if the temperature field is known because terms like k, pc, ok/oT, and O(pc)/oT can be treated as known functions of position. However, the differential equation for 111 is no longer the same as that for T. Although the parameters (PI' P2, P3) d o not explicitly appear in Eqs. (6.5.8) - (6.5.11), they appear implicitly in the predicted temperature field T(P., P l' PJ). I f the heat flux model is nonlinear in the parameters, then the parameters appear explicitly in the sensitivity coefficient boundary condition. F or example, using Eq. (4.4.4) as the heat flux model, (4.4.4) then the boundary condition for 11, is 11 -(111 oT Ox x=o - (k 0ox ')1 x=o =Plte- fl ,' ok T)I o (6.5.12) which is nonlinear in the parameters. The foregoing procedure can be readily extended to heat flux models with a large number o f parameters. In such a case, it is beneficial to use matrix notation; additional details on matrix calculus can be found in Chapter 6 of Beck a nd Arnold. 4 Let Y be the vector of the discrete temperature measurements and T be the vector of the corresponding computed temperatures at the given sensor location, (6.5.13) The subscripts on T and Y denote time and there are r discrete times. The matrix form o f the least squares function, Eq. (6.5.2), is (6.5.14) 236 C HAP 6. O NE-DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM T he vector o f p parameters, I, is SEC . 6.6 237 SINGLE T EMPERATURE SENSOR where (6.5.24) (6.5.15) (with p <r) a nd it is t o be determined so t hat S is a minimum. This procedure will involve the matrix derivative o f Eq. (6.5.14); this derivative is defined as a a pi VI=- a ap2 (6.5.16) 6.6 SINGLE T EMPERATURE SENSOR, F UNCTION S PECIFICATION ( q=C)" F UTURE T IME STEPS a app Differentiating Eq. (6.5.14) with respect to F or a linear problem, X (I) is i ndependent o f IJ a nd only one iteration o f Eq. (6.5.22) is required. I f t he sensitivity coefficients are calculated using differences as in E q. (6 .5.7), a p -parameter problem will require p + 1 direct problem solutions for each iteration. F or p large and approaching r, as was the case for the sequential method o f Section 6.3, the number o f calculations required for the whole domain estimation becomes very large, even for the linear problem. Consequently, the whole domain estimation procedure is n ot recommended except when the functional form o f t he heat flux model is strongly suggested by the physics o f t he problem a nd the number o f p arameters is reasonably small. I, VIS = - 2(VI TT)(y - T)=O (6.5.17) T he sensitivity coefficient matrix is defined as X (I) = (V, TT)T (6.5.18) or, I -- aTI a T ap, ap2 a TI app aT, oT, -api ap2 aT, ap p (6.5.19) X= tr Substituting Eq. (6.5.18) i nto Eq. (6.5.17), X T(Y-T)=O T he discussion in Section 1.6 p ointed o ut t hat the sensitivity coefficient associated with a heat flux pulse o f finite duration can be nonzero after the heat flux returns to zero. Consequently, q(t)=q"" over the interval t"'_I~I~t", a nd zero otherwise, will influence sensor measurements Y"" Y",+ I ' . ... Therefore, powerful procedures for estimating q", s hould use future temperature measurements. BeckS was the first "to recognize the importance o f future temperature information and apply it t o the Duhamel's theorem solution o f the I HCP. Beck and Wolf' were the first t o c ombine difference methods a nd future temPeratures for the solution o f the I HCP. Beck et al. 6 applied sensitivity coefficient concepts and substantially reduced the number o f c omputations required for difference methods with future temperatures. This section considers a single temperature sensor for the function specification method t hat assumes a constant heat flux and an arbitrary number o f future time steps. Suppose the I HCP h as been solved up to time t ", - I; the estimated heat flux li", - I a nd the entire temperature field - I a re known. Next. the time is advanced one step to t ", a nd a n estimate o f q", is calculated. A single temperature sensor is located a t a d epth x . below the active surface. A temporary a ssumption is m ade that the heat flux is c onstant over r future times. The I HCP is shown schematically in Figure 6.3. An estimate is s ought o f t he value o f q"" c onstant over r future time steps. t hat minimizes the least squares e rror between the computed and measured sensor temperatures. The least squares function is (6.5.20) T he matrix form o f t he Taylor series expansion o f the temperatures in terms of parameter space is (6.5.21) E quation (6.5.21) is a matrix form o f Eq. (6.5.5). Retaining only the linear terms in Eq. (6.5.21) a nd substituting into Eq. (6.5.20) gives the normal equations (6.5.22) which is a matrix form o fEq. (6.5.6). T he solution can be symbolically written as , s= iI [ Y(x., t "'+i- d- T(x•• t "'+i-I. t "'_I' q "'_I' q",)]2 =l Differentiating S with respect to q",. replacing q", by l i",. a nd setting as/aq", equal to zero yields, as o=~= - 2 u q", (6.5.23) (6.6.1) r I [ Y(x •• t "'+i-1) i =1 - T(x., t "'+i-I' t "'_I' q "'_I, l i",)]Z(x., t "'+I_I, t "'_I' q "'_I) (6.6.2) 238 C HAP.6 O NE·DIMENSIONAL INVERSE HEAT C ONDUCTION P ROBLEM S EC.6.6 239 SINGLE T EMPERATURE S ENSOR Substituting Eq. (6.6 .5) into Eq. (6.6.2) and solving for qM gives. q (t) , T (t) L (¥t.M+I - 1q M=qt, +1= I f:'+I-' ) Z•. I (6.6.6) , LI (Z.Y j= E quation (6.6.6) can be written in terms of a gain coefficient as was Eq. (5.3.4a) (6.6.7a) I i M -IM M +2 T M +r-l M+ 1 M +r-2 Time 11 I I where M -IM M +2 M +r-l M+ 1 M +r-2 Time Z •. i K •.1 =-,-- F IGURE 1 .3 Pictorial representation of I HCP using r future temperatures. L (6.6.7b) Z l.j j= I where the sensitivity coefficient is defined by Z . '.1 = Z(xto I M +.-I. tM -I, qM -I )-_aT(XtotM+i-l.tM-I.qM-I.qM) aqM (6 .6.3) Since qM is assumed constant over r future time steps, the Z sensitivity coefficient is used instead o f X ; see Figure 6.1 for a comparison of the two sensitivity coefficients. Applying a Taylor series expansion similar t o Eq. (6.3.4) gives T (Xl. t M+i - l • ' M- I • q M-I' q M)= T (x • • t M+/- I , t M - I , q M-I. q *) + (qM - q*)Z(x• • tM + 1- 1, t M - I , q M-.) (6.6.4) Figure 6.4 presents a sketch o f the Taylor series expansion. In terms of index notation, Eq. (6.6.4) c an be written as T M+i-1 ~fM+I-1 (A * )Z • -. + qM-q '.10 .1 1 = . 2. ... . r (6.6.5) I f thermal properties are independent of temperature. then the unit step function solution ljJ used in Chapter 5 is identical to the sensitivity coefficient Z solution of this chapter and the Z values need t o be calculated only once for a given problem. If the thermal property variation with temperature is treated in a quasi-linear manner, then Z =1= ljJ a nd Z m ust be recalculated for each value o f qAl' T he subscript k in Eq. (6.6.4) is kept as a reminder that the sensitivity coefficients depend o n t he depth x . o f the sensor below the active surface; in the next section where the analysis is extended t o include multiple sensors, the subscript k is even more important. I t is necessary to remember that the consta{lt heat ftux assumption indicated in Figures 6.3 and 6.4 is a temporary assumption. The computed heat ftux qM is retained only over the time interval tM - 1 ~ t ~ t M . F or each subsequent time interval, a new heat ftux is calculated. The detailed computational aspects of the algorithm given by Eq. (6.6.7) for the constant heat ftux assumption. single temperature sensor, and r future times are summarized below: 1. Assume q (t)=q* over tM - 1~t~tM+'-1 and perform the forward elimination portion of the TDMA. Starting with the last node N, T *M+r-l • 1 w N=b ' N eN=wNcN, f j=wj(dj+aJj+.), f N=wNdN e j=wjcj. j =N-I.N-2 . .... 2 f l = W.(dl + atf2). d l =dl(q*) M -IM t M+2 ~r-l M+ 1 M +r-2 Time FIGURE 1 .4 M -IM t M+2 I tM+,-l M+ 1 M +r-2 Time Reference condition in the calculation of q ., by differences. Note that all a j. bj , cj • ej ' a nd Wj are fixed during the calculation of both f 'j'+i - I Z j . M + 1_ I . i = 0, 1, . .. • r - 1 because all thermal properties are held fixed a t ~nd T'j'. 240 C HAP.6 O NE-OIMENSIONAL I NVERSE H EAT C ONDUCTION P ROBLEM OM ° ' . .. 2. C aIcu Iate T j ' TOM+I , . .• , T M+, _ I usmg the back substitutIOn portion j j o ftheTDMA, o T M +i-I_ f I.i I - M+ t M+I-I-e) t ) -1 1- 1+ J'),i' I ) - j =2,3, . .. , N} Repeat for ; =O,I, . . . , r-l T he ; subscript on h,i is a reminder that this term must be recalculated for each value o f future time (;). However, the e/s are identical to those o f step 1 and do not change for ; = I, 2, . .. , r. 3. Calculate the step function sensitivity coefficients ( Z) using the difference equations developed from Eq. (6.2.7). The TDMA coefficients a j' bj , Cj ' e j' and W j have already been calculated in step 1. Therefore only dj a nd j j need be recalculated and the back substitution must be performed. Z I,M+I-I } Repeat for = f t.i Z j,M+i-l=ejZj-I,M+i-l+h,i> j =2,3, . .. , N S EC.6.7 6.7 M ULTIPLE T EMPERATURE S ENSORS, F UNCTION S PECIFICATION ( q=C), , F UTURE T IME STEPS Multiple temperature sensors were considered in Section 6.4 for a single future time step. In this section, the analysis is extended to include an arbitrary number of future time steps with the temporary assumption of constant heat flux. T he least squares error function must be modified to include a summation over the number of sensors. F or J sensors and r future times, define a sum of squares function by J , t= I i= I LL S= [ y",M+j-I-T(xt,t M+i - I ,IM-I,QM-I,qM)]2 oS N ote that ! he dj a nd h i terms for the calculation of Z j,M+i - 1 are different from those for T f + i - I. Also, the boundary conditions o n Z are different from those on 4. Calculate qM from (6.6.7a) t. , (6.7.1) The value of qM, constant over r future time steps, that minimizes S is sought, J, LL , ,-=0=2 uqM t =1 ; =O,I, . .. , r-1 241 M ULTIPLE T EMPERATURE S ENSORS [ T(Xl>tM+j-l,tM-I,qM-I,qM)-y",M+j-I]Zt,j (6.7.2) j =1 where the step function sensitivity coefficient is defined by - Z( Zt,j- Xl> t M+j - l , t M-I' ) q M-I o T(Xt,tM+j-l,tM-I,qM-I,qM) OqM ( 673 .. ) Expanding the temperature field in a Taylor series about an assumed heat flux " * '" q M=qM+ I... (Y.t ,M+I-I- TOM+i-I)K t ,i t q*, i= I 5. Knowing qM, reevaluate d l a nd I I a nd repeat the back substitution portion o f T DMA for t (x, tM, tM - I ' q M - I , qM) T :+j-l=t:+J-I+(qM-q*)Zt,j, j =I,2, . .. , r-1 (6,7 .4) substituting Eq. (6.7.4) into Eq. (6.7.2) a nd solving for the desired heat flux, t~=fl t f = e)Tf-1 + f i' j =2, 3, . . . , N (6.7.5.) EXAMPLE 6.2. Using the information in Table 1.1, calculate the K factors o f Eq. (6.6.7) for r =2,&I+ = (J.!u/L2 =0.05, a nd a temperature sensorlocated a t x+ = x/L=0.5 . Repeat the calculation for x+ = x/L= 1.0. Equation (6.7.5) can be written more compactly as J So/ulion. &1+ = 0.05, x+ = 0 .5 Z . = 0 .0153659 Z l =0.0593109 0.0153659 K - ------~--------1 - (0.01 5 365W + (0.0593109)2 = 4.0933 K 0.0593109 2 (0.01 5 365W +(0.0593109)2 =15.80 &1+ =0.05, q M=q*+ x+ = 1.0 , LL jl Kt,j(Y,..M+j - l-t:+ - ) (6.7.6a) t= I j= I Z . =0.0002693 Z2 =0.0078853 where K_ 0 .0002693 • - (0.0002693)2 + (0.0078853)2 (6.7.6b) =4.3261 0.0078853 K - ~----~------~ 2 - (0 .0002693)2 + (0.00785W EXAMPLE 6.3. Repeat Example 6.2 with the modification that two sensors ( J = 2) a re located at x+ = 0.5 a nd x+ = 1.0. F rom Eq. (6 .7.6b), = 126.7 The x+ = 1.0 a nd 1+ = 0.05 case is also included in Table 4.1. o Z •. i K •. i = Z2 Z2 Z2 Z2 •.• + 2.• + •. 2 + 2.2 b Zt,j SEC. 6 .9 CHAP. 6 242 Note that it is necessary to use the pulse sensitivity coefficient X .,I_ j + 1 = because the heat flux is allowed to vary in steps over the r future time steps. Substituting Eq. (6.8.1) i nto Eq. (6.8 .3) gives Solution. F rom E xample 6.2. A Z.,I _ j + I Z I.I = 0.01537 Z2.1 = 0 .0002693 ZI.2=0.05931 b =262.0 Z 2.2 = 0.007885 K I.I = 4 .027 K 1 • 2 =15 .54 K 2. 1= 0.07057 T :'+ j -I = t :,+ j -I + q",I/I.,j- qM-I I /I.,j-l - q*Z.,i K 2 • 2 =2.066 L i i 1 o T he presence o f a second sensor a t x + = 1 does not add significant information as evidenced by the fact t hat K 1 .1 a nd K I . 2 for this problem are very close to KI a nd K2 for Example 6.2 (for x + =0.5). T his is because o f t he relative size o f the sensitivity coefficients for x + =0.5 a nd 1.0. A Z• . I _j+I=Z•. i > L j =1 j AZ•. I _ j + l =L Z .,j=I/I.,j (6.8.5) j =1 Substituting Eq. (6.8.4) i nto Eq. (6.8.2) yields s= t ( t:'+I-I- Yt."'+I-1 + q",I/I.,I-qM_II/I•. I _I-q*Z.,1)2 (6.8.6) i =1 The heat flux component q", _ I is treated as a quantity known from the previous time step. T he value o f qM t hat minimizes S is found t o be Repeat Example 6.3 for t wo identical sensors. b oth a t x+ = 0.5. Solution. (6.8.4) where the following relationships are used: j= I EXAMPLE 6 .4. 243 SECOND-ORDER SEQUENTIAL REGULARIZATION M ETHODS O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM b = 133 .2 K 1 • 1 = K 2 • 1 = 2.047 o K 1 • 2 =K 2 • 2 =7.900 6 .8 SINGLE T EMPERATURE SENSOR. FUNCTION SPECIFICATION, L INEAR HEAT FLUX (CONNECTED S EGMENTS) (6.8.1) U p t o this point. all the function specification methods in Chapter 6 have utilized the (temporary) assumption o f a c onstant heat flux. O ther assumed functional forms were considered in Chapters 4 a nd 5. All o fthe h eat flux variations utilized with the convolution integral methods in Chapter 5 c an also be utilized with the difference methods o f this chapter. An example is the connected linear heat flux functional form shown in Fig. 4.8. T his approach is based on a straight line e~tra~olation in time through the two points q"'_1 a nd q",. F or equally spaced time mcrements. the heat flux variation can be expressed as q "'+I-I=(l-i)qM-I+iq",. i =O.I • ...• r (6.8.1) N ote t hat Eq. (6.8.1) c ontains only a single unknown. qM' Additional details can be found in Section T he following least squares e rror is minimized , s= L [ Y(x •• t M+1_I)-T(x•• t "'+I_I. t "'-I' q M-I. q", •. ..• qM+,_1)]2 (6.8.2) 1 =1 A comparison of Eqs. (6.8.2) a nd (6.6.1) indicates t hat t he temperature field for the linear heat flux case depends on the r heat flux values q", • .. . • q M+,-I; however. there is only one independent heat flux (qM) since they are related through Eq. (6.8.1). Using the standard Taylor series expansion. the temperature field can be written as 1 T M+1-I-T"M+I-I • -. + " AZ. ,I-j+1 (q M+j-I-q *) ~a i -I (6.8 .3) A c omparison o f the constant heat flux result. Eq. (6.6.1). s hould be m ade with the linear heat flux result Eq. (6.8.1). Both are o f t he same general form b ut the coefficients are different a nd the latter depends explicitly on q", - I ' T he linear heat flux variation yields a more complicated equation t han does the constant heat flux result. Experience has shown that the linear heat flux algorithm requires larger time steps for stability t han does the constant heat flux a lgorithm for the same value o f r. 6.9 S ECOND-ORDER SEQUENTIAL REGULARIZATION M ETHODS The general aspects o f regularization methods were discussed in Section 4.5 a nd the convolution integral method was applied to first-order sequential regularization in Section 5.4.3. In this section. the equations for the secondorder sequential regularization will be developed for the case o f r = 3 a nd then generalized by using matrix notation for an arbitrary r. T he matrix derivation also applies t o o ther orders o f regularization. The second-order regularization simulates the continuous term r'"'' j,.,-1 d2 -I ( i)2 dt (6.9.1) dt This integral is used in the cubic interpolatory spline 7 a nd the minimization o f it yields the so-called minimum curvature property. 244 CHAP. 6 O NE-DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM Following the material in Section 4.5, the standard least squares function is augmented as follows: SEC . 6.9 2 45 SECOND-ORDER SEQUENTIAL REGULARIZATION M ETHODS 3 L [ f:'+I-I- Y"'+i-I 1 =3 r- 2 , S= L ( Y"'+I-I- T:'+I-I)2+CX 1 =1 L ( q"'+I_I-2q"'+I+q"'+I+d 2 (6.9.2) 1=1 i + L X 1,I-j+l(q",+j-l-q·)]X1,i-2+ CX (q",-2q"'+1 + q"'+2)=O (6.9.6c) J= 1 with r~ 3. N ote that cx is dimensional and has units of T/q2. I n order to keep the development simple, a value of r =3 is assumed; therefore, the second summation in Eq. (6.9.2) contains only one term. Differentiating Eq. (6.9.2) with respect t o the three unknown heat flux components q"" q", + 1 , a nd q", + 2 a nd setting the results equal t o zero gives, After straightforward but lengthy algebra, Eq. (6.9.6) can be written as C tl X t'iXl'i+ CX ) C tl X1,,+1 X t,i- 2CX ) C tl Xt ,i+2 Xl,i+ CX) C tl X 1,i X l,i+I-2CX) C tl X 1"Xt ,i+ 4CX ) C tl X t ,i+ I X 1,,-2CX) (.± (.± (.± 1 =1 X l,iXl ,i+2+CX) p =l X1,iXl,i+I-2CX) . .= ) X1"X1'i+ CX) 3 L X 1,;(Y"'+I_I-f:'+i-l) 1= 1 2 (6.9.3c) L • M +· X1,;(Y"'+i - T t (6.9.7) ') i =1 1 T he lower limit o n the summation does not always start with unity because a sensitivity coefficient is identically zero for all times prior to the initiation of the particular heat flux component. Sensitivity coefficients are defined through i JT:'+I-1 X1•1- j+ 1 iJq"'+j_1 (6.9.4) The usual Taylor series expansion is written as 1 T I"'+I-I = T IO"'+I-I+'" X t .l-j+1 (q "'+j-I-q . ) L... (6.9.5) j= 1 Substituting Eqs. (6.9.4) a nd (6.9.5) into Eq. (6.9.3), one obtains 3 ' " [TO " '+/-1 L... t - i-I Y"'+/-1 . L X t ,l-j+l(q"'+j-I-q·)]X1 ,I+CX(q",-2q"'+1 + q"'+2)=O i -=1 The symmetrical coefficient matrix affords some computational savings. The regularization parameter cx improves the diagonal dominance o f the coefficient matrix and reduces the sensitivity to measurement errors. Since the equations are linear in q j-q., the value o f q . is completely arbitrary. If q ·=O is used, then f is replaced by Tlq=o. I n the sequential procedure, generally only q", found from Eq. (6.9.7) is retained. The procedure is then repeated for the next time step. F or the second-order regularization, the smallest allowable value o f future times is r = 3. Even for this simplest case, the algebra is rather lengthy. In order to extend this procedure to arbitrary r, matrix methods are used. The matrix analog of Eq. (6.9.2) is S =(Y - T)T(y - T)+cx(H 2 Q)T(H 2 q) 1 + 'L" Xt,i (Y." '+/+1- TO"'+i+l) ... t (6.9.6a) j= 1 The Y and T matrices have already been defined in Eq. (6.5.13). Since the procedure being developed is sequential, the heat flux vector contains only r elements, 3 ' " [TO I" '+/-1 L... 1 =2 - Y"'+/-1 . Q= L1 X 1,I-J+I(q"'+j-I-q·)]X1,i-I-2cx(q",-2q"'+1 + q"'+2)=O j= - q", q"'+1 : [ q "'+r-2 1 + (6.9.6b) (6.9.8) q "'+r-\ ] (6.9.9) 246 C HAP.6 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM H 2= 0 0 0 0 0 1 -2 1 0 00 10 -2 0 0 0 0 1 r -2 0 0 0 0 -2 1 00 00 01 00 00 . T =T+X(q-q*) (6.9.10) 1 (6.9.20) (6.9 .11) for T, and q* is a vector of r elements all equal to q*, and X is the sensitivity coefficient matrix, Xl,r X l,l 0 0 0 X1,r - 1 0 0 0 (6.9.12) 0 X1. 2 Xl,. Note that X is lower triangular and each element is defined by Eq. (6.9.4) . The matrix derivative of S with respect to the heat flux vector q is given by VqS = 2[Vq(Y - Tn(y - T) + 2exVq(H 2q)T(H 2q) ==0 Equation (6.9.19) is a matrix form o fEq. (6.9.7) a nd similarities can be observed. t is replaced by Tlq=o. With matrix manipulation subroutines available, Eq. (6.9.19) is much preferred over a generalization of Eq. (6.9 .7) t o arbitrary r. The solution to Eq. (6 ,9.19) can be written symbolically as I f q*=O, then . where the elements of T are defined in a m anner completely analogous to that 0 0 247 (6.9.19) F or r = 3, only the first row o f H2 is nonzero. The matrix equivalent of the Taylor series expansion, Eq. (6.8.3), is Xu 0 Xl.2 Xl,l X = X u Xl ,2 SPACE M ARCHING TECHNIQUES The matrix product H21 is the vector, the elements of which are the sum of the elements of each corresponding row o f H 2' T he elements of H 21 indiv,idually sum to zero, Equation (6.9.16) can be written as the matrix normal equation By analogy with Eq. (4.5.16c), the H2 matrix is -2 SEC. 6.10 (6.9 . 13) F rom Eq. (6 .9.11), I t is c ommon in the sequential procedure t o retain only the qM c omponent o f q and repeat the process for each successive time step. In the linear version o f the sequential regularization, the thermal properties are maintained constant at values corresponding to the temperature field at time t M _ I ' I f t he thermal properties are independent of temperature over the range o f interest, the whole domain regularization equations are identical to those o f Eq. (6.9.20) provided r is the total number of heat flux components being estimated. Since the sequential regularization method requires the solution of a system of r algebraic equations for each value of qM, it is computationally more costly than the other sequential methods presented in this chapter but for small values of r such as r~ 5 the increase in computation may not be large. The individual elements of X are calculated by means similar t o those discussed in Sections 6.2 a nd 6 .3. Namely, X 1,j-l+ 1 is determined from the difference solution of a direct problem in which the initial profile is zero, the heat flux is unity during the time interval tM+i - 2~t~tM+i-l a nd zero thereafter, and the inactive boundary is perfectly insulated. F or a discussion of the choice o fthe regularization parameter ex, see Chapter 5. Although Eq. (6.9.19) is derived for a second-order regularization, it also applies to other orders. F or a zeroth-order regularization, for example, H2 is replaced by Ho which is developed in Section 4.5 . (6.9. 14) Also, Vq(H 2q)T = vq(qTHIl=HI (6.9 . 15) Substituting Eqs. (6 .9.11), (6.9 . 14), and (6.9.15) into Eq . (6.9 . 13) and simplifying: XTX(q - q*)+exHIH 2q =X1"(y - t) (6 .9.16) Equation (6.9.16) can be put in a different form by noting that exH~ ' H2q*=0 (6.9.17) I f I is a vector with all elements equal unity, then for q* = q*1, H 2q*=q*H 21 = 0 (6.9. 18) 6.10 SPACE M ARCHING T ECHNIQUES FOR O NE-DIMENSIONAL P ROBLEMS T he difference methods presented in the previous sections can be termed "time marching" because starting with the initial temperature profile, the entire temperature field is calculated at the next time step. O ne marches ahead in time. F or the inverse problem, two boundary conditions are known at the sensor location; namely, temperature, and heat flux. The knowledge of two conditions at the same spatial location suggests that the I HCP is related to an initial value problem instead of a boundary value problem. This section explores several space marching techniques which are applicable to "initial" value problems. 248 CHAP. 6 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM I n the analysis that follows, a planar geometry with constant thermal properties is assumed; these restrictions are for simplicity only. F or the IHCP, the geometry can be divided into inverse and direct regions, as shown in Figure 1.4. T he direct region is the classical boundary value problem; the experimental temperature measurements are the boundary conditions for the interface between the direct and inverse regions. The boundary conditions are known a t the remaining boundary. Within the direct region, it is possible to solve this well-posed problem. Once the temperature field is known for the direct region, the heat flux a t the sensor location can be determined by differentiating the temperature profile. In the discussions that follow, all temperatures and heat fluxes are assumed known within the direct region; hence, both temperature and heat flux are known at the inverse boundary. 6.10.1 A nalytical S olution SEC.6.10 Equation (6.10.3) demonstrates the significance of the dimensionless time step for inverse problems. .' I f an error is made in the measurement of Y , there wIll be errors m both I I q£,' and qO , I' A lower bound o f this error can be . estimate~ by re~lacing Y by Y, + r, YI a nd ignoring the error in qli,l' The resultmg error m qO,1 IS found t o be t:.t£=~6t/E2 r,QI = _1_ sinh ( _1_) b YI (6.10.1) I f the time derivative is replaced by a backwards difference, Eq. (6.10.1) can be converted from a partial differential equation t o a n ordinary differential equation (6.10.2) where the superscript i denotes time index. The boundary conditions for Eq. (6.10.2) a re I x =1i Equation (6.10.2) has the appearance of an initial value problem because the function and its first derivative are both specified at x = E. I fthe initial temperature profile (To) is uniform throughout the body, a very simple solution of Eq. (6.10.2) exists for i = l (t = 6t) T '(z)- To=(YI - Yo) cosh ( ~)+qli",J;Xi sinh ( _ z _ ) ~~6t , J;Xi where z =E-x. The unknown flux a t x =O, ( z=E), is found by differentiating this equation with respect to x t o o btain O,1 E qli" E QI =q-k-=-k- cosh ( E ) + (Y, - Yo) C A. ~~6t ~ +q~1 ~sinh(~)J An analytical solution of this equation is possible; however, due to the complexity o f the inhomogeneous terms, the algebra is more c0'!lplicated an~ solutions a t subsequent times even more so. Consequently, this approach IS not pursued further. The important points t o remember from this analysis are that for a given value of time, it is possible to march in space, and that the important length scale in the dimensionless time step 6 tE is the sensor distance below the heated surface, E. 6.10.2 M ethod o f D 'Souza D'Souza 8 ,9 used the pure implicit difference approximation t o Eq. (6.10.~); the time derivative was replaced by a backward difference a nd the spatIal T i(x=E)= Y; q£.i= - k dTi(x) dx ~ which is tabulated in Table 6.6; these results are consistent with earlier observations that small values of dimensionless time 6 tli can produce large errors in surface flux. The foregoing process can be repeated for t = 2 6t by solving 2 d T2(Z) _ _1_ T2(Z) = _ _1_ [ To+(Y, - Yo) cosh ( ~) dz 2 ~6t ~6t ~~6t Consider the constant property heat conduction equation for a planar geometry 02T(x, t) 1 o T(x, t) iJx2 =~ - o-t- 2 49 SPACE M ARCHING TECHNIQUES E C A. ~~6t . smh ( E ) (6.10.3) C A. ~~6t T ABLE 6 .6 D imensionless H eat F lux E rror f or S ingle T emperature E rror «SY, d "tE cSQd{)Y, 10 1 0.5 0.25 0.125 0.0625 0.03125 0.1017 1.175 2.737 7.254 23.843 109.160 809.618 250 C HAP. 6 O NE-DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM S EC.6 .10 261 S PACE MARCHING TECHNIQUES (6.10.4) T he heat flux qh is unknown a nd is t o be estimated while q~ is the heat flux leaving t he inverse region shown in Figure 1.4; this latter heat flux is e ntering the direct region a nd is found from calculations for the direct region. I n Eq. (6.1O.6b) t he temperatures T~-' a nd T~ a re replaced by l j_1 a nd l j, respectively. E quation (6.1O.6b) is solved for T~ - I for i = 1, 2, . . . , n t o g et (6.10.5) (6.10.7) where p =aAt/Ax 2 • U sing t he space-time grid shown in Figure 6.5, t he calculation procedure is s tarted a t n ode j = K a nd sequentially steps t hrough time i = 1, 2, . .. , n. Because t he t emperature field is known along b oth t he s patial grid lines K a nd K + 1, Eq. (6.10.5) is a n explicit relationship even t hough it was derived from t he fully implicit equations for direct problems. After t he t emperature is c alculated all along space line K - 1, t he procedure is repeated for K - 2, . .. , 1. A finite difference e quation for node j = 1 yields a n e quation for calculating t he 'surface heat flux. T he D 'Souza p rocedure is intended t o b e implemented o n a c omputer as just outlined. Insight c an b e gained, however, by obtaining some algebraic expressions for a small n umber o f nodes, K , from the heated surface t o the sensor. U sing t he p ure implicit procedure, equations a re given for each node. F or n ode 1, t he h eated surface, the e quation a t time t i is T i _ Ti qi + k 2 I (6.10.6a) o Ax F or K > 2, Eq. (6.10.7) is i ntroduced into E q. (6.10.4) w ritten for j = K - 1 t o obtain a n expression for T~-2' i = 1, 2, . . .. F or K = 2, Eq. (6.10.7) gives expressions for T~ a nd (by replacing i by i -I) T~-' which a re i ntroduced into Eq. (6. 1O.6a) t o o btain d erivative by a central difference T~_,-2T}+ T}+, Ax 2 Solving for T~_, in Eq. (6.10.4) yields i T J -I = (2 + !) T i._TI J p J+ I _! P Ti. - I { j=K.K-l • ...• 2 J . 1 = 1• 2, ... , n T he e quation for a n i nterior n ode j is given by Eq. (6.10.4) a nd t he e quation for n ode K is (6.10.6b) (6.10.8a) i _q~E QE =T (6.10.8b) since Ax = E for this case. T he s ymbol V is the backward difference o perator, V lj= l j-lj_I' V llj= l j-2lj_, + l j-2 (6.10.8e) I f t he foregoing procedure is used for K =3, a nd t hus A x=E/2, t he heat flux expression is i V lj 3 V 2lj 1 V3lj i 1 VQ~ 1 V2Q~ Qo= AtE + 16 ( A't'd+ 128 (A't'E)3+QE+ 2 A't'E + 32 (V't'£l2 0 Direct region Inverse region 0 -1 • D"Souzal8.91 molecule • r -- ~ OJ • JI~---"-- OJ)" - -"\ 0 -2 " II E i= .) • ( i. I i -2 2 • j -l r I I I \ !I j • j +l Known conditions • • I nitial conditions j =O 1 • • 2 3 X -2 X -I Space index x X +l FIGURE 6 .6 Space-time grid for space marching methods. (6.10.8c,d) (6.10.9) Several c omments regarding Eqs. (6.10.8a) a nd (6.10.9) a re n ow made. First, this is a procedure t hat utilizes exact matching o f t he calculated temperature with the measured temperature a nd t hus is sensitive t o m easurement e rrors. Second, these equations a re n umerical analogs o f t he e xact solution which is given by t he B urggraf equation, Eq. (2.5.20). T he coefficients o f Eqs. (6.10.8a) a nd (6.10.9) a pproach t hose o f Eq. (2.5.20) as the n umber o f n odes increases; the first difference coefficients in b oth o f t he former equations a re the same as each o ther a nd t hose in Eq. (2.5.20) whereas the V2lj s tarts a t 0.25 [given by Eq. (6.10.8a)], decreases t o 0.1875 in Eq. (6.10.9), a nd a dditional nodes bring the coefficient closer t o t he exact value o f 0.1 667. Third. only past temperatures are used t o o btain t he p resent heat flux estimates. I t h as repeatedly been noted that future temperatures a re needed t o o btain accurate estimates with small time steps. F ourth, if a n a nalogous procedure were utilized starting with the explicit equation [i.e., left side o f Eq. (6.10.4) e valuated a t i -I], forward difference o perators w ould replace the backward difference operators o f Eqs. (6.10.8a) a nd (6.10.9), e nsuring t hat o nly future temperatures a re used. Fifth, as m ore I r 262 CHAP. 6 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM SEC.6.10 nodes are used, higher-order differences are needed, increasing sensitivity to temperature measurement errors. 6 .10.3 253 SPACE M ARCHING TECHNIQUES :0--~ n = 1 0· I 9· )( Weber worked with the hyperbolic form o f the constant property heat conduction equation molecule (R&Blll,12 7· ) 6• X 0 0 X X 0 )( 0 X )( X 0 0 X X ~ 3 ' --;'I I X \ I "5.I ~ _ __ ~ X r - _J X 1 L-(Ti-I_2Ti+I+Ti+I)+_(Ti.+I_Ti-I)=~(Ti _ 2Ti+Ti. ) .&t2 j j j 2.&t J j .&X2 j + I J rI )( X )( ~ __ .0 I -2 X (6.10.11a) i =00 j= 1 Solving for T }_I, Eq. (6.10.11 a) yields, Cl.&t ( 1= y -2 (.&X)2 IX .&t K Future times )( where y is a normalized thermal wave speed. Central differences in both time a nd space were used, p = .&x 2 ' 0 X _~'1 5 •_ II i= (6.10.1lb) 0 X )( E n -K+l=4(x _ __ X I (6.10.10) 0 0 x lO I o Computation 8· M ethod o f W eber r -- ! r-0=-'-~ x oo 3 2 S pace,j 0 Computation '..£;: molecule (Weber10) , __~ , ~ 0 oooo 7 o X 4 x 5 )( 6 K -l K 8 K +l FIGURE 6.6 Space-time grid for Weber'o and Raynaud and B ransier"· (node points), n = 10 (temperature measurements), 0 known. x calculated. 11 (6.10.11c,d) This algorithm is explicit; it can be started a t node j = K a nd time i = 2, 3, . .. , n - 1. T he value i = n is not permissible because T~+ I is n ot defined. Consequently the maximum time index for the grid locationj is one less than that for the grid location j + 1. This grid pattern is shown in Figure 6.6. I f the surface temperature is desired a t time M, then the temperature must be measured out to time i =M+K. T he algorithm given by Eq. (6.10.11b) has the appearance o f using K future times where K is the number of spatial grid points. Refining the grid has the effect of increasing the number offuture times. For the methods presented in Sections 6.3 a nd 6.4, the number o f future times is dependent o n the physics o f the problem and should be relatively independent o f the number o f spatial grid points. Weber indicated that the parameter (1 should be chosen small, without indicating how small. In fact, the algorithm works with (1=0. Algebraic expressions similar to those for the D'Souza procedure can also be found for Weber's. Following the same procedure as in Section 6.10.2 one can derive for two nodes; that is, K = 2, a nd with ( 1=0 (6.10.12b) where Q~ is a set equal t o zero. These equations are quite similar to Eqs. (6.10.8a) and (6.10.9), found for D'Souza's method, and comments one, two, and five made in connection with the latter equations also apply for Eq. (6.10.12). However, Eq. (6.10.12) has two advantages over Eqs. (6.10.8a) .and (6.10.9): (I) both past and future temperatures are used and (2) central differences are used that are much less sensitive t o measurement errors. 6.10.4 M ethod o f R aynaud a nd B ransier The method o f Raynaud and Bransier ll ,I2 consists of a two-ste~ marching procedure. The first consists of a space marching step followed by a tIme marching step; the final result is an average of the two procedures. F or the space march, the heat conduction equation is differenced as follows: .&x - . -. k -. -i k p c-(T,.+I_Tj)=--(Tj-Tj-d+ "A (T-i+t l - f i+l) j+ j (6.10.12a) ... methods; K =7 .&t J .&x (explicit) aX (implicit) ( 61013) .. 264 CHAP. 6 O NE-DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM ' Note t hat the conduction terms are mixed explicit/implicit. Solving for t i. I from Eq. (6.10.13), J- t~_1 = (1- !) t i+(1 +P t J.+ t~+ !) i J+ P J 1:...- I I the same as the Weber algorithm; see Figure 6.6 for details. The second part o f the algorithm approximates the heat conduction equation as follows: (implicit) (explicit) Again, the conduction terms are mixed explicit/implicit. Solving for from Eq. (6.10.15), ... T~:!:~=(I+~) T~+I+(I-~) T~-T~+I of this approach, see References 13 a nd 14. Some o f the basic features of this method can be understood by working with the combined form of Eq. (6.10.18). t he so-called leapfrog difference scheme was used, .1x .+ I .I e x. .. - (T'. - T' - )=-(T"+1-2T'+T"_d j 2M J . 1xJ J J TI+ I j -I (6.10.16) . 6.10.5 oq ox (6.10.18a) oT ox (6.10.20) ~----- -~ )( 9x )( x )( x )( 8x x )( )( x x G> G> 7 )( X )( X X x G> G> 6 )( )( X X X )( G> G> 5 )( )( )( X )( x I ., (6.10.21) 0 __ /i>, _~ _ )( I \@ G> 'C!)I X X X x 3 )( x x x x )( X X X I I I ~--, r--? x X X )( G> G> 0 G> I X '@ X ~ G> ,, __ J 1 )( G> , rx--' G> ' - 0' I i=OG> 0 I I " '- _ _ " ~--§._% 2 3 4 5 (6.10.18b) G> G> 0 K -l j= 1 and then converted them to finite difference equations individually; the two unknowns were T(x, t) a nd q(x, t). F or a discussion of the purported advantages . - Tj+1 )( and q =-k- ' +1 )( 2 )( equations oT ot 1 . 4 )( Hens~l a nd Hills 1 3 wrote the heat conduction equation as a pair o f first order p c-=- - I n = 1 0 )( E i= M ethod o f H ills a nd H ensel 1. 2p T r + 2Tj+ 2p T j . I ex _1 -1x (T.2 - Tj)=-(T jI + l - 2T jI + TjI ) M J .1x (6.10.17) Although there is a time marching component to the Raynaud and Bransier algorithm, the com~ined e~ect ~fth~ two steps has the same overall appearance as the Weber algorithm given 1 0 Figure 6.6 but it is even less sensitive to the ~eas~rement e rrors than Weber's method. The apparent number of future tImes IS the same as the number o f spatial grid points ( K). Raynaud and Bransier also considered temperature-dependent properties. =- See Figure 6.7 for the computational molecule for Eq. (6.10.20). (A "computational molecule" is a schematic representation o f a difference equation such as Eq. (6.10.20).) T he preceding algorithm is explicit and is identical to the Weber algorithm, Eq. (6.10.11 b), for u = 0. Although Eq. (6.10.20) is valid for j = 1, Hensel and Hills l3 chose to replace the centered time difference in Eq. (6.10.19) with a forward time difference T he Ti:!: ~ al.gorithm is explicit and time marching. Starting with j = 1, the ~pace mdex} steps through K - 1, K - 2, . .. , 1. T he final step o f the algorithm IS a n average o f t he foregoing two procedures. T 'j =Z(tlj +Tj) ~ •. (6.10.19) Although Richtmyer and Morton lS indicate that the foregoing difference technique is unconditionally unstable for direct problems, it can still be useful for inverse problems. I f Eq. (6.10.19) is solved for T~_I' Tj_1 (6.10.15) 266 SPACE M ARCHING TECHNIQUES (6.10.14) ~he e~plicit algorithm for t~_1 s tarts with n odej = K a nd steps through the time hnes 1 =2, 3, . .. , n - 1. h i general, the maximum value o f j is n +j - K which is . 1x... •. -k. • k .. pc - (T'.+ 1 _ T j)=- (Ti.+ 1 _ T i+ 1 )+_ ( T i T i) .1t J .1x J j -l.1 x j +l- j SEC.6.10 K K +l 6 7 8 Space F IGURE6.7 Space-time grid for Hills and Henscl'3.'4 method; K =7 (node points), n =10 (temperature measurementtimcs), 0 's known, x 's to be calculated. 2 68 C HAP.6 O NE-DIMENSIONAL I NVERSE HEAT C ONDUCTION P ROBLEM Solving for 1'/_ 1 , 1 1 ')-1 = ( 2 - 1) P 1 1 1') - 1')+ 1 1) +p 1'2 (6.10.22) The computational molecule for Eq. (6.10.22) is shown in Figure 6.7 . Equation (6.10.20) is n ot applicable for i = n. Instead of using the approach of Weber and generating the solution up to the diagonal line of Figure 6.6, Hensel and Hills l3 replaced the time derivative with a backward difference (6.10.23) Solving for 1 'j-1o 1 '- -1= ) (2 +p1) )- p) 1 1 '- 1 '--1 - 1'" )+1 (6.10.24) which is the same equation used by D'Souza. The computational molecule for Eq. (6.10.24) is shown in Figure 6.7. The depth of the discussion presented here does not do justice to the analysis presented by Hensel and Hills l3 and Hills and Hensel. I4 They extended the analysis to consider temperature-dependent thermal properties, nonuniform space-and-time grid, planar, cylindrical, and spherical geometries, and an estimate ofthe variance of the calculated surface heat ftux provided the variances of the temperature measurements and heat ftux a t the inverse boundary are known. Also, their algorithm has the capability of considering a digital filter of the measured temperatures. 6.10.6 C omparison w ith P rior M ethods A major weakness of most of the methods in Section 6.10 is the lack of a statistical basis. Another is that the methods do not readily extend to multiple interior sensors and to two- and three-dimensional cases. The function specification and regularization methods can have a statistical basis and the same formalism applied to one-dimensional problems can be directly extended to two- and three-dimensional cases. Moreover, quite different problems such as those involving integral equations o r ablation can be treated by the latter methods whereas the space marching techniques m ayor may not be appropriate. S EC.6.11 I t is appropriate to ask what difference method is the most accurate in the calculation of sensitivity coefficients. Although an exhaustive study has not been pedormed, Figure 6.8 provides some insight into this question; this figure considers a sensor at the midplane of a planar slab. The sensitivity coefficient is evaluated at a dimensionless time equal to the dimensionless time step of the difference method; that is, t + = At+. The backward difference ( 8= 1) method is superior for At+ near unity and larger whereas the central difference method is superior for At + < 0.1. As pointed out earlier, the sensitivity coefficient calculation time step A tx can be smaller than the problem time step At. This should produce more accurate sensitivity coefficients. Chapter 5 presented several example calculations related to a triangular heat ftux pulse. This same problem will be utilized here for a number of examples. The triangular heat ftuX pulse is maximum at t+ = 0.6 and then falls linearly to q+ = 0 a t t + = 1.2. The analytical solution of this problem was presented in Chapter 5 [see Eqs. (5.2.4) - (5.2.6). Table 5.1 tabulates the exact dimensionless temperature response (to six decimal places) for a sensor located 1 .o1.----------------__ 0.1 0.01 o 6.11 2 67 N UMERICAL CALCULATIONS N UMERICAL C ALCULATIONS - Analytical 0 6 = I , Lumped capacitance o 6 = 112, Lumped capacitance [ )'I/[)'Ix This section presents numerical calculations that bring out the salient features of several difference methods presented in this chapter. Most methods in this chapter require the calculation of various heat ftux sensitivity coefficients (iJ1'jiJq). The most powerful of these use differences. =1 o 0.OO~-::: . 0c:-1-~-----;:0";. 1---~-~---t1.0 [ )'I+ = a [).tlL2 F IGURE 6 .8 Calculation o f p lanar slab sensitivity coefficient by difference methods for x /elL =0.5. 268 C HAP. 8 O NE - DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM T ABLE 6 .7 T emperature R espon.e t o a T riangular H eat F lux P ulse. A dditive N ormal U ncorrelated R andom E rror. w ith S tandard D eviation o f 0 .0017 - 0.24 - 0.18 - 0.12 - 0.06 0.0 0.0600 0.1200 0.1800 0.2400 0.3000 0.3600 0.4200 0.4800 0.5400 0.6000 0.6600 0.7200 0.7800 0.8400 0.9000 0.9600 1.0200 1.0800 1.1400 1.2000 1.2600 1.3200 1.3800 1.4400 1.5000 0.00034 0.00281 0.00135 -0.00090 - 0.00020 - 0.004241 -0.000639 0.001307 0.007161 0.012874 0.021764 0.038954 0.054821 0.072952 0.098381 0.127722 0.155529 0.193275 0.223066 0.248541 0 .275499 0.293646 0.311805 0.330811 0.342215 0.348997 0.350227 0.355722 0.359318 0.361609 a t the insulated inactive surface Te by using a random numbe . mpe~ature measurement errors simulated deviation o f 11+ - 00017 r ge~erat~r With normal distribution and standard r- . are given I D Table 6 7Th ' approximately 5% o f the m ' d' ... IS standard deviation is . 0 aXlmum Imenslonless t . tnangular pulse. All difference I I . emperature nse for the ca cu atlons that follow are performed using S EC.6.11 259 NUMERICAL CALCULATIONS ( /= I and lumped capacitance. Similar results are expected for (J = 1/2 and/or other capacitance distribution schemes. Figure 6.9 compares results for exact and inexact sensor data, both with , = 3 for the function specification method, q = C. The results are quite good for the exact «(1; < 10 - 6) temperature data. However, the inclusion o f temperature errors causes scatter in the computed heat flux. The results for (1; =0.0017 are similar to those using Duhamel's theorem that are given in Figure 5.13; the (1; /11; ratio is 23.2 for Figure 6.9, which is close to the 25.4 value o f Figure 5.13. Figure 6.10 compares the function specification method ( q=C) for r =4 and 5 for data containing errors. Both values of r yield satisfactory values of heat flux considering that the d ata contain errors. As r increases, the ability to follow sharp changes in q is impaired, as evidenced by the calculations near 1+ = 0 ,0.6, a nd 1.2. T he optimal choices o f r a nd at are discussed in Section 5.6.1. The results shown in Figures 6.9 and 6.10 are nearly identical to those obtained using the convolution approach, Figures 5.12 a nd 5.13. This is expected because the I HCP algorithms are the same; the method o f solving the heat conduction equation and the evaluation of the sensitivity coefficients . are different. The random errors are also different. The same example problem was also solved using both first- and secondorder sequential regularization with CXI =CXl = 0.001 ; these results are shown in 0.6 0.5 0.4 q+ 0.3 0.2 o ,,~ < 10-6 '" ,,~ == 0.0017 0.1 '" 0 .0 - 0.1 - 0.3 " : == 23.2 for ,,~+ == 0.0017 "'''' "y 00 '" '" '" '" 0 .0 0.3 0.6 0.9 1.2 '" 1.5 t+ FIGURE 1 .9 Comparison or exact (11; < 1 0- 6 ) and inexact (11; =O .OOt 7) d ata ror runction specification method, q =C, 20 clements, r =3,lu/Al x =2, A I+ = 0.06,6= 1,Iumped capacitance, x /L= 1.0. SEC. 8.11 0.6 0.5 0.4 q+ 0.3 0.2 0.1 &6 & 0.0 0 0 6 - 0.1 - 0.3 0.0 0.3 0.6 0.9 1.2 1.5 t+ FI~URE 8.10 Function specification mcthod, q =C, 20 clemcnts, a ; =0.0017, A t/lux =2, At = 0.06,8= 1,Iumped capacitance, x /L= 1.0. 281 N UMERICAL C ALCULATIONS Figure 6.11. The first-order regularization appears superior to the secondorder by a considerable margin. The choice of (X was based on the results in Chapter 5 for the convolution integral approach. The analysis for the optimum choice of (X2 has not been performed. However, Figure 6.12 compares results for (X2 = 5 x 1 0- 3 with (X2 = 5 X 1 0- 4 and shows very little difference. Numerous examples were presented in Chapter 5 illustrating how a single temperature error affects the heat flux calculation at times both before and after the temperature error. Similar results are shown in Figure 6.13 and were determined by setting Y = 1.0 a nd all other temperature measurements identit cally zero. The second-order sequential regularization exhibits much greater swings than either the function specification (q = C) o r first-order sequential regularization. This partially explains why the second-order regularization performed poorly in Figure 6.11. The function specification method has the smallest swings. I t has been demonstrated that comparable results can be obtained for both function specification (q = C) and first-order sequential regularization methods. The presence of the regularization parameter (X requires that trial calculations be performed for each new class of problems solved by the regularization method. Both the function specification and first-order regularization methods have given reasonable results. 0.6 0.6 0.5 0.5 0.4 0.4 q+ 0.3 q+ 0.3 0.2 0.2 o 1st order, 0.1 6 6 a 0 0.0 - 0 III 6 0.1 2nd order, 0 2 = 0.00 1 o 8 & 0.0 112 = 5X 1 0-3 0 = 0.001 02 = 5X 1 0-4 • 8 III fI f! ~~.3--~~0~.0n-----rO~.3'-----~0.~6----~0~.9~----~I~.2~----~1.5 t+ FIGURE 8.11+ Comparison of first- and second-ordcr sequcntial rcgularization with «=0.001, 20 clemcnts, a r =0.0017, A t/Atx =2, At+ =0.06, r=5, 8 = 1,Iumped capacitance, x /L= 1.0. 280 0.3 0.6 0.9 1.2 1.5 t+ FIGURE 8.12 Influcnce of «2 on second-ordcr sequcntial regularization, 20 clcmcnts, 0.0017, A t/Atx =2, At+ =0.06, r =5, 8 = I, lumped capacitance, x /L= 1.0. a; = 2 62 CHAP. 6 ONE-DIMENSIONAL INVERSE HEAT CONDUCTION PROBLEM SEC. 6 .12 263 COMPUTER PROGRAMS T ABLE 6 .8 L ist o f C omputer P rograms f or t he I nverse H eat C onduction P roblem O q=C 1st order sequential regularization 10 Program Name Developer/ Organization C ONTA 6 o 2nd order sequential regularization Beck; Michigan State Univ.Sandia, Albuq. Muzzy, Avila, R oot;GESan Jose H CODE HEATINV v=[l.i= 1 . 10. ti = all other i INVERT 1.0 itJ.t Bell, Wardlaw; Naval Surf. Weap. Center Snider; E GG I daho O RMDIM Bass, Drake, O tt;ORNL O RINC O tt, Hedrick; O RNL S ODDIT Blackwell ; Sandia, Albuq. SMICC Hills, Hensel; New Mexico State Drake; ORNL Time index. i FIGURE 6 .13 Dimensionless heat ftux error due to a single measurement error at t + = At + = 0.06. r =S. 20 elements, At/At x =2.8= l,lumpedcapacitance, x /L= 1.0. 6.12 C OMPUTER P ROGRAMS A great many computer programs have been written for the inverse heat conduction problem. The purpose of this section is t o mention some of these. A review o f some computer programs is given by Beck. t 7 A list of some computer programs is given in Table 6.8. All o f the computer programs, except the last one listed, utilize some difference procedure for solving the transient heat conduction equation and thus can treat temperaturevariable properties. Each program has unique features and capabilities. Brief descriptions o f the programs listed in Table 6.8 are given next. C ONTA was developed by Beck 18 and was supported by both Michigan State University and Sandia National Laboratories. I t uses the function specification procedure with the temporary assumption o f a constant q. T he sensitivity coefficients are calculated in a very efficient manner.6 The program can treat multiple internal sensors in composite plates. H CODE was developed by General Electric 19 in connection with nuclear reactor safety analyses. I t can treat composite cylindrical rods. A modification o f the function specification method is used. HEATlNV 20 was developed a t the Naval Surface Weapons Center and is the only difference-based program listed that uses whole domain regularization. ... ARIES I HCP Algorithm Sequential function specification Modified sequentiai function specification Oth- + l st-order whole domain regularization Sequential function specification Sequential function specification Uses temperatures only a t time tAl Sequential function specification Space marching Oth-, 1st-, and 2nd-order whole domain regularization Difference Method Integral Method Crank-Nicolson No Crank-Nicolson No Crank-Nicolson No Generalized Crank-Nicolson No Finite element No Backward difference No Finite control volume, backward difference Yes No No Yes A combination o f zeroth- and first-order regularization is employed. A sequential regularization procedure, instead of the whole domain method, would be more computationally efficient. INVERT 1.0 was written by Snider 21 a t EG&G, Idaho, and uses the function specification method. I t can treat composite, cylindrical rods. O RMDIM was developed at O ak Ridge National Laboratory by Bass et at. 2 2 I t is a two-dimensional, finite element-based program. O RINC 23 was also developed a t O ak Ridge National Laboratory. I t is for composite cylinders. The fully implicit method o f solution of the difference equations is employed. The temperatures only up to time t are used to calculate the flux a t time t . 2 64 C HAP . 6 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM S ODDIT ~:s deve~oped a t Sand~a N ational .Laboratories, Albuquerque, by Blackwell. The difference equations are denved using the finite control v?lume. method. T he sequential function method is used a nd a rbitrary onedimenSional geometries are possible. i SMI<?C is a ~rog~am de~eloped by Hilts a nd Hensel 1 4 t hat uses space marchng a nd IS d escnbed I n Section 6.10.5. I t also provides estimates for the resolving power and accur~cy o f the calculated heat flux components. Hilts a nd Hensel are a~ New. MeXICO S tate University and the research was funded in part b Sandia NatIOnal Laboratories. Y ARIES 1 ' is t~e ~nlr progr~ listed in Table 6.8 t hat is based o n a n integral m~e~ ~nd t hus IS h mlted t o h near problems. The integral model does give the flex~bJllty t o tr~at a variety o f inverse problems, however. T he user has the OptlO~ o f choos~ng ~roth, first o r second (or some combination thereof) whole d om.un regulanzatlon. Copies o f some o f t he reports may be obtained from: National Technical Information Service U.S. D epartment of Commerce 5285 P ort Royal Road Springfield, VA 22161 a nd possibly from the authors. In some cases the reports contain listings of the programs. PROBLEMS 2 65 I I. Raynaud, M ., D etermination d u F lux Surfacique Traversant U ne P aroi Soumise a U n Incendie a u Moyen D 'Une M ethods D'inversion, Laboratoire D'Aerothermique G roupe ' Echanges Thermiques' Universite Pierre e t M arie C urie, Paris, F rance, August 1983. 12. Raynaud, M., a nd Bransier, J., A New Finite Difference Method for N on L inear Inverse Heat Conduction P roblem, t o be published in Numerical Heat Trans/er. 13. Hensel, E. C. a nd Hills, R. G., A Space M arching Finite Difference Algorithm for the O ne Dimensional Inverse Conduction H eat T ransfer P roblem, ASME P aper No. 84-HT-48, presented a t 22nd N ational H eat T ransfer Conference, Niagara Falls, NY, August 6 - 8,1984. 14. Hills, R. G . a nd Hensel, E. c., S MICC, t he Space M arching Inverse Conduction Code, S AND 84-1563, Sandia N ational L aboratory, Albuquerque, N M, 1985. IS . Richtmyer, R. D . a nd M orton, K.W., Difference M ethods/or Initial Value Problems, 2nd ed., Interscience Publishers, New Y ork, 1967. 16. Alifanov, O. M., Inverse Boundary Value P roblems o f H eat C onduction, J . Eng. P hys.l9(I), 8 21-830 (1975). 17. Beck, J. V., Review o f Six Inverse H eat C onduction C omputer Codes, A NL/RAS/LWR 81-1, Argonne National Laboratory, Argonne, Illinios, February 1981. 18. Beck, J. V., U ser's M anual for C ONTA-Program for Calculating Surface H eat F luxes F rom T ransient Temperatures Inside Solids, S andia N ational Laboratory, C ontractor Report, SAND83-7134, December, 1983. 19. Muzzy, R. J., Avila, J. H., a nd R oot, R. E., Topical Report: Determination o f T ransient H eat T ransfer Coefficients a nd t he Resultant Surface H eat F lux from Internal Temperature Measurements, General Electric, GEAP-20731, 1975. 20 . Bell, J . B. a nd W ardlaw, A. B., Numerical S olution o f a n I ll-Posed Problem Arising in Wind Tunnel H eat T ransfer D ata R eduction, N SWC T R 82-32, N aval Surface W eapons Center, Dahlgren, Virginia 22448, D ecember 1981. 21. Snider, D . M., Invert 1 .00A P rogram for Solving t he N onlinear Inverse H eat C onduction P roblem for One-Dimensional Solids, E GG-2068, E G&G I daho, Inc., I daho Falls, I daho 83415, February 1981. Bass, B. R., Drake, J. B., O tt, L . J., O RMDlN: A F inite Element P rogram for Two-Dimensional Nonlinear Inverse H eat C onduction Analysis, N UREG/CR-1709, O RNL/NUREG/CSD/ T M-17, U .S. N uclear Regulatory Commission, Washington, D .C., December 1980. 23 . O tt, L. J . a nd H edrick, R. A., O RINC: A O ne D imensional Implicit Approach t o t he Inverse H eat C onduction P roblem, N TIS R eport O RNLjNUREG-23, O ak Ridge National L abora- 22. REFERENCES ~k, J. V. a nd Wolf, H., T he N onlinear Inverse H eat C onduction P roblem, ASME Paper S-HT-40, p resented a t t he A SME/A IChE H eat T ransfer Conference a nd Exhibit, L os Angel'" CA, August 8 -11, 1965. ...., 2. ~OC 7600/CYB~R ~O M odel 76, C omputer Systems H ardward Reference Manual, Control a ta C orp., P ubbcauon no. 60367200, Minneapolis, 1972. 3. Blackwell, B. F., Efficient Techni~ue for the Numerical Solution o f t he One-Dimensional Inverse P roblem o f H eat C onductIon, Numer. Heat Trans/er 4 ,229 - 238 (1981). 4. Beck, J. V. a nd Arnold, K. J., Parameter Estimation in Engineering and Science Wil~y New ' , York, 1977. S. Beck, J. V., Surface H eat F lux Determination Using a n I ntegral M ethod N ucl Eng D 7 170- 178(1968). ' . . es. , I. 6. Beck, J. V., L itkouhi, B., a nd St. Clair, C. R., Effective Sequential S olution o f t he Nonlinear ......... ~ ... Inverse H eat Cond~ction P roblem, Numer. H eat Trans/er 5 ,275 - 286 (1982). 7. d~Boor, C ., A Pracllca/ GUIde to Splines, Springer-Verlag, New Y ork, 1978. 8. ~ Souza, N., Inverse H eat ~onduction P roblem for P rediction o f Surface T emperatures a nd e at Transfer from I ntenor T emperature Measurements, Report No. SRC-R-74 S ace Research C orporation, Montreal, D ecember 1973 'P D 'So.~ N:- Numerical S olution o f O ne-Dimensional Inverse Transient H eat C onduction . 9. b y FI~Jte DIfference Method, A SME P aper N o. 75-WA/HT-81, presented a t W inter Annual M cetlOg, H ouston, TX, Nov. 3O- Occ. 4, 1975. 10. Weber, C. F~ Analysis a nd S olution o f t he Ill-Posed Inverse H eat C onduction P roblem I nt. J . H eal Mass 1/'ans/er, 24(11), 1783- 1792 (1981). ' tories, O ak Ridge, T N ., 1977. Blackwell, B. F., User's M anual for the S andia O ne-Dimensional Direct a nd Inverse Thermal ( SODDlT) C ode, S andia N ational Laboratories internal report o f Div. 7537, 1983. 25. Drake, J . B~ ARIES : A C omputer P rogram for t he S olution o f F irst Kind Integral Equations with Noisy D ata, K /CSDfTM-43, C omputer Sciences a t O ak Ridge Gaseous Diffusion P lant, Post Office Box P , O ak Ridge, Tennessee 37830, O ctober 1983. 24. P ROBLEMS 6 .1. P rove that E q. (6.2.5) is valid for constant properties. 6 .2. Write the difference equations for the solution of Eq. (6.2.7) for the step sensitivity coefficient Z for constant properties and one-dimensional planar geometry. Use t he time integration of your choice. Demonstrate that the same difference equations can be obtained by differentiating Eq. (3.3 .45) with respect to q,.,. 6 .3. I n the algorithm for a single temperature sensor that exactly matches 286 CHAP. 8 O NE-DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM the temperature data, can the X sensitivity coefficient be replaced by Z in Eq. (6.3.8)? Why? 6 .4. Write an efficient computer subroutine. t o solve systems o f linear algebraic equations that are tridiagonal in form; see Eqs. (6.3.10) and (6.3.11). 6 .5. Verify the algebra for &is. (6.3.12) - (6.3.14), the algorithm for calculating the sensitivity coefficients. 6.6. Calculate the sensitivity coefficients Zl.1 a nd Z l.2 for x 1/L=0.5 and 1 for L1t+ = 0.05 a nd 0.1 respectively and for planar geometry with a perfectly insulated surface. Use the pure implicit method. Compare your results with the tabular values given in Table 1.1. CHAPTER 7 M ULTIPLE H EAT FLUX E STIMATION N 6.7. Prove that I W1 = 1 in Eq. (6.4.11). 1 =1 6 .S. Repeat Example 6.1 for x 1/L= 1 a nd for M + =0.05 a nd 0.5. 6 .9 . I n t he analysis of Section 6.6, the Z (step) sensitivity coefficient was used in the Taylor series expansion Eq. (6.6.4). Why Z and not X ? 7.1 I NTRODUCTION 6 .10. Repeat Example 6.2 for L1t+ =0.1. 6 .11. Extend the analysis o f Section 6.9 t o include multiple temperature sensors. 6 .12. C ompare the computational molecules for all the space marching methods o f Section 6.10. 6 .13. As a h eat transfer consultant, you have been asked t o comment on the d ata quality from the following inverse problem : A thermocouple is located 0.001 m below the surface o f a 0 .005 m thick copper slab. The known back face b oundary condition is a specified heat flux of 50 kW/ m 2 • T he person describing the experiment to you " thinks" the unknown heat flux you are trying to estimate is o f the order o f magnitude of 1 kW/m2. Will this experiment yield meaningful information? In the previous chapters the case o f the single unknown surface heat flux history was considered. I n this chapter the case o fthe multiple heat flux I HCP is treated. A multiple heat flux case arises when both surfaces of a one-dimensional body are exposed t o unknown heat flux histories. See Figure 7.1a for a plate heated o n b oth sides. The unknown heat flux histories are q1(t) a nd Q2(t). Figure 7.1b depicts either a hollow cylinder o r sphere. Again there are two unknown heat flux histories. In both Figures 7.1a a nd 7.1b there must be at 6 .14. a. F or t he D'Souza procedure, derive a two-node expression for the surface heat flux for a solid cylinder with the sensor a t the center. Derive the difference equations from an energy balance. b. Repeat p art (a) for three nodes. c. C ompare the expressions with that obtained using the exact ex pression given by Burggraf. (a) 6 .15. Solve Problems 6.14 for a solid sphere. 6 .16. F or the Weber procedure, derive a two-node expression for the surface f iGURE 7 .1 Some one-dimensional cases o f bodies exposed t o two unknown heat lIux histories, q I (t) and qa(/). (a), plate; (b), hollow cylinder o r sphere. heat flux for a solid cylinder with the sensor a t the center. The wave speed is infinite. Derive the difference equations from an energy balance. Compare the equation with the exact expression given by Burggraf. 6 .17. Solve Problem 6.16 for a solid sphere. (b) 287