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}Ui}UI 4M =OqU.:\t/>jO)2
J
L
+
j =l
( l}.U+li}.u+d 4M =4M+I=OQU.:\t/>jlqU+l.:\t/>jO)2
+ aWo(Qi,+Qi,+d+aW1(qU+IQU)2
(4.5.29)
T he estimator for qu is given by
d lC22d2CI2
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 zerothorder case and Fig. 4.9b for the firstorder 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+lf),u+d 4M =4M+I=O).:\t/>jO
j=1
+I
(4.5.32a)
=O).:\t/>jl
l~r'~rr'~
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 08
1 06
1 04
0.01
...J....l~I
a
FIGURE 4 .9a Gain coefficients for zerothorder 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~rr~''~~
Based o n this limited investigation, the firstorder regularization a ppears t o be
superior t o t he zerothorder 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
100
106
1 0'"
0.Ql
a
FIGURE 4.9b Gain coefficients for firstorder regularization algorithm given by Eq. (4.5 .33).
r =2. Wo=O . W ,=I .
T he z erothorder method has Wo = 1 a nd Wt = 0 a nd t he f irstorder 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 zerothorder (Figure 4.9a), and
the firstorder (Figure 4.9b) m ethods, b ut K I decreases with increasing ex and
approaches zero for t he z erothorder case. N ear ex = 10 6 , t he z erothorder KI
c urve h as a n inflection b ut t he firstorder 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 firstorder 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 erothorder K2 d ecreasing for ex> 1 0 6 ,
whereas t he firstorder K2 r emains c onstant a t 127 for ex> IO 9 .1t is r emarkable
t hat t he firstorder 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 erothorder
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 zerothorder regularization method. This c orroborates the
previous assertion t hat t he z erothorder method reduces the qM values for
large ex's. Since, for ex g reater t han 1 0 6, t he firstorder K I a nd K 2 values approach t he s ame values as in t he function specification method, the firstorder
m ethod is seen to be less sensitive t o t he choice o f ex t han t he zerothorder 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 + rl 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(lB)]TWo[Ho(l B)] + [ HI(IB)yW I [HI(I B)]})q
= X T",  I(y  Tlq=o)+a{[Ho(IB)YWoHoql + [HI(IBWWIHlq/}
(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 firstorder term is d ropped
by setting W I = 0. T hen symbolically, II c an be given by
q =[(12XTX + a(IB) T(IB)r I [(11XT(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 erothOrder R egularization M ethod
The z erothorder 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
twodi~~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>ICX][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 firstorder
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 reformulated 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 online method of analysis.
Heat flux measuring devices can incorporate the digital filter concept and
immediate visual digital output can be provided. (By online 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 +iItM+itlq.,=" , =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"'+r1
t5q2

t5Y)
t5q)
t5q",
t5Y4
t5Y"'+r2
=  = . .. =
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 < jr+ 1
t5lj
Y =1
3
1 2
I I
10
11
12
13
I"
Y I=l
12
1 1
10
11
12
13
t5Y,.
t5q",
t5q2
t5Y"'+r2
Hr
f  I = t5q", = (jqr  I
hY",+1
h Yr
A
(jq",
(jq",
(jq",
q "'=W(YITo )+ ( jy(Y2 To)+ . .. + (jy'
A
q ",=
(jq"' +r I
(jy'
r
( YITo)+
(jq"'+r2
(jy'
r
(jql
( Y2 To)+"'+(jy'(Y",+rI To)
r
+ f l(Y",1 
+ fo(Y",  To)+ f  I(Y",+ 1  To)+ . ..
(4.7.3b)
To)
+ f 2r(Y"'+rl 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 +rl
q ",=
~y.
(jq",
(4.7.2)
"'+r i I (jy"
Note that f lr is not dependent on M but on the difference in the subscripts of
(jq", a nd (jY"'+rI ' Notice also that the first equal sign in each equation of
f l  r+i
( Y"'+r_ITo)(4.7.3a)
I
2
" '+rI
which looks like a Taylor series expansion with the higherorder 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
12
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"'+r1
1 2
1 1
10
11
12
U
L
i= 1
f",  i(llTo )
(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/ mC, 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(3030)29.9(3030)+31.8(30.26930)=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)+. (YJTo)+(Y. To)
UJ
UJ
uYJ
oYJ
= 5 .4(30 30)12.1(30 30) 29.9(30.26930)+ 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 illposed 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 illposed 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 GaussMarkov estimator (Reference 6, p. 232). F or illposed
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{[qME(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 illposed 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,=[qME(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 illposed 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 /mC,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,(OO.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 /mC.
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 IIIPosed Inverse Heat Conduction Problem
I nt.j. H eat Mass Transfer 2 4,1783 1792 (1981).
'
5. O lmsted, J. M. H., Advanced Calculus, PrenticeHall, 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, SpringerVerlag, 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 RADTR76038 ( 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 2HT 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 OneDimensional
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
HeatConduction 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 IIIPosed Problem Arising in Wind
Tunnel Heat Transfer D ata Reduction, Naval Surface W eapons Center, N SWC T R 8232,
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. AIAA840491, 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 Nonlinear 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, PrenticeHall, 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 " IIIConditioned" 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., PrenticeHall, 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/mC 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 tTo
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 jAt)
2
k'
J
W hat expression is used for
T M+iIlqM='"
=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 semiinfinite 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
em2,,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 volumeenergy 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 zerothorder 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=aYj1 + bYj+aYj+1
Show t hat Eq. (4.7 .6) c an be written as
M + r1
qM =
L
r.<"
F Mi(Yj  To)+a[frf1(YM+r To)tl7fM(Y1  To)
i =2
where
F M i=a!Mi1 + b!Mi+a!Mi+ 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 heatconducting solid. I t is
an illposed 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
