Back to Inverse Heat Conduction: Ill-Posed Problems Home
Table of Contents
142 C HAP.4 I NVERSE HEAT C ONDUCTION E STIMATION P ROCEDURES S EC.4.5 regularization; the sum of squares function is J s= L j= 1 ( l}U-i}UI 4M =O-qU.:\t/>jO)2 J L + j =l ( l}.U+l-i}.u+d 4M =4M+I=O-QU.:\t/>jl-qU+l.:\t/>jO)2 + aWo(Qi,+Qi,+d+aW1(qU+I-QU)2 (4.5.29) T he estimator for qu is given by d lC22-d2CI2 Q u=C • • C 2 2- C 2 12 A (4.5.30) • J LL (4.5.3la) .:\t/>ll l = 0 J= 1 J C 22 =a(Wo + W1 )+ L j= .:\t/>lo (4.5.3lb) L .:\t/>jO.:\t/>jl j-. (4.5.3Ic) 1 J C l l = - aWl + Equation (4.5.33) is in exactly the same form as the r =2 sequential function specification method, Eq. (4.4.25c), and can be used in the same manner. The only difference is in the values o f the gain coefficients. See Example 4.2 for an illustration o f the sequential manner in which Eq. (4.5.33) can b e employed. Small o r large values of a are defined in relation to the .:\t/>~ a nd .:\t/>: values. Figures 4.9a a nd 4.9b display the gain coefficients for Eq. (4.5.33) for the case of a plate heated a t x = 0 a nd insulated a t x = L . T he sensor is a t x = L and the dimensionless time step is . :\t+ =0.05. T he t/>i values are taken from Table 1.1; the .:\t/>~ a nd .:\t/>: values are .:\t/>~ =t/>: = (0.00026W = 7.2 x 10 - 8 .:\t/>: =(t/>2 - t/>d 2 =(0.007885 -0.000269)2 = 5.8 x 10 - s where C l l = a(Wo+ W1 )+ 143 REGULARIZATION M ETHOD which suggests that 1 0- 7 to 1 0- 4 is a n important range for a values. An inspection of Figures 4.9a a nd 4.9b reveals that this is a range where interesting changes occur in K 1 a nd K 2' Figures 4.9a a nd 4.9b depict the gain coefficients for three different algorithms; Figure 4.9a is for the zeroth-order case and Fig. 4.9b for the first-order case. The function specification results are shown in both figures. The coefficients for the function specification method are independent of a; the K . value is 4.33 and K2 is 127. These values are also given in Table 4.1 for .:\t+ =0.05. J d. = L L (l}.UH -1J.uHI4M=4M l=O j= 1 d2= L ( l},U+l-f),u+d 4M =4M+I=O).:\t/>jO j=1 +I (4.5.32a) =O).:\t/>jl l~r----'----~-----r-----r----'---~ J F or the case of a single temperature sensor, J = I, the estimator for by Eq. (4.5.30) can be written as (4.5.32b) qu given 1 001=----- (4.5.33a) [cx(WO+ Wd+.:\t/>~].:\t/>o .:\ K =[a(Wo+ Wd.:\t/>. +aWl.:\t/>O] 2 .:\ (4.5 .33b) (4.5.33c) 0.1 .:\ = [cx(Wo + W d+ .:\t/>~ +.:\t/>n[a(Wo + Wd+.:\t/>~] - (aW. - .:\t/>o.:\t/>.)2 (4.5.33d) where the K . a nd K 2 expressions are gain coefficients. I t can be shown that this algorithm reduces to the Stolz algorithm for a ..... O a nd to the Q =C function specification algorithm for a ..... 0 0. (See Prob. 4.16.) - 0.01 L -._ _- --J'--__- I ._ _ _ _......J.._ ___- L-____ 1 0- 12 1 0- 10 1 0-8 1 0-6 1 0-4 0.01 ...J....l~---I a FIGURE 4 .9a Gain coefficients for zeroth-order sequential regularization algorithm given by Eq . (4.5 .33). ~/·=O.05 for sensor at x =L in a plate insulated at x =L. r =2, W o=I, WI=O. 144 C HAP.4 INVERSE HEAT C ONDUCTION ESTIMATION PROCEDURES 1~r----r----~---'----'-----~--~ Based o n this limited investigation, the first-order regularization a ppears t o be superior t o t he zeroth-order procedure b ut m ore analysis is needed t o verify this conclusion. 4.6 T RIAL F UNCTION M ETHOD 4.6.1 K ,. (Wo=O. W, = 1) 0.1 0.01 '-:-:--...J,.~_....L.-:-_""""'~---I_ _- -'-_ _. .J 10-'2 10-'0 10-0 10-6 1 0'" 0.Ql a FIGURE 4.9b Gain coefficients for first-order regularization algorithm given by Eq. (4.5 .33). r =2. Wo=O . W ,=I . T he z eroth-order method has Wo = 1 a nd Wt = 0 a nd t he f irst-order method h as Wo = 0 a nd W I = 1. F or t he s equential regularization method, the K I value for ex-+O s tarts a t 3720, which is the same as the S tolz m ethod value ( r= 1, also see T able 4.1). T his is t he s tarting value for the zeroth-order (Figure 4.9a), and the first-order (Figure 4.9b) m ethods, b ut K I decreases with increasing ex and approaches zero for t he z eroth-order case. N ear ex = 10- 6 , t he z eroth-order KI c urve h as a n inflection b ut t he first-order curve approaches 4.33, t he s ame value as given by t he function specification method. T he K 2 coefficients for both the zeroth- and first-order cases s tart n ear z ero for small ex's (less t han 1 0- 12 ) and then increase to a bout 127 with t he z eroth-order K2 d ecreasing for ex> 1 0- 6 , whereas t he first-order K2 r emains c onstant a t 127 for ex> IO- 9 .1t is r emarkable t hat t he first-order sequential regularization procedure (Figure 4.9b) for r =2 a nd At + = 0.05 yields almost identical results as the q = c onstant function specification procedure for the large range o f ex> 10 - 7. F or t he z eroth-order regularization method similar results a re o btained for ex between 10 - 7 a nd 1 0- s. F or large ex values (ex g reater t han 1 0- 4 in this case) b oth K I a nd K2 decrease toward zero for the zeroth-order regularization method. This c orroborates the previous assertion t hat t he z eroth-order method reduces the qM values for large ex's. Since, for ex g reater t han 1 0- 6, t he first-order K I a nd K 2 values approach t he s ame values as in t he function specification method, the first-order m ethod is seen to be less sensitive t o t he choice o f ex t han t he zeroth-order method. - 145 TRIAL FUNCTION M ETHOD SEC. 4 .6 I ntroduction A p rocedure t hat c ombines the function specification a nd r egularization approaches is called the trial solution method by Twomey. 33 .34 We shall call it the trial function method. T he m ethod as presented here is n ot identical t o that given by Twomey but the basic concept is the same. T he basic function t o be minimized is a s um o f s quares criterion plus a n a dditional term t hat is a generalization o f t he regularization term in Eq. (4.5.14); t he q vector is replaced by q m inus a predetermined " trial" function. T he t rial function c an i ncorporate prior information regarding t he s hape o f t he expected heat flux o r it c an b e a simple function, such as a constant. The trial function m ethod is potentially very powerful b ut a c hoice o f ex and o f t he f unctional form o f q m ust be made. F inding satisfactory combinations may r equire numerical experimentation. The derivation is given in matrix n otation a nd includes many interpretations. Sequential a nd w hole d omain (Twomey's interest) estimations a re i ncluded as is the possibility o f mUltiple t emperature sensors. A generalized least squares formulation is given t o p ermit n onconstant v ariance a nd/or c orrelated e rrors in the measured temperatures. The plan o f t his section is t o p resent a general matrix a pproach a nd t hen to discuss special cases, some o f which a re given in algebraic form. 4.6.2 M atrix Analysis The function t o minimize is s = (Y - Tf'" - I (Y - T) +ex{[Ho(q - q)]l'Wo[Ho(q - q)] + [ HJ(q _ q)]TW I [HI(q - q)]} (4.6.1) which is similar t o t he S function given by Eq. (4.5 . 14) for the regularization method. F or simplicity in presentation, the term analogous to the secondorder regularization term is n ot i ncluded in Eq. (4.6.1). T he s ymmetric square matrix '" is the covariance matrix o f t he r andom m easurement e rrors in Y. (See Reference 6, C hapter 6.) I t is a k nown matrix; if the first four s tandard statistical assumptions are valid, ' " a nd '" - I a re given by '" = where CJ2 is the variance o f Y; . C J2), '" - I = CJ - 2) (4.6.2) 146 C HAP.4 INVERSE H EAT C ONDUCTION ESTIMATION PROCEDURES E quation (4.6.1) c an b e used for b oth t he w hole d omain a nd sequential m ethods. F or t he case o f a single sensor, the whole d omain p rocedure has Y, T , q, a nd q v ectors with n elements, Y~[hl T~[n q-[!:l.tJ (4.~3J w here q is a trial function vector. F or t he sequential p rocedure t he subscripts in E q. (4.6.3) g o from M t o M + r-l a s s hown in Eq. (3.2.20). T he s quare m atrix is e qual t o I as given by Eq. (4.5.16a) a nd H I simulates first time differences a nd is displayed in E q . (4.5.16b). T he h eat flux vector q in Eq. (4.6.1) is u nknown a nd is t o be estimated. The symbol q d enotes t he t rial function vector which is a linear function o f q a nd a nother f unction q/, "0 (4.6.4) T he m atrix B is s quare a nd is selected in accordance with w hat t he analyst seeks t o a ccomplish; B is discussed later. T he v ector ql c an b e the p rior e stimate o f q before the algorithm is used t o e stimate q. I f, for example, the heat flux is k nown t o b e nearly c onstant n ear a c ertain value, t hen ql w ould be chosen to be this value. A nother e xample is when q is k nown t o v ary in a certain way with time such as a decaying exponential; ql c ould then be this decaying exponential function. I n t he r egularization p rocedure b oth B and ql a re z ero matrices. T he g eneral trial function e stimator for q with q r eplaced by Eq. (4.6.4) requires t he m atrix derivative o f E q . (4.6.1) with respect t o q. I n a ddition, an expression for T is needed; for t he l inear I HCP, T c an b e written as (4.6.5) As pointed o ut in C hapter 3, this expression is found for b oth t he D uhamel's n umerical a nd finite difference a pproaches. E quation (4.6.5) c an b e used for b oth t he w hole d omain a nd s equential procedures. F or t he whole d omain m ethod (4.6.6) where To is the initial temperature. F or t he sequential method, see Eq. (3.2.22). T he m atrix derivative o f Eq. (4.6.1) modified by the use o f Eqs. (4.6.4) a nd (4.6.5) yields, after collecting terms a nd s etting equal t o zero, ( X T" ,-IX + a{[Ho(l-B)]TWo[Ho(l- B)] + [ HI(I-B)yW I [HI(I- B)]})q = X T", - I(y - Tlq=o)+a{[Ho(I-B)YWoHoql + [HI(I-BWWIHlq/} (4.6.7) This matrix e quation r epresents a set o f n s imultaneous, linear algebraic e quations for q l" . . , q. for the whole d omain m ethod. F or t he s equential method, SEC. 4.6 1 47 T RIAl. FUNCTION M ETHOD the unknowns a re qM,' . . , qM +, - 1 b ut only qM is used. E quation (4.6.7) is valid for b oth single a nd m ultiple sensors. In o rder t o e xpedite t he e xamination o f E q . (4 .6.7), several simplifications are used; ' " is given by (121, Ho by I, Wo by I, a nd t he first-order term is d ropped by setting W I = 0. T hen symbolically, II c an be given by q =[(1-2XTX + a(I-B) T(I-B)r I [(1-1XT(y - Tlq=o)+a(I _ B)Tq/] (4.6.8) Many special cases c an be considered in connection with this e quation, several of which a re briefly discussed in the following. 4 .6.3 Z eroth-Order R egularization M ethod The z eroth-order r egularization m ethod c an b e o btained from Eq. (4.6.8) by setting ql = 0 a nd B =O (4.6.9) F or t he whole d omain m ethod, E q . (4.6.8) is used with E q . (4.6.6) . See also Eq. (4.5.20). E quation (4.6.8) c an a lso be used for the sequential. r~gularization m ethod; if there a re J t emperature s ensors then Y c an be partitIOned so t hat YM) [ YI.M+i) Y = ! M+I , Y M+ i = ~l.M + i [ J Y M+,-I Y •M+i (4.6.10) where y . M + . is u nderstood t o b e a n e lement o f t he Y v ector a nd n ot p art o f a two-di~~ns;onal a rray; t hat is, Y h as dimensions o f r J x I . An example o f t he resulting e quations is given for r =2 by Eqs. (4.5.30) a nd (4.5.31) where a is replaced by au 2 i n Eq. (4.6.8) a nd Wo = 1 a nd WI = 0 . 4.6.4 G eneralized S equential F unction S pecification M ethod T he trial function m ethod c an be used t o yield a generalization o f t he sequential function specification m ethod . F or t he q =constant t emporary a pproximation, ql = 0 (4.6.11) 1 0 0 0 .. .) 1 0 0 0 ... B = 1 0 0 0 . .. (· . . . (4.6.12) ·... ·... which results in the trial function q being the c onstant vector o f [ see Eq. (4.6.4)] q =qMI (4.6.13) 148 C HAP.4 INVERSE HEAT CONDUCTION ESTIMATION PROCEDURES F or the case of r =2 and a single sensor, Eq. (4.6.8) is equivalent to solving the equation, a - 1[(A4>o)1+(A4>I)1]+cx a - 1A4>oA4>I-CX][qM ] [ a - zA4>oA4> I - cx a - 1(A4>o)1 + cx qM + I (4.6.14) 1 = [a-:(YM - TM I)A4>o+a- (yM + 1 - t M + II)A4>I] a - (YM+ I -tM+ I I)A4>o This equation is solved only for qM because qM + I is not used in the sequential procedure. Solving this equation for qM and putting the equation in the gain coefficient form yields Eq. (4.5.33) if Wo=O, WI = 1, and cx-+cxa 2. I t is surprising that for r =2 this trial function method gives the same result as the first-order regularization method. Notice that the Wo=O, WI = 1 curves in Figure 4.9b (which is equivalent to the trial function r =2 case) approach the q=constant function specification result for (Xa 2 greater than 10- 1 , which is a large range of values indeed! The trial function method with B given by Eq. (4.6.12) is shown by this example to be a generalization of the q = constant function specification method; this conclusion is also valid for r greater than two. 4 .7 4.7.1 FILTER F ORM O F L INEAR I HCP I ntroduction The inverse heat conduction algorithms given in previous sections of this chapter can be re-formulated as digital filters. The only restriction is that the problem be linear; the sequential and whole domain procedures can be considered and also single and multiple sensors can be included. Moreover, the method o f approximating the heat conduction problem (using Duhamel's theorem, finite differences o r other) does not affect the validity of the digital filter approach. The digital filter approach is important because it can be demonstrated that it is much more computationally efficient than other methods. Due to its efficiency, it can be readily implemented in an on-line method of analysis. Heat flux measuring devices can incorporate the digital filter concept and immediate visual digital output can be provided. (By on-line is meant a measurement that is given in real time, but the I HCP algorithm may necessitate a delay of a few time steps.) 4 .7.2 S equential F ilter A lgorithm All the sequential algorithms derived in this chapter for the linear I HCP can be given in terms of gain coefficients. See, for example, Eq. (4.4.25) which can be SEC.4.7 149 FILTER F ORM OF LINEAR IHCP written for a single sensor as qM= i: Ki(YM +i-I-tM+i-tlq.,=" , =o) (4.7.1) i= l where the gain coefficients, K;. depend on the specific algorithm such as those for the function specification, regularization, o r trial function procedures. F or present purposes the Ki are assumed known. . . . A filter form of Eq. (4.7.1) is based on each estimated qM bemg linear m the Y; values. This linearity was demonstrated several times in previous sections a~d in connection with Eqs. (4.3.6) - (4.3.9). Linearity implies that the effects of various components of 1'; can be individually determined and the total effect can be obtained by superposition. . To illustrate the linearity, the case of a finite plate heated at x = 0 and msulated at x = L is considered. The sensor is a t x = L and the dimensionless time step of A t + = 0.05 is used with r = 3 future times in the function specificati~n method with q temporarily held constant. Equation (4.4.25) is used and the gam coefficients from Table 4.1 are K I =O.292, K z =8 .56, and K )=31.8 (more significant figures were used in the computations). T~e initial temperature, To , is set equal to zero and a series o f separate problems IS solved. The first of these problems has Y = 1 and the subsequent 1'; values are equal to zero. The secof!d I problem has 1'; = 0 except Y = 1 and the j th problem has 1) = 1 an~ 1'; ~ 0 for 2 i h. The associated heat fluxes for each of these problems are gIven m the respective columns of Table 4.4. There are several important characteristics of Table 4.4. Except for the first two columns the calculated heat flux components are the same in each column but each s~cceeding column is shifted down one line. In order to describe T ABLE 4 .4 H eat F lux C omponents f or Y j =1 a nd Yj=O f or i n: V alues a re G iven i n C olumns f or j =1 t o 6 . ( For , =3. Al+=O.05) All Y; V alues E qual Z ero E xcept: Y =1 Y =1 Y =1 4 3 2 I 2 3 4 5 6 7 0 .29 - 0.35 - 0.02 0.06 0.03 O.OOOS S 9 o S.6 - 10.1 - 0.9 1.7 0.9 0.05 - 0.15 - O.OS - 0.003 0 Y s=l Y =1 6 31.S - 29.9 - 12.1 5.4 4.7 0.97 - 0.51 - 0.41 0.006 0 31.S - 29.9 - 12.1 5.4 4.7 0.97 -0.51 - 0.41 0 0 31.S - 29.9 - 12.1 5.4 4.7 0.97 - 0.51 0 0 0 31.S - 29.9 - 12.1 5.4 4.7 0 .97 0 0 0 0 1 60 C HAP.4 INVERSE HEAT C ONDUCTION ESTIMATION PROCEDURES SEC. 4.1 151 FILTER FORM OF LINEAR I HCP the results of Table 4.4, a notation is needed. Let TABLE4.S P attern f or D igital F ilter C oefficients f or T able 4 .4 t5q", t5lj M represent the entries of the j th column of Table 4.4. T he notation denotes the change in q", when l j is increased one unit. With this notation, observe that for this case with r = 3, t5q", = 0 t5ql t5q2 t5q", t5Y) = t5Y4 = . .. = t5Y"'+r-1 t5q2 - t5Y) t5q) t5q", t5Y4 t5Y"'+r-2 = - = . .. = Thus for j~r, the same entries in a given column of Table 4.4 a ppear in the other columns. Consequently if one column is calculated (for j ~ r), the entries for the subsequent columns can be inferred from the one calculated. The first two columns o f Table 4.4 are noted t o be different from the subsequent ones. This is true in general for the columns 1,2, . .. , r -1. I n order to treat each d ata p oint in the same manner, the calculations for sequential IHCP algorithms should start a t least r - I time steps before heating starts. If this is done, the anomalous first few steps have negligible effect on the calculations. The entries in Table 4.4 for columns 3, 4, . .. can be given as the filter coefficients f -2' f -I> f o, fl> f 2,'" as shown in Table 4.5. These f values are related to hq",/t5lj by t5q", Y,,=l Y s=l Y =1 6 0 0 0 0 0 0 5 6 7 31.8 - 29.9 Y =1 2 1 2 3 4 for M < j-r+ 1 t5lj Y =1 3 1- 2 I- I 10 11 12 13 I" Y I=l 1-2 1- 1 10 11 12 13 t5Y,. t5q", t5q2 t5Y"'+r-2 Hr f - I = t5q", = (jqr - I hY",+1 h Yr A (jq", (jq", (jq", q "'=W(YI-To )+ ( jy-(Y2 -To)+ . .. + (jy' A q ",= (jq"' +r- I (jy' r ( YI-To)+ (jq"'+r-2 (jy' r (jql ( Y2 -To)+"'+(jy'(Y",+r-I- To) r + f l(Y",-1 - + fo(Y", - To)+ f - I(Y",+ 1 - To)+ . .. (4.7.3b) To) + f 2-r(Y"'+r-l- To) (4.7.4) + f l - r(Y"'+r - , - To) Notice that the subscripts of f plus those of Y in Eq. (4.7.4) always equal M . Provided the calculations for q", start a t least r - 1 steps before heating starts, the I HCP algorithm given in the gain coefficient form by Eq. (4.7.1) can be equivalently written in the digital filter forms as (4.7.5) = (jq", = (jqr+ I (jY",-1 bY,. M +r-l q ",= ~y. (jq", (4.7.2) "'+r- i- I (jy" Note that f l-r is not dependent on M but on the difference in the subscripts of (jq", a nd (jY"'+r-I ' Notice also that the first equal sign in each equation of f l - r+i ( Y"'+r_I-To)(4.7.3a) I 2 " '+r-I which looks like a Taylor series expansion with the higher-order terms neglected. However, Eq. (4.7.3a) is an exact result for linear problems. Substituting Eq. (4 .7.2) into Eq. (4.7.3a) gives f o= t5q", = (jqr (jY", (jy" fl 1-2 1- 1 10 11 Eq. (4.7.2) is for a row in Table 4.5, t hat is, fixed M, a nd the second equal sign is for a column of Table 4.5, i.e., fixed r. Using these symbols and the principle of superposition, the linear relation for q", can also be written as = f", - I(YI - To)+ f ",-2(Y2 - To)+ . .. t5Y"'+r-1 1- 2 1- 1 10 11 12 U L i= 1 f", - i(ll-To ) (4.7.6) Note that these equations indicate that q", depends on the M - 1 previous temperature measurements plus the r future temperature measurements. However, in practice the lower limit in the summation need not start at 1 because the associated filter coefficients become negligible as demonstrated in Table 4.4. 1 62 CHAP. 4 INVERSE H EAT C ONDUCTION ESTIMATION PROCEDURES T he digital fi~ter coefficients given by f o r {)q/{) Y c an be calculated by using the I HCP algorIthm such as Eq. (4.7.1) with SEC. 4.8 T WO CONFLICTING OBJECTIVES 1 53 These are the same values as given by the gain coefficient form of the original I HCP algorithm, Eq. (4.7.1). 0 Y ,-To=l, 1 ';- To=O for i fr (4.7.7) T he resulting values of liM a re (4.7.8) This equation can be verified by using the Y = 1 columns o f Tables 4.4 and 4.5. 3 The fore~oing discussion o~ a digital filter is in connection with the sequential method. I t IS d emonstrated I n C hapter 5 that the whole domain estimation ~Igorithms c an be written in exactly the same form as Eq. (4.7.5), provided r IS chosen to be a sufficiently large integer. EXAMPLE 4.5. F or the case o f a plate insulated a t x ~ L a nd heated at x = 0, the measu ~ed temperatures at x~L for times 0.05, 0.1, a nd 0.15 s a re 30.269, 37.885, and 59 .306 C. Rel~vant v~lues a re L =I m, cx=1 m 2/ s, k =1 W/ m-C, and T o=30°C. The value o f r =3 IS used I n t he function specification method. Use t he filter form o f the algorithm to calculate the first three values o f heat flux . 'L.., o"~e"1) Solution. T he equation used is Eq. (4.7.5) b ut the calculational starting time must be moved a t least 2 times early. Let the d ata be arranged as : >; Time(s) -O/!J S- --@ 0 0.05 0.1 0.15 2 3 4 5 30 30 30.269 37.885 59.306 where the first calculation is a t t = - 0.05 s r ather than at t=0.05 s. F or t he first time o f M = I , Eq. (4.7.5) gives (using column 3 o f T able 4.4) ~ oqJ Oq2 Mi, q, = . y (Y, - To) + <Y ( Y2 - To)+ -. - (YJ - To) UJ UJ oYJ ~ - 12.1(30-30)-29.9(30-30)+31.8(30.269-30)=8 .55 W/m2 which corresponds to time t = - 0.025 s. T he value for q2 is o btained using (4.7 .5) t o get t>( I ~ o/i. .5/iJ OQ2 .5Q, q 2= . y ( Y,-To)+.-y ( Y2 -To)+.- (YJ-To)+-(Y.- To) UJ UJ uYJ oYJ = 5 .4(30- 30)-12.1(30- 30)- 29.9(30.269-30)+ 31.8(37.885 - 30)= 243 W/m2 / iJ=693 W/m2 4.7.3 P refiltering T emperature M easurements The subject of digital filters 35 is i mportant in the field of signal processing. Some o f the digital filter algorithms can be applied to the temperature measurements before the I HCP algorithms are employed. There are many such filters possible and their discussion is beyond the scope of this book. I f the statistics of the temperature measurement errors are known, then it is possible to design filters to admit only "low" frequency signals. F or more discussion of digital filter design, see Reference 37. One possible prefilter (Reference 37, p. 57) is 1';-1 + 21';+ 1';+1 4 (4.7.9) That is, each 1'; value is replaced by Yi before the I HCP algorithm is applied. More specifically, Yi o f Eq. (4.7.9) replaces 1'; in Eq. (4.7.5) o r Eq. (4.7.6). Notice that the coefficients of 1';-1, 1';, a nd 1';+ 1 sum to unity. See Problem 4.15. 4.8 4.8.1 T WO C ONFLICTING O BJECTIVES M inimum D eterministic B ias One o f the characteristics of good estimators is that of minimum bias; in fact, in parameter estimation problems, an unbiased estimator is usually sought. Mathematically, it is desirable for the estimator, p, o f a parameter, /1, t o be the expected vall!e o f /1, (4.8.1) In o ther words, the expected value o f the estimate is the true value. I f Eq. (4.8.1) is true, the estimator fJ is said to be unbiased and " on the average" the estimate given by would be near /1. I f Pwere a biased estimator of /1, would on the average tend to be either higher o r lower than the true value o f /1. In ill-posed problems it is preferable for the estimators for the heat flux components, l il' ... , li., t o have low bias. I t is not required that the bias be zero, however. In ill-posed problems, the bias and variability (variance) of the estimators for the unknown parameters are related. Consequently the bias is only one aspect of the estimators that must be considered. The bias of a n estimator can be investigated when the random measurement errors are set equal to zero. This bias o r e rror in the estimator can be called the deterministic bias. I t is advantageous to make it as small as possible and yet not make the variance unacceptably large. p p 164 C HAP. 4 INVERSE HEAT C ONDUCTION ESTIMATION PROCEDURES SEC. 4 .8 4 .8.2 M inimum S ensitivity t o R andom E rrors 1 55 T WO CONFLICTING OBJECTIVES side of Eq. (4 .8.4) gives Another characteristic of a good estimator in addition to minimum bias is that o f minimum variance. I f Pis the estimator of the parameter p, then this means that the variance of p, f/~ = E({[qM - E(qM)] - [iJM - E(qMW) (4.8.2) - 2E{[qM - E(qM)][iJM- E(qM)]} should be a minimum. F or linear parameter estimation problems with the first, second, seventh, and eighth standard assumptions satisfied, the estimator that is a linear function of the measurements and has unbiased and minimum variance is called the Gauss-Markov estimator (Reference 6, p. 232). F or ill-posed problems, better estimates are obtained if the requirement of unbiased estimators is not imposed. In the latter case, it is necessary to specify an alternative objective which is discussed next. 4 .8.3 = E{[qM - E(qM)]2} (4.8.5b) + E{[qM-E(qMW} The first term on the right side of Eq. (4.8.5b) is the variance of the estimator, q~h (4.8.6) The second term on the right side of Eq. (4.8.5b) can be shown to be zero because iJM - E(qM) is not a random variable and also E[qM - E(qM)] = E(qM)-E(qM)=O M ean S quared E rror F or ill-posed problems the common requirements of zero bias and minimum variance do not yield satisfactory estimators. F or zero bias, the estimators are very sensitive to measurement errors; that is, the variance o f Pis large. There are ways to reduce the variance for the I HCP such as requiring that all the components of q be equal, or 9 i,=[qM-E(qM)]2 (4.8.8) Hence the mean squared error of qM is composed of two parts, a variance and the square of a deterministic error, and the relationship between them is f /i,= V (qM)+9i, (4.8.9) An expression for V(qM) is developed in Section 4.8.4 and the deterministic error is discussed in Section 4.8.5. The significance of the components of the mean squared error given by Eq. (4.8.9) can be illustrated by using Figure 4.10. T he continuous solid line is a q",\+ (4.8.4) In general, estimators that minimize f / ~ for all values of M, M = I, 2, . .. , n, are sought. I t is i mportant to observe that Eq. (4.8.4) does not give the variance of qM, that is, f /! 4= V(qM), except for the special case when the expected value of qAt is the true value, iJM; this only occurs if qM is an unbiased estimator of iJAt. Since for ill-posed problems the estimator is usually biased, Eq. (4.8.4) is not, in general, the variance of qM. I t can be shown that Eq. (4.8.4) includes two parts: deterministic and stochastic. Adding and subtracting the expected value of qAt inside the right (4.8.7) The last term in Eq. (4.8.5b) is the square of a bias and the outer expected value symbol can be dropped; the resulting expression is given the symbol 9 ~, (4.8.3) I f a requirement such as Eq. (4.8 .3) is introduced, the variance of qi is relatively small but it usually has a large bias. Hence, the deterministic bias and variance of the estimator are related and an optimal strategy should consider both aspects. A function that considers both bias and variability is the mean squared error. Let qM be a component of the estimated q vector and let qM be the true value of the component. The mean squared error of qM is (4.8.5a) + q + • -+ • • •••• + t!JJ", •I + o •• + • + I + - I I q", E(q",l FIGURE 4 .10 The true heat flux q", a t time '",. the estimated value fj.,. and the mean o f the estimated value. E(fj.,). 168 C HAP.4 INVERSE HEAT C ONDUCTION ESTIMATION PROCEDURES r epresentation o f a " true" h eat flux t hat s tarts a t time zero a nd a t time to j umps t o a c onstant value. T he t rue heat flux a t time tM is shQwn by a s quare in Figure 4.10 a nd is denoted qM' T he d ots i llustrate estimated values o f t he h eat flux curve for errorless measurements o r equivalently t he e xpected (i.e., theoretical average) value o f t he estimated heat flux, E(qj). N ote t hat t he E(qj) values are biased since a t m ost times the values a re e ither greater o r less t han t he true values. F or example, a t time t M , t he m ean e stimated value is low. T he difference between qM a nd E(qM) is t he d eterministic bias, 9 M. T he crosses in Figure 4.10 show simulated a nd c alculated heat fluxes with measurement errors; t he value a t t ime tM is denoted qM' D ata from o ther sets o f m easurements would yield a set o f e stimated heat fluxes a t tM t hat w ould center a bout E(qM), r ather t han the true value qM' A m easure o f t he variability o f t he qM values with respect to E(qM) is called the s tandard d eviation o f qM a nd t he s quare o f t he s tandard d eviation is t he v ariance o f qM which is denoted V(qM) ' T he m ean s quared e rror given by Eq. (4.8.9) contains b oth t he effect o f m easurement e rrors which is described by t he variance, V(qM), a nd t he effect o f t he d eterministic bias, 9 M, which is d ue t o a b iased I HCP a lgorithm. In t he inverse h eat c onduction problem s ome b ias is a ccepted in o rder t o r educe t he e xtreme sensitivity t o m easurement errors, particularly a s t he t ime steps a re m ade small. 4 .8.4 V ariance o f E stimated H eat F lux C omponent I n Sections 4.3- 4.6 two types o f e stimators are derived: sequential a nd w hole domain. F or b oth types t he a lgorithms c an be written in the form o f a d igital filter. An example is Eq. (4.7.5) which is convenient for determining a n expression for the variance o f qM because t he l j c omponents a re explicitly given, in c ontrast t o t he implicit dependence in t he g ain coefficient equation, Eq. (4.7.1). T o derive a n e xpression for t he v ariance o f qM' s ome statistical assumptions a re necessary. T he first four s tandard s tatistical assumptions a re a ssumed valid; t hat is, t he r andom e rrors, £j, in l j a re additive, have zero mean, have constant variance, ( 12, a nd a re u ncorrelated, SEC. 4 .8 T WO CONFLICTING OBJECTIVES = (12 L M +,-I j=1 1 57 (b"i)2 3 bY,. (4.8.llb) N"otice t hat .the equality o f E q. (4.8. 11 a) a nd (4.8.1Ib) does n ot im I h bqM+,-i/bY,. IS e qual t o bqi/b Y,.. P Y t at EXAMPLE 4.6 Calculate the. variance of tJ", for a fiat plate insulated at x = L and heated at x = O. The sensor IS a t x = L. Use At = 0.05 s for L = I _ I 2; k W /m-C,andr=3. m ,ll- m s, = 1 Solution. T~e Mi/lJ Y, values are given in the Yl = 1 column o f Table 4 4 F M - 1 EQ. (4.8.11) gives . . or - , ql V(qd=0'2 [(lJ )2 + (lJtJ2)2 + (lJ ql )2] . lJYl lJYl lJYl =0'2[(31.8)2+( - 29.9)2+( -12.1)2J=20520'2 W2/m4 and for M = 2, one finds V(Q2) = 0'2[(31.8)2 + ( - 29.9)2 + ( -12.1)2 + (5.4)2J =20800'2 W2/m4 ~or Ql' q4,.q~, q6' and q, the variances are 21030'2, 21040'2, 21040'2, and 21040'2 res _ lively. Notice that the varian~ approach a constant value as M increases. ' pee The SQ.uare ro~t o f the vanance is the standard deviation which can be compared ~.the eSllmated q", values. For O' 2=2(C)2 the standard deviation of q is 65 W/m 2 IS value can be compared to the value of q l=693 W/m2 that was ~alculated f o' Example 4.5. r o 4.8.5 F lux E stimate o f D eterministic E rror i n S urface H eat ~here a re several ways ? f e stimating the deterministic e rror in the surface heat ux. T wo o f these a re discussed here. Suppose t hat t he true surface heat fluxes are the values o f q j=O . q ;=bq, for i fr (4.8.12a) (4.8.10b) (4.8.12b) a re calculated . E (3 2 2 . I at times t t E ( 32 1 ,.2,'" u smg q. . .1) with To e qual t o z ero. An expanded form o f q. . .12) IS (4.8.lOc) (4.8.13a) (4.8.10a) cov(£j, £ )=0 for i fj for i =r Th~ c orrespondmg temperature rises a t t he interior location x F or these assumptions, the variance o f qM given by Eq. (4.7.5) is 7; = L!¢lql +L!¢oq2 (4.8 .13b) T3 =:L!¢2ql +L!¢lq2 +L!¢oq3 (4.8.1Od) (4.8.13c) M (4.8.1la) TM = L L!¢M_jqj i =l (4.8.13d) 1 58 CHAP.4 INVERSE HEAT C ONDUCTION ESTIMATION PROCEDURES 169 REFERENCES T ABLE 4 .6 R esults f or E xample 4 .7 F or simplicity, consider the case o f r = 2 and use Eq. (4.8.12) in Eq. (4.8.13) to get t he simulated temperature measurements of (4.8.14a) M bq",/bqr S um f) (4.8.14b) 1 2 3 4 5 6 7 8 9 10 0.00856 0.23427 0.45060 0.29232 0.06487 - 0.02623 - 0.02382 - 0.00513 0 .00246 0.00209 0.00042 - 0.00020 - 0.00021 0.00856 0.24283 0 .69343 0.98575 1.05062 1.02439 1.00056 0.99544 0 .99790 0.99998 1.00040 1.00020 0 .99999 0.00856 0.23443 0.59733 0.66502 0.66818 0.66869 0.66912 0.66913 0.66914 0.66914 0.66914 0.66914 0.66914 (4.8.14c) (4.8.14d) These values of Y are then introduced into I HCP algorithms such as Eq. (4.4.25) t o calculate the heat flux estimates, [,qM' T he ratio of [,fi.M/[,qr is independent of [,q,; this notation o f M iM/['q, makes clear that the calculation yields changes in the surface heat flux for temperatures associated with a change in the heat flux a t time t r • Because [,qM is the change in qM given by an I HCP algorithm and b q, is the change in the input heat flux, b qM/bq, with M = r is not equal to unity (except for the Stolz algorithm). The maximum deviation from the input qi values given by Eq. (4.8.12) is o ne estimate of the deterministic error, 9 , q., [,q. 9 = max ,~ - -.!.. bqr i [,qr [,qr (4.8.15) C I) 9= [ i= I bqr bqr gain coefficients a re K,=0.292046. K 2=8 .56054. a nd · K)=31.8168 (see T able 4.1). T he bq,/bq) v alue is o btained f rom Eq. (4.4.25) a s bq, F or most cases, t he; value in Eq. (4.8.15) t hat yields a maximum for 9 corresponds to ; = r. Another estimate o f the deterministic error 9 is obtained from the square root of the sum of squares of the differences between MiM a nd the input values given by Eq. (4.8.12). Mathematically, this estimate of the deterministic error is q L ([,qi)2 + (b , -1 )2J 1/2 ['q, - 11 12 13 - bq3 ~ a nd Mi2/bq3 is equal t o bq2 - bq3 = K )Aq,o = 0.00856 5 b ~ =K,(O-O.oo@jAtP,) + K 2(AtPo - 0.00856 AtP2) + K )(AtP, (4.8.16) = 0.23427 if" F urther v alues a re d isplayed in t he s econd c olumn o f T able 4 .6. A r unning s um o f t he Surprisingly this definition of 9 frequently yields values only 20 o r 30% larger than that given by Eq. (4.8.15). Hence the deterministic error 9 is usually nearly the same value for both definitions, Eq. (4.8.15) o r Eq. (4.8.16). This is important because the two measures emphasize different aspects. The first one, Eq. (4.8.15), o btains its sole contribution at ; = r o r very near r, whereas the second, Eq.· (4.8.16), obtains contributions for an extended range of ;'s. Even so, the results tend t o be close. In later chapters, Eq. (4.8.16) is used rather than Eq. (4.8.15). lJii", /bq) values. given in t he t hird c olumn o f T able 4.6, a pproaches unity, a s it s hould for c onservation o f e nergy. T he f ourth c olumn gives t he d eterministic e rror a s defined by E q . (4.8.16); t he e ntries a re r unning s ums o f t he t erms in Eq. (4.8.16). N ote t hat t he v alues quickly converge t o t he value o f 0.66914. T his v alue c an b e c ompared t o t hat o btained f rom Eq. (4.8.15) which is f ) = 1 -0.45060=0.54940 T he v alue given by Eq. (4.8.16) is only 2 0% l arger t han t his volume. Since a n a pproximate m easure o f t he d eterministic e rror is satisfactory, b oth e quations c an b e used; 0 Eq. (4 .8.16) is p referred in t he r emainder o f t his book. F or a flat p late i nsulated at x = L. a nd h eated at x =O a nd t he s ensor at x = L . calculate t he d eterministic e rror for t he f unction specification a lgorithm with r=3. U se AI = 0.05 s for L = 1 m. cc= 1 m 2/s. a nd k = 1 W /m-C. ExAMPLE 4 .7 So/ulion. T he d eterministic e rror c an b e found from b oth Eq. (4.8.15) a nd E q. (4.8.16). I n b oth c ases t he lJijdbqr values a re n eeded which a re c alculated using t he function specification a lgorithm w ith t he i nput t emperatures o f Y, = Y = 0. Y )=AtPo=0.000269; 2 Y4 =AtP, = 0.007099. .... See T able 1.1. T he a lgorithm is Eq. (4.4.25) with r =3. T he 0.00856 AtP 3) REFERENCES I. H adamard. J .• Lee/ures on CauL'hy"s Problem in Linear Partial DijJeren/ial Equa/ions. Yale University Press. New Haven. cr. 1923. 160 C HAP.4 INVERSE H EAT C ONDUCTION ESTIMATION PROCEDURES 2. T ikhonov, A. N. a nd Arsenin, V. Y., Solutions o f Il/-Posed Problems, V. H . Winston & Sons, Washington, D.C., 1977. 3. Widder, D. V., T he Heat Equation, Academic Press, New York, 1975. 4. Weber, C. F ., Analysis a nd S olution o f t he III-Posed Inverse Heat Conduction Problem I nt.j. H eat Mass Transfer 2 4,1783 1792 (1981). ' 5. O lmsted, J. M. H., Advanced Calculus, Prentice-Hall, Englewood Cliffs, NJ, 1961. 6. Beck, J. V. a nd Arnold, K. J., Parameter Estimation in Engineering and Science, Wiley, NY, 1977. 7. Beck, J . V., C riteria for Comparison o f M ethods o f S olution o f t he Inverse Heat Conduction Problem, N I 'd. Eng. Des. 53, 11 - 22 (1979). 8. d eBoor, Carl, A Practical Guide to Splines, Springer-Verlag, New York, 1978. 9. Stolz, G ., Numerical Solutions t o a n Inverse Problem o f H eat Conduction for Simple Shapes, J . Heat Transfer 8 2,20 - 26 (1960). 10. F rank, I ~ An Application o f Least Square Methods t o t he Solution o f t he Inverse Problem o f H cat C onduction, J . H eat Transfer 8 5,378 - 379 (1963). 11. Davies, J. M., Input Power Determined F rom T emperatures in Simulated S kin P rotected Against Thermal Radiation, J . Heat Transfer 8 8,154 - 160 (1966). 12. Beck, J. V., C orrection o f T ransient Thermocouple Temperature Measurements in HeatConducting Solids, P art II, T he C alculation o f T ransient H eat Fluxes Using the Inverse Convolution, A VCO C orp ., Res. a nd Adv. Dev. Div., Wilmington, MA, Tech. Report RADTR-7-60-38 ( Part II), March 30, 1961. 13. Beck, J . V., C alculation o f Surface H eat Flux from a n I nternal Temperature History, ASME P aper 6 2-HT -46 (1962). 14. Beck, J. V., Surface H eat Flux Determination Using an Integral Method, N ud. Eng. Des. 7, 170- 178 (1968). 15. Beck, J. V., N onlinear Estimation Applied t o the Nonlinear H eat C onduction Problem, Int. J . Heat Mass Transfer 13, 703 - 716 (1970). 16. Blackwell,. B. F., An Efficient Technique for the Numerical Solution o f t he One-Dimensional Inverse Problem o f H eat Conduction, Numer. Heat Transfer 4, 229 - 239 (1981). 17 . Alifanov, O . M., Inverse Boundary Value P roblems o f H eat Conduction, J . Eng. Phys. 25 (1975). 18. Alifanov, O . M . a nd A rtyukhin F. A. "Regularized N umerical Solution o f N onlinear Inverse Heat-Conduction P roblem," J . Eng. Phy. 29, 934- 938,1975. 19. Alifanov, O . M ., Identification o f Processes o f Heat Container Apparatus. A n Introduction to the Theory o f Inverse Problem o f H eat Transfer, M achinery Publisher, Moscow 1979 (in Russian). 20. Cozdoba, L. A. a nd Crykowsky, P. G., Methods o f Solution to the Inverse Problem, Scientific Publisher, Kiev, 1982 (in Russian). 21. Bell, J. B. a nd Wardlaw, A. B., Numerical Solution o f a n III-Posed Problem Arising in Wind Tunnel Heat Transfer D ata Reduction, Naval Surface W eapons Center, N SWC T R 82-32, Dec. 1981. 22. Bell, J. B., T he Noncharacteristic Cauchy Problem for Class o f E quations with Time Dependence. l . Problems in O ne Space Dimensions, S I A M J . Math. Anal. 12, 759- 777 (1981). 23. Beck, J. V. a nd M urio, D., Combined Function Specification Regularization Procedure for Solution o f Inverse Heat Conduction Problem, AIAA P aper No. AIAA-84-0491, January 1984. 24. Hoerl, A. E. a nd K ennard, R. W., Ridge Regression Biased Estimation for Nonorthogonal Problems, Technometrics 12, 55 - 67 (1970). 25. Draper, N. R. a nd Van N ostrand, R. c., Ridge Regression - Is I t Worthwhile, Univ. o f Wisconsin Department o f Statistics Technical Report 501, M arch 1978. 26 . . Levenberg. K., A Method for the Solution o f C ertain Non-linear Problems in Least Squares, Q. Appl. Math. 2, 164- 168 (1944). 27. M arquardt, D. W., An Algorithm for Least Squares Estimation for Nonlinear Parameters, J . Soc. Ind. Appl. Math. 11,431 - 441 (1963). 28. M arquardt, D. W., Generalized Inverses, Ridge Regression, Biased Linear Estimation, a nd N onlinear Estimation, Technometrics 12, 591 - 612 (1970). 161 PROBLEMS 29. 30. 31. 32. 33. 34. 35 . 36. 37. Lawson, C. L. a nd H anson, R. J ., Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs, NJ, 1974. T he I MSL Library, International Mathematical a nd Statistical Libraries, Inc., Houston, TX. G raham, N. Y., S moothing with Periodic Cubic Splines, Bell System Tech. J. 62, 101 110 (1983). Reinsch, C. H. J., Smoothing by Spline F unction, Numerische Muthematik 10, 177 183 (1967). Twomey, S., O n the Numerical Solution o f F redholm Integral Equations o f the First Kind by the Inversion o f the Linear System P roduced by Q uadrature, J . Assoc. Camp. Mach. 10, 97 101 (1963). Twomey, S., T he Application o f N umerical Filtering t o t he Solution o f Integral Equations Encountered in Indirect Sensing Measurements, J . Franklin Ins/. 279, 95 - 109 (1965). Blackwell, B. F ., Some C omments o n Beck's Solution o f t he Inverse P roblem o f H eat C onduction T hrough t he Use o f D uhamel's Theorem, Int. J . Heat Mass Transfer 2 6,302 305 (1983). Alliney, S. a nd Sgallari, F., An " III-Conditioned" V olterra Integral Equation Related t o the Reconstruction o f Images from Projections, S IAM J . Appl. Math . 44, 627- 645 (1984). Hamming. R. W., Digital Fill.-rs, 2nd ed., Prentice-Hall, Englewood Cliffs, NJ, 1983. PROBLEMS 4 ,1. F or exact matching o f t he temperatures given below a t x =O.25L in a steel slab, find t he h eat flux c omponents ql, q2. a nd q3 i n W/m2 . 2 2 26.6 3 3 29 28 T he initial temperature is 25 °C; t he slab is 1 cm thick, insulated a t x =L a nd h eated a t x =O; a nd t he t hermal properties a re k =40W/m-C a nd ex = 1 0- s m 2 /s. 4 2. T he surface temperature o fa body, T IM' is t o b e estimated from t emperature m easurements, Y2M a t a n i nterior location o f t he body. Exact matching is t o be used. Derive tiM = tIMlu.=o + ( Y2M - t2M1u.= 0) ( :::) a nd use this e quation t o find t he h eated surface temperature a t t = 1 ,2 a nd 3 s o f t he slab o f P roblem 4.1. 4 .3. F or t he S tolz I HCP a lgorithm a nd for Y 1 d erive Y t-To q l= - - A ¢I 6¢1 q 2= - Y1 - 10) - 2 ¢I A ( . ,.. +To and 11 = To, i = 2, 3,4, . .. , 1 63 PROBLEMS 162 C HAP.4 INVERSE HEAT C ONDUCTION ESTIMATION PROCEDURES tion algorithm for the constant q temporary assumption. G'.= a GjAt Gj = G ( x, t j-At) 2 k' J W hat expression is used for T M+i-IlqM='" =o? 4 .9. Calculate the gain coefficients, K ;, for At+ =0.05, 0.2, a nd 0.5 for the algorithm given in Problem 4.8. T he Green's function for x = L in a plate heated a t x =O a nd insulated at x =L is 1 - A<pI ~ =-----;pr-'''' f o= <PI' a. b. c. Show for measurements a t x =O for a semi-infinite body that fo, flo . .. , f4 are proportional to (At)-1/2. Calculate and plot <Pdi values for i =O, 1 ,2,3, and 4. F or a lumped body, <Pi=iAt/pcL, investigate the / ; values for this case. Compare the results o f parts a and b and give conclusions. 4 .5. Write a computer o r a programmable calculator program for Eq. (4.3.5) a nd verify values obtained for Example 4.1. 4 .6. F or a linear problem with the heat flux history for t>Oapproximated by q (t)=PI + P2 t + . . · + Pr(-I give a whole domain estimation algorithm for estimating P I"'" Pro Use least squares in matrix form for a single temperature sensor with temperatures measured a t n discrete times where n > r. Show the components o f the matrices and let <plj) be the temperature rise for q(t) = tj. 4 .7. Using the least squares method, derive estimators in algebraic form for PI and P2 in Eq. (4.4.18). 4 .8. An alternative integral to Duhamel's theorem is obtained by using Green's functions, T (x, t) = To + ~ E q (l)G(x, t - l)dl where G (x, t) is the Green's function, a =thermal ditTusivity, and k = thermal conductivity. Using the numerical approximation of this equation given in Problem 3.8, derive the following function specifica- G (L,t)=![1+2 L f e-m2,,2allL2(_I)m] m =1 Let a and k be unity. Compare the values with those in Table 4.1. a. Calculate for r = 1, 2, and 3. b. Calculate for r =4. C. Calculate for r = 5. 4.10. Derive a function specification algorithm for estimating g (t) from two interior temperature measurement histories where g (t) is the timevariable volume-energy generation term in a solid cylinder. The differential equation is aT) aT k -a ( r - +g(/)=pcvr Dr al (a) Use the temporary assumption of g(/) equal a constant for r future time steps. The solution of (a) is T (r, I) = To + f~ g (l) VO(r~: - l) d l where O(r, I) is the temperature rise for a unit step increase in g(/) at time 1=0. 4.11. Give the matrix elements for the whole domain regularization equation, Eq. (4.5.18), when Wo= WI = W2 = 1. 4 .12. The backward heat problem is the estimation of the initial temperature distribution in a body knowing one o r more internal temperature histories and the boundary conditions. F or the case of a flat plate with T = To a t x = 0 and L a nd the initial temperature distribution T (x, 0) = 1.. + o 164 CHAP 4 INVERSE H EAT C ONDUCTION ESTIMATION PROCEDURES F(x), derive a zeroth-order whole domain regularization algorithm for estimating n c omponents o f F(x), F i=F(ax- a;), ax=~ C onsider the case o f three interior sensors and m equally spaced time steps. T he describing integral equation is T(X~o + f: F (x')G(x, t, ,x')dx' where G(x, t, x ') is a G reen's function. Modify the notation o f Problem 4.8 t o p ermit multiple sensors. 4 .13. Modify Eq. (4.6.14) for r = 3 a nd calculate the gain coefficients for C HAPTER 5 I NVERSE C ONVOLUTION P ROCEDURES FOR A SINGLE S URFACE H EAT FLUX = 0.05 for the insulated surface o f a h ot plate. Vary a (12 o ver a large range a nd p lot the results. Give conclusions. a t+ 4 .14. Develop . a digital filter procedure for multiple sensors. Give your results in the form o f Eqs. (4.7.6)-(4.7.8). 4 .15. F or the digital prefilter o f 5.1 Y i=aYj-1 + bYj+aYj+1 Show t hat Eq. (4.7 .6) c an be written as M + r-1 qM = L r.<" F M-i(Yj - To)+a[f-rf1(YM+r- To)tl7fM(Y1 - To) i =2 where F M -i=a!M-i-1 + b!M-i+a!M-i+ 1 4 .16. a. Show that Eq. (4 .5.33) reduces t o the Stolz algorithm for a-+O. b. Show t hat Eq. (4.5.33) becomes Eqs. (4.4.25c, d) for a-+ OO A-+'\r;( . Wo= O. I NTRODUCTION In previous chapters, general procedures were given for treating the inverse heat conduction problem ( lHCP) a nd for mathematically modeling the physical problem. The I HCP c an be viewed as the estimation o f t he surface heat flux from transient temperature measurements inside a heat-conducting solid. I t is an ill-posed problem which is characterized by extreme sensitivity o f the surface heat flux to small variations in the interior temperatures. Methods o f reducing this sensitivity were given in Chapter 4. Some o f these methods are used in this chapter. O f t he various ways o f modeling the transient heat conduction in solid bodies, the one used in this chapter employs a convolution integral e quation based o n D uhamel's theorem. This method requires the problem to be linear; t hat is, t he thermal properties (k, p, a nd c) are not functions o f t emperature but can be functions o f position. Numerical methods o f t reating the convolution equations were discussed in Chapter 3. Advantages o f t he Duhamel's theorem approach are that the body can have an arbitrary shape a nd t he thermal properties can be functions o f position (Figure 5.1). T he temperature distribution can be one-, two-, o r threedimensional. T he only requirement is t hat the influence function q,(x, t) be known. [q,(x, t) is the temperature rise a t x d ue to a unit step increase in the surface heat flux a t time zero.] In Figures 5 .la a nd 5.1b the temperature distributions a re functions o f only one spatial independent variable; in Figure 5.1a x becomes x a nd in Figure 5.1b x becomes r, the radial coordinate for a cylinder o r for a sphere. The interfaces between dissimilar materials can have perfect o r imperfect contact characterized by he [see Eq. (1.2.7)]. 1 65 -