SEC . 6.2
C HAPTER
6
D IFFERENCE M ETHODS FOR
T HE S OLUTION OF T HE
O NEDIMENSIONAL I NVERSE
H EAT C ONDUCTION P ROBLEM
6.1

~
I NTRODUCTION
Chapter 3 indicated that difference methods offer considerable potential for
solving a wide variety o f nonlinear direct heat conduction problems. This
statement is also true for the inverse heat conduction problem. In particular,
the nonlinearity associated with temperature dependent thermal properties
can be easily accommodated with difference methods. The recommended
procedure is to locally linearize the problem by evaluating all thermal properties
at temperatures corresponding to the previous time step. Thermal properties
are seldom known to sufficient accuracy to justify iteration. This approach is
called quasilinearization.
Chapter 6 focuses entirely on difference methods for the solution of the
IHCP. The techniques presented are equally applicable to lumped capacitance,
finite control volume, and finite element techniques. Section 6.2 expands earlier
discussions on sensitivity coefficients and how they can be efficiently calculated
by difference methods. Section 6.3 considers exactly matching the calculated
temperatures with measured values for a single sensor for each time step;
Section 6.4 considers mUltiple sensors. Section 6.5 presents techniques for
simultaneously estimating several heat flux components. Sections 6.6, 6.7 and
6.8 apply difference techniques to the function specification method. Section 6.9
considers the sequential regularization method. Section 6.10 introduces space
marching techniques. Section 6.11 gives several examples and Section 6.12
lists some computer programs.
218
2 19
S ENSITIVITY COEFFICIENTS
6.2 S ENSITIVITY COEFFICIENTS A ND THEIR
CALCULATION BY DIFFERENCE M ETHODS
The concept of sensitivity coefficients was introduced in C hapter 1 and has been
used extensively throughout this book. Since the key to understanding the
material in this chapter is a thorough understanding of sensitivity coefficients,
some concepts about them will be expanded. Two sensitivity coefficients occur
repeatedly in this chapter: (1) that due to a step qM in heat flux for infinite time
duration and (2) t hat due to a pulse qM in heat flux for time duration t M  tM  I '
B~th the heat flux time variation and the corresponding sensitivity coefficient
variation with time are shown schematically in Figure 6.1. Actual step function
sensitivity coefficients are shown in Figures 1.7 and 1.11 for planar slabs with an
insulated surface and slabs that are semiinfinite, respectively. Since the duration
of the step function is infinite, the step function sensitivity coefficient grows
without bound. The curve labeled A in Figure 6.1 is for a step that begins at
time t MIo t hat labeled C is for the same step but beginning at time t M' The
symbol Z is used to denote the step function sensitivity coefficient; this Z is
the same as the 4> in Chapters 3 5 i fthe thermal properties are constant. F or a
sensor a t depth X l, the temperature response due to the step in heat flux can
be written as T(Xko t , t M  I , q MI, q M); t M  I denotes the time a t which the
heat flux step o r pulse begins. The components of the heat flux vector qM  I ,
which are <1 I, <12' ... , <1M  I' contain al\ the information about the time variation
of the heat flux prior to tM  10 and hence contain all the (initial) temperature
information necessary to continue the temperature calculations to t M , tM + I , • • • .
The step function sensitivity coefficient is defined by
OT(Xko t , tM 
10
11M  I, q M)
(6.2.1)
OqM
aT
aq
q (l)
1 .,1
I .,
1 .,+1
Time
FIGURE 6.1 Schematic o f sensitivity coefficient for step in heat ftux a nd pulse in heat ftux.
A, sensitivity coefficient for step in heat ftux at time 1 .,_ , . B, sensitivity coefficient for pulse in heat
ftux, beginning at time 1., _,. C , sensitivity coefficient for step in heat ftux at time I.,.
SEC. 6.2
220
C HAP.6
221
SENSITIVITY COEFFICIENTS
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
profile at tM  1 • Using Eq. (6.2.2) and the notation of Eq. (6.2.4), one can show
that for constant properties,
is sensor location
is arbitrary time
is the time at the beginning of the heat flux step
is the previous heat flux history
is the magnitude o fthe heat flux step
where:
r
Z l,r=
1,
q M_.)=Z(Xl, t, t Ml' ~I)
 Z[Xlo t (tM tM 1 ), t Ml, q Mt]
which simply states that for the linear problem the response due to a step is
the sum of the response due to a series o f pulses of constant magnitude and
distributed over the same time interval.
F or nonlinear problems, the sensitivity coefficients can be calculated from
differences of the temperatures at the same time and location for two values
of the heat flux. F or example, the temperature distribution is calculated for two
different values of heat flux, q . and (1 + e)q·, where e is a small parameter o f
ord~r 0.001. Then the sensitivity coefficient is given by
X ( Xl, t, t M 
A
T.=.[x...::.:.,...:.t,...:tM
=.:,:I,...:q=::M:......:I.:.'q..:...·....:.(_I_+'e):,=.]__T_(:...X::..:1':...t:..'.=M:......!.I.:.'q..=:M:!.._I!..:.'...!.q....:,·)
t
1
eq·
) __
1 , qM  1 
(6.2.5)
The only difference in the X and Z sensitivity coefficient calculations is the
boundary condition:
Sensitivity Coefficient
Boundary Condition
X(pulse)
q(t) = { q.,
Z(step)
q (t)=q·,
(6.2.2)
The two terms on the right of Eq. (6.2.2) are shown in Figure 6.1 as curves A
and C ' for constant thermal properties, curve C is identical to curve A except
it is di~placed in time by an amount tM tM  1 • Curve B is ( AC) and is the
pulse sensitivity coefficient X which is analogous to AljJ o f Chapters.3 to 5.
When using difference methods, discrete values of X and Z are Important.
At time t =tM'
(6.2.3)
because Z (Xl' t M 1 , t M 1 , qM .)=0. F or the constant property case, a convenient subscript notation can be used.
Z l.jM+l = Z(Xl, t j t M
X1,j
j= 1
F or a linear problem, Z is independent of qM; that is, T is linear in qM' There
are applications in which the thermal properties are held constant .ov.er a few
(typically 1 10) future time steps while allowing some property ~anatlon over
the entire problem time. F or this case, Z is independent of qM b ut It does depend
on t  and qMl through the temperature profile at t Mt· I f the thermal
M1
properties are independent of temperature over the entire problem time, then
Z is independent of qM a nd q Ml and the important time var~a~l.e is t  t M. 1 ,
the time from initiation of the step. These three important sensltlVlty coefficient
cases are summarized in Table 6.1. By definition, a sensitivity coefficient is
identically zero prior to the initiation of the heat flux step o r pulse.
The pulse sensitivity coefficient will be denoted by X with th~ s~me arguments
as those for the step function sensitivity coefficient. F or quaslhnear and constant property cases, superposition can be used to calculate X from Z.
X (Xl' t, t M
L
(6.2.4)
1)
This same notation will also be used for the quasilinear property case where it
is implicitly understood that all properties are evaluated at the temperature
0,
Case
II
I II
None, general case
Quasilinear thermal properties
Constant thermal properties
 k a T!

=q(t)
. <=0
az = ~ (k az)
at
= qdt)
. <=L
T (x, tM 
d
(6.2.6)
 k a z!
=1
ax . <=0
ax
ax
qMI'
qM
Z (Xh t , t MI'
q M)
Z (Xh t , t M1t
d
Z (Xh t t M  d=4>(Xh t t M
t~tMl
Sensitivity Coefficient, Z
 k a T!
Notation
~t~tM
Temperature, T
N otation f or S tep F unction S ensitivity C oefficients
Restrictions
1
The evaluation of Eq. (6.2.5) requires two different temperature calculations
(or direct solutions). I t is computationally more efficient to calculate X directly
from its own partial differential equation, which can be derived from the
transient heat conduction equation. F or example, consider the problem discussed in Section 3.3 for the quasilinear case· for t > tM  1 ;
pc
T ABLE 6.1
tM 
all other t
d = F (x)
ax
i kJZ!
aX
ax
(6.2.7)
=0
. <=L
Z (x, t MI)=O
· The quasilinear case implies that k, p, a nd c are evaluated for the temperature profile at
1 ",1'
222
CHAP. 6
ONE DIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
A comparison o fthe above two sets of equations for T a nd Z indicates that they
both satisfy the same partial differential equation. The boundary and initial
conditions are different and those for Z are considerably simpler. If X is to
be calculated instead o f Z, then the foregoing x =o b oundary condition must
be modified as follows
 k OXI
= {1 for pulse duration
o x .. = 0
0 for all other times
(6.2.8)
The differential equation and the remaining boundary and initial conditions
for X are identical to those for Z. D ue to the similarity o f the T a nd Z differential
equations. considerable computational savings are possible. This point is
explored in detail in a subsequent section.
6.3 S INGLE TEMPERATURE SENSOR, FUNCTION
SPECIFICATION ( q=C), S INGLE F UTURE T IME STEP
( EXACT M ATCHING O F D ATA)
A single temperature sensor is considered to be located at a depth X l below
the active surface. If the heat flux qM is constant over the time interval tM  I ~
t~tM' the value of qM t hat forces a matching of the computed temperature at
X l with 'the measured temperature can be calculated. F or realistic temperature
d ata containing errors and small dimensionless times, this approach produces
significant oscillations in the computed heat flux and would probably not be
used in practice. The method is studied here because it aids in the understanding
of inverse problems and also because an exact matching of temperature data
is a building block that can be used in the development of other techniques.
6.3.1 M odification o f D ifference E quations o f t he D irect
H eat C onduction P roblem f or t he S olution o f t he I HCP
Difference methods for solving direct heat conduction problems were introduced
in Section 3.3. By analogy with Eq. (3.3.45), the difference equations for a onedimensional direct heat conduction problem with a nonuniform grid can be
written in symbolic form as
SEC . 6.3
2 23
SINGLE T EMPERATURE SENSOR
The coefficients ai' bi , Ci ' ai, bi, a nd ci depend on grid spacing and thermal
properties a nd a re independent of the constant val~e o f ~eat flux qM for the
linear and quasilinear cases. F or the purpose o f diSCUSSion, only five nodes
are used in what follows. Equation (6.3 .1) can be written more compactly as
bl
 C2
o
o
o
al
b2
 c]
0
0
0
 a2
0
0
0
b]
 a]
 C4
0
b4
0
0
 a4
C s
bs
Tf
T 't
d; +glqM
d2
T~
d]
T~
d4
ds
T 't
(6.3.2)
The d's represent all of the righthand side o f Eq. (6.3.1) except the. qM term. In
the direct problem the T~,j=I, . .. , 5 values are found from a Simultaneous
solution of Eq. (6.3.2). N ote the explicit dependence of the temperature o n the
heat flux qM'
.
Instead o f solving the direct problem given by Eq. (6.3 .2), suppose qM IS a n
unknown and a temperature sensor is located at node 3 ( x./L=O.5). T he direct
problem unknown T~ is now a known measured temperature YM • F or t he
IHCP, Eq. (6.3.2) becomes
bl
 C2
0
0
0
 al
b2
 C]
0
0
gl
0
0
0
0
0
0
 a]
b4
 Cs
0
0
0
 a4
bs
Tf
T 't
qM
T~
T 't
d;
d 2+a2 YM
d ]b]YM
d4 +C4YM
ds
(6.3.3)
Equation (6.3.3) represents a system of five linear algebraic equations for .the
five unknowns ( Tf, T 't, qM, T~, T 't) a nd can be solved by any approprIate
technique for linear algebraic equations. As indicated in previous chapte.rs,
this formulation is linear in the unknown heat flux a nd s hould n ot reqUIre
iteration. I t h as been demonstrated by Beck and Wolfl t hat the foregoing
formulation is unstable for 0 = 1/2 a nd lumped thermal capacitance. Hence,
when the d ata from a single temperature sensor are matched exactly, 0 = 1/2
should be avoided. Alifanov l6 used the same technique with 0 = 1. (See Sect.
3.3.3 for definition of 0.)
6.3.2 S ensitivity C oefficient A pproach f or E xactly
M atching D ata F rom a S ingle S ensor
(6.3.1)
In Section 6.3.1 it was demonstrated that matching d ata exactly from a single
temperature sensor produces a system of linear algebraic equations. ~~wever,
the system of equations represented by Eq. (6 .3.3) is n o longer of the trIdiagonal
form. (See Section 3.3 .4.) Hence, a direct solution of Eq. (6.3.3) is n ot the most
computationally 'efficient approach. Because o f the linearity of the I HCP, it is
possible to utilize the tridiagonal structure of the algebraic equations (for onedimensional problems) and develop a very efficient algorithm.
2 24
C HAP.6
O NEDIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM
The temperature field T (x, t) depends in a continuous manner on the unknown heat flux qM (constant over the interval t M I ~t~ tM)' This dependence
is written as T (x, t, tM_ I , qM _ I , qM) where qM _ I is the vector of all previous
heat flux values and tM  I indicates the time that" the heat flux step begins.
Because the temperature field is continuous in qM, it can be expanded in a
Taylor series about an arbitrary but known value o f heat flux q *,
T (x, t, t MI' qMI , q M)= T (x, t, t M I , ~lo q *)
+ (q Mq
S EC.6 .3
known, the complete temperature field can be calculated from Eq. (6.3.5).
Equation (6.3.8) is equivalent to the Stolz equation given by Eq. (5.3.1) if q * =0.
An efficient way o f calculating sensitivity coefficients is to solve the differential equation for X that is analogous to Eq. (6.2.7) . Since the T and X
calculations are similar, considerable computational savings are possible. The
development o f a n efficient algorithm requires a thorough understanding of
the following tridiagonal matrix algorithm (TDMA):
Linear eqUIJtion lorm :
b l ula l u2=d l
*) o T(x, t, t M1! qMI ' qM)\
OqM
t M='.
 Cjuj_l+bjujajuj+l=dj , j =2,3, . .. , Nl
+ ( qMq*)2 o2T(x, t, tM  I , q Mlo qM)\
+
2!
o ql,
t M=,.
...
(6 .3.4)
Forward elimination (starting a t node N ) :
1
WN= b ' eN=WNcN,
N
T (x, t M , t M1! qMl> q M)= T (x, t M , tM I , ~I' q *)
W j=b
(6.3.5)
'
t
M,
M I' MI 
!l
uqM
'M='·
(6.3.6)
Note that X does not depend on qM for a linear problem.
Equation (6.3.5) is a prescription for calculating the temperature field
corresponding to qM provided the temperature field is known for a given q*
a nd a means exists for calculating the sensitivity coefficient X. I f the problem
is amenable to an analytical solution, X can be calculated as indicated in
Section 1.6. Difference methods for obtaining X are discussed in Section 6.2
in connection with Eq. (6.2.5).
F or the I HCP in which the sensor data are matched exactly, the lefthand
side ofEq. (6.3.5) is replaced by the experimental temperature; hence, Eq. (6.3.5)
becomes
y M=tt' + (qMq*)X t . 1
(6.3.7)
Solving Eq. (6.3.7) for the estimate of the unknown heat flux, qM, gives:
A
q M=q
IN=wNdN
.
J =Nl,N2, . .. , 1
e j=wjcj , j =Nl,N2, . .. , 2
(6.3.10)
I I = wI(d l + ad2)
) _OT(X,tM,tM  I,qMI,qM)\
q
1
'
j ajej+1
J j=wj(dj+aJj+d,
where the sensitivity coefficient is defined by
X (x t
(6.3.9)
CNUNI +bNuN=dN
F or linear problems, only the first derivative in Eq. (6.3.4) is nonzero, thus the
following is a n exact result for location x a t time tM
+ (qMq*)X(x, t M , t M I , q Md
226
SINGLE TEMPERATURE SENSOR
* + y Mtt'
X
(6.3.8)
t .1
where t t' is the temperature at tM for the sensor node k with q = q * over tM  I ~
t ~ t M • T he calculational procedure is to assume an arbitrary value of heat
flux q *, calculate the temperature field T (x, tM, tM lo q MI' q*), and knowing
the sensitivity coefficients X (x, tM, t M lo qMI), calculate the heat flux qM
from Eq. (6.3.8) t hat exactly matches the temperature data YM • Once qM is
Back substitution :
u l=il
u j=ejuj_I+Jj, j =2,3, . .. , N
(6.3.11)
In comparing the two problems for T and X, the following conclusions can
be drawn:
T
aj' bj ,
Wj' ej
Cj
X
+same for b oth+same for both because t hey_
depend only o n aj' bj, C j
aj, bj ,
Wj ' ej
Cj
dl = gl
dj =O,j=2, 3, . . . , N
d l =d'i + glqM
dj"/=O,j=2, 3, . .. , N
Note that the difference equations for X can be derived by differentiating the
difference equations for T with respect to qM ' G oing through the forward
elimination portion of the algorithm for X, it can be shown that all the J.'s
are identically zero except for I I a nd
•
I I = glwi
(6.3.12)
Consequently, the sensitivity coefficients are calculated as follows :
X I.I = /I=glwl
(6.3.13)
226
C HAP.6
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
SEC. 6.3
X j.l=ejXj _ I.I,
(6.3.14)
j =2,3, . .. , N
The algorithm for exactly matching the temperature data from a single
sensor can be summarized as follows:
1. Assume an arbitrary value o f q* and calculate the entire temperature
field T =T(x , tM ' tM1> q MI, q*). q *=qM1 is a common choice; however,
q* = 0 is also perfectly acceptable for the linear problem and will reduce the
total number of calculations.
2. Using the same matrix coefficients a n, bn , Cn , W n , and e n as were used for
the T calculation in step 1, calculate the sensitivity coefficients from Eqs. (6.3.13)
a nd (6.3.14)
(6.3.13)
(6.3.14)
3. Calculate the heat flux t hat exactly matches the experimental temperature d ata by using the Taylor series expansion
OM
l
q M=q * + YM  T
X
( 638)
..
A
227
SINGLE T EMPERATURE SENSOR
N a , N m , and Nd represent the cycle times required for a single addition, multiplication, and division, respectively. Using the information in Tables 6.2 and
6.3, the approximate fractional increase in computation time for the inverse
problem relative to the direct problem is given by:
~T
2 (N + l )+ 2N(Tm/Ta) + (Td/Ta)
T
3(N  1)+5(N 1)(Tm/Ta)+N(T,,/To)
"Computation index" is probably a more appropriate term for ~T/T because
Eq . (6 .3.16) does not account for the fact that some cycle time is lost between
the end of one operation and the beginning o f the next operation. F or a particular computer, ~T/T depends only on the total number o f nodes N. Using
the computation times given in Table 6.3, ~T/T for various values of N are
presented in Table 6.4. The I HCP causes less than a 40% increase in computations in comparison to the direct problem.
Note that Eq. (6.3.16) does not depend on the node number of the temperature
sensor; a temperature sensor a t node 1 requires the same number of computations as a backface sensor at node N. Some computational improvement is
l .1
4.
sion
Calculate the complete temperature field from the Taylor series expan(6.3.15)
Table 6.2 presents an arithmetic operation count to determine the relative
increase in computations o f the inverse problem over the direct problem. The
operational counts can be used to evaluate the computational efficiency for a
particular algorithm. Table 6.3 presents the number of cycles on the C DC 7fIXJ
computer required for the four arithmetic operations; one cycle is 27 ns. Let
T ABLE 6 .2 O perational C ounts f or O ne T ime S tep o f I nverse
A lgorithm t hat E xactlv M atches S ingle T emperature S ensor D ata
( N=total n umber o f n odes)
E quation
+ o r
x
(6.3.10)
(6.3.11 )
S tep
3(Nl)
S (NI)
C omments
N
D irect
calculation
T ABLE 6 .3 R elative S peeds o f V arious
A rithmetic O perations o n t he C DC
7 600 C omputer 2
O peration
N o . o f M achine Cycles
Relative T ime ( f/r.)
4
5
20
1.0
1.25
5.0
+ o rx
T ABLE 6 .4 C omputation
I ndex f or E xactlv
M atching D ata f rom
S ingle T emperature S ensor
( for C DC 7 600 c omputer.
T able 6 .3). N o. o f
e lements=N1
N
2

~
~
(6.3.13)
(6.3.14)
3
(6.3.8)
2
4
S um(2+3+4)
(6 .3.15)
2N
2 (N+ I )
4 .5N+7
t n/f=    14.25N  9 .25
11
21
31
41
81
0 .383
0.350
0.339
0.333
0.324
N
N
2N
Inverse
calculations
only
(6.3.16)
228
CHAP. 6
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
possible provided the algorithm is restructured. This restructuring is possible
because it is not necessary t o k now the sensitivity coefficients X i , l a nd t~ for
x > X x (the sensor depth). I f t he sensor is n ear node 1, t he active surface, a restructuring c an offer considerable savings. T he c omputational algorithm is
r eordered in the following manner:
1. Assume a value for q* a nd perform the forward elimination portion
o f t he tridiagonal matrix algorithm. (Note t hat t he only calculation in this step
t hat d epends o n t he heat flux is f l ')
2. C alculate t~ from the back substitution portion o f t he tridiagonal
matrix algorithm.
t il = fl
OM
T J =eJ t MI+fj,
J
j =2,3, . .. , K
(6.3.17)
X i,l
229
SINGLE TEMPERATURE SENSOR
T ABLE 6.5 O perational C ount f or R estructured A lgorithm t hat
E xactly M atches S ingle T emperature S ensor D ata ( K=sensor n ode
n umber. N =number o f n odes. K :r;,N)
Step
E quation
1
5
1 +5
3
(6.3.13)
(6.3.14)
(6.3.8)
C omments
x
+ o r
(6 .3.11)
(6.3.10)
O ne equivalent direct problem
solution
5 (N1)
3 (NI)
N
;(6~;. ;;) ~ ~; ~ ~I t~ c~~~~I~;i~~           4
t (xx, t M)=t:
3. C alculate the sensitivity coefficient
the previous algorithm
SEC . 6.3
5
Sum
(inverse)
K
X . ,I calculation
2
I I = W.(dl +a.J2)
E xtra inverse calculations
2
K +2
2K + 1
in exactly the same way as for
(6.3.13)
(6.3.14)
4. Calculate the heat flux t hat exactly matches the experimental d ata by
using the Taylor series expansion,
(6.3.8)
5. With the heat flux qM calculated from step 4, recalculate f l = wl(d l + ad2)
a nd repeat the back substitution portion o f t he tridiagonal matrix algorithm.
T~=fl
(6.3.18)
Figure 6.2 presents the results of Eq. (6.3.19) for the C DC 7600 computer. I f
the temperature sensor is located near the active surface (small numerical
value o f K), t hen the restructured algorithm is very efficient a nd is faster t han
the originally proposed algorithm.
The foregoing type o f minimization o f calculations for the I HCP was first
developed by Blackwell.3 However, the original approach did n ot m ake explicit
use o f t he sensitivity coefficient concept.
I n realworld problems, there will be errors in the temperature measurements. I f t he temperature d ata a re matched exactly, the temperature errors
produce heat flux errors. Consequently, exact matching o f t he temperature d ata
0 . 40r.....~~~,
~I. 0.30
'
I
!' ~ .M._
~ ..
N ote t hat it is n ot possible t o calculate the temperature field from the Taylor
series result, Eq. (6.3.15), because the sensitivity coefficients for depths greater
than the temperature sensor depth were n ot calculated in step 2.
An operational c ount d emonstrates that the restructured algorithm is m ore
efficient than the one first proposed. Table 6.5 summarizes the operational count
for the new algorithm. N ote t hat when steps 1 a nd 5 a re combined, they are
equivalent t o o ne direct problem solution plus one additional evaluation of
f l = wl(d l + ati2)'
/
/
~T For ~f.W!i~5
/
00
20
40
60
80
K, Node number of temperature sensor
T he c omputation index for the restructured algorithm is
M
t

8 0.10
K + 2+(2K + 1)(r",/tQ)+(td/tQ)
3 (N  1)+ 5 (N  1)(t",/tQ)+ N(td/tQ)
FIGURE 6 .2 C omputation index for algorithm that exactly matches single temperature sensor
data ( K = sensor node number, N = number o f nodes, K ~ N ) .
(6.3. 19)
230
C HAP. 6
ONEDIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM
is n ot recommended unless the temperature data are smoothed prior to calculating the heat flux and/or relatively large time steps are taken (a.lit/xi. > 1).
The function specification method using future temperatures (Chapter 4) can
be used with difference methods and allows smaller computational time steps
to be used without encountering stability problems; this method is discussed
in Sections 6.6 and 6.7.
6 .4 M ULTIPLE T EMPERATURE SENSORS. FUNCTION
SPECIFICATION ( q=C). SINGLE FUTURE T IME STEP
Multiple temperature sensors are recommended in order t o o btain as much
experimental information as possible. These include sensors at different depths
as well as several sensors at the same depth. F or multiple sensors, only one
sensor can be matched exactly. An alternative is t o determine the value of qM,
c onstant over the interval tM  1 ~ t~ tM • t hat minimizes the least squares error
for all sensors considered. I f there is a total of J sensors, then the least squares
error is defined as
SEC . 6.4
M ULTIPLE TEMPERATURE SENSORS
I f J = 1, then Eq. (6.4.5) reduces to Eq. (6.3 .8). Again, note that q* = 0 is a valid
assumption.
The computation index in Figure 6.2 is very nearly valid for the multiple
sensor case provided K is the node number of the sensor located the greatest
distance below the active surface. The only additional calculation of the multiplesensor problem compared to the singlesensor problem (with the sensor located
at node K ) is the evaluation of Eq. (6.4 .5) instead of Eq. (6.3 .8). This additional
calculation is generally trivial in comparison to the total number of calculations. Hence, Figure 6.2 is a n approximate indicator for the multiple temperature sensor case.
An alternative to the use of Eq. (6.4.5) for the multiplesensor problem for
a single future time step is to calculate the heat flux Zit .u t hat exactly matches
the individual (kth) measured temperature. > t,u, a nd then take a weighted
average o f all Zit.M'S . This average can be derived by using a Taylor series slightly
different from Eq. (6.4.4). Expanding the temperature field for q u in terms of
ilt.M
T(xto t M , t "'_I, q UI' q u)= T (xt, t u, t Ulo q UI, Zit.u)
_ ) oT(Xtotu.tul,qUI,qu)1
+ (q uqt.u
0
qu
9. . 4....
(6.4.6)
J
s=
L
t =l
[>t . MT(x to t u ,tU  lo qUloqU)]2
(6.4.1)
where k is the sensor index and M is the time index in >t . M'* I t is i mportant to
remember that k is the sensor index, not the finite difference node index. Differentiating Eq. (6.4.1) with respect to qM, gives
oS
 0 = 0=2
qu
J
L [ >t.uT(Xt, tM.tU 
t =l
1•
q UloqU)][X(xtotU,tU1,qUI)]
(6.4.2)
But exactly matching the kth sensor at t u gives
T (xt. t u, t UI' q Ulo Zit.u) = >t.M
T (xt, t u, t UI' q UI' q",)= f r + (quq*)Xt . 1
(6.4 .4)
Substituting Eq. (6.4.4) into Eq. (6.4.2), replacing q u by q u and solving for
q u, gives
J
L
X t.I(>t.ufr)
t
q u = q* + c..==..I . J    
OT(Xl' t u , t"' _ I' q Ul, q u )1
o
= X(xl , t u, t UI, q ud=Xt.1
qM
9. .  4. ...
The minimum least squares condition, Eq. (6.4.2), becomes
L Xf.l
t =1
° The double subscript notation f (x, . t ",)= Yo.", corresponds t o the classical (x. t) notation.
 
(6.4.8)
J
L X t.l(quZit.u)Xt . 1= 0
t =1
(6.4.9)
Solving for q u,
J
q u=
L WtZit .u
t= I
(6.4.10)
where the heat flux weighting factors are related t o the sensitivity coefficients
through
Xl. I
(6.4.5)
(6.4.7)
and
(6.4.3)
As in Section 6.3, a n arbitrary value o f heat flux q* is assumed and the Taylor
series expansion given by Eq. (6.3.5) is applied:
231
J
(6.4.11)
L x li
t =1
Equation (6.4.10) provides an alternative interpretation to this solution of
the multiple sensor, single future time step I HCP. T he procedure is to treat
232
C HAP.6
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
each sensor independently (over a single time step) and determine the heat
flux qt.M that exactly matches the single data point l't.M ' This. heat flu~ is
related t o the assumed arbitrary heat flux q* through the Taylor !.enes expansion
~
_
qt,M=q • + l 't,Mf: , k 1,2, . .. , J
After evaluating Eq. (6.4.12), these qt.M heat flux values are weighted (averaged)
according to Eq. (6.4.10). F rom the discussion of sensitivity coefficients in
Section 1.6, the sensors closest to the active surface will have the largest sensitivity coefficients and hence will automatically be weighted more heavily. F or a
problem with constant thermal properties, the sensitivity coefficients X t,l
a nd hence heat flux weighting coefficients Wt need be calculated only once. I f
the thermal properties are treated as being quasilinear in temperature, both
Xt • 1 and Wt must be recalculated a t each problem time step.
6.1 . Using the information in Table 1.1 for the response o f a p lanar slab
with a n insulated inactive surface to a step in heat ftux, calculate the heat ftux weighting
coefficients o f E q. (6.4.11) for two sensors a t dimensionless depths o f x t!L=O a nd
X2/L=0.s a nd for two dimensionless time steps o f d t+ = alit/L2 = 0.05 a nd 0 .5. N ote
t hat T able 1.1 is also shown in graphical form in Fig. 1.7.
E XAMPLE
Solution. Although Table 1.1 c ontains the dimensionless temperature response, the
dimensionless sensitivity coefficients iJT/iJ(qL/k) a re identical t o T+.
a lit
0 =0.5
X t,l = 0.252313
X i,I =0.458333
W I =0.7671
w2=0.2329
response for t > tM' One approach that makes use of all tempera.tur~ informat.ion
simultaneously is the whole domain estimation procedure which IS the subject
of the next section.
6.5 W HOLE D OMAIN E STIMATION W ITH DIFFERENCE
M ETHODS
The whole domain procedure has been discussed in Sections 4 .4.2. This procedure is most applicable when the physics o f the problem suggests that the
heat flux variation with time is o f a particular analytical form, for example,
forms such as given by Eqs. (4.4.1) (4.4.5). This list is by no means exhaustive.
A definite advantage o f the whole domain estimation procedure is t hat stability
problems are minimized. However, more computations are required.
.
The whole domain procedure will be introduced first by means o f a speCific
example and then will be generalized, Suppose a parabola adequately represents
the time variation o f the heat flux,
(6.5.1)
The parameters P I' Pz, P3 are to be estimated so that the computed temperatures agree (within certain limits) with the experimentally measured temperatures. The extent of the fit o r agreement is determined by the ordinary least
squares criteria. I f Yi is the experimentally measured temperature and 7; is the
predicted temperature, both at time t ; a nd a t the same location, then the least
squares error is given by
X t,l =0.831876
X i,I =0.0153659
W I = 0.9963 (x/L=O)
W 2 = 0 .0037 ( x/L=0.5)
233
W HOLE D OMAIN E STIMATION WITH DIFFERENCE METHODS
(6412)
..
X t,l
+ a lit
d t = 2 = 0.05
L
SEC. 6.5
r
S=
L (Yi _7;)2
(6.5.2)
i =1
o
The predicted temperatures depend implicitly on the parameters (P I , P2, P3)
7;=7;(PI,P2,P3)
In comparing the W t'S for the foregoing example, several important conclusions
can be drawn. F or at+ = 0.05 the sensor at x+ = 0.5 yields very little information
because its sensitivity coefficient is small relative to that for the surfacemounted
sensor. F or at+ = 0.5, the time step is sufficiently large that sensor 2 has time to
respond; however, sensor 1 is still weighted more heavily because o fthe square
of the sensitivity coefficients.
Because temperature measurements invariably contain errors, the heat flux
components determined by exactly matching a single sensor contain amplified
errors. Even for the multiple sensors considered in Example 6.1, sensors located
far apart tend to olTer very little smoothing o f the heat flux; consequently,
other approaches should be considered. Smoothing o f the temperature data
prior to the solution o f the I HCP certainly reduces the oscillations in the heat
flux; however, this approach does not make use of the fact that a pulse of
magnitude qM over the interval tMI~t~tM has an influence on the sensor
(6.5.3)
The object is to choose (Ph Pl, P3) such that the ordinary least squares error S
is minimized. This condition is determined by
oS
r
07;
.
(6.5.4)
L (Yi1;) U!>pj = 0, J = 1 ,2,3
oPj
; =1
Assuming that the computed temperature is continuous in the parameters
(PI' P l, P3), the temperature can b e expanded in a Taylor series about trial
values o fthe parameters. I f (PI. Pl, pj) represent the trial values of the parameters and (P~+ I, Pl+ I, pj+ I) represent new o r improved estimates o f the parameters, then
 = 2
t ;+I=t ; +(Pi+IPD07;1
OPI
,=, +(P'2+
I _Pl)o7;1
OP2
,=, +( Pj+l_p;)o7;1,=,OP3
(6.5.5)
_ _ _ __ 0
_. _ _ _ _ _ _ • _ _ _ _ _
SEC. 6.5
234
C HAP.6
W HOLE D OMAIN E STIMATION WITH DIFFERENCE METHODS
235
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
sensitivitycoefficient (11) differential equation is
Note that higherorder terms are ignored in the Taylor series expansion;
however, the expansion is exact for linear problems. Substituting Eq. (6.5.5)
into each o f the three equations represented by (6.5.4) and simplifying,
±OoT; OoT; ±OoP2 OoT; ±aopJ OoT; • IW'Pi
T;
T;
PI PI
PI
PI
±OoT; OoT; ±OoP2 0oT; ±ooPJ OoP2 P't Pi
T;
T; T;
PI P2
P2
±OoT; ooPJ ±OoP2 0oT;
T;
T;
Pt
PJ
j= I
i= I
ok OTI
 1I1"T"
u
u x x =o
011 / 1
 k"
u x x =o
= t j  I jl,2,3
, ·_
(6.5.8)
i= I
I
i ='
j =.
i= I
i ='
_ 1I
i= I
r
oT.
r
oT.
i =1
OP2
r
oT.
L (~1;)'
i~1
OPI
L (~1;)'
OkOT
I
_k01l/1 = 0
loT ox x =L
ox x =L
j =1,2, 3
v
(6.5.10)
(6.5.11)
(6.5.6)
L (~1;)'
OP3
i =1
The superscripts v o n the coefficient matrix and the righthandside vector
indicate that all terms should be evaluated a t P= p'. Equation (6.5.6) is a linear
system of three algebraic equations with three unknowns; P't I  Pi is treated
as an unknown correction factor. F or a linear problem, the sensitivity coefficients
are independent o f the parameters, the coefficient matrix is in turn independent
of parameters, and only one step (iteration) is required. Even though the heat
flux in this example is linear in the parameters (PI, P2, PJ), temperaturedependent thermal properties make the problem nonlinear and require iteration.
F or the nonlinear problem, one approach to the sensitivity coefficient calculation is to use differences. F or example,
(6.5.7)
where e is a small number, say 0.001. The evaluation of each sensitivity coefficient requires the solution o f two direct problems but one is the standard
condition; for the foregoing example with three parameters, a total of four
direct problems must be solved for each iteration in Eq. (6.5.6). Each iteration
will also require the simultaneous solution o f three linear algebraic equations,
Eq. (6.5.6). (See the comments below (4.5.19).)
At first glance, it might appear that some computational savings are possible
if the sensitivity coefficients are calculated directly from a differential equation
similar to Eq. (6.2.7). However, this is not the case for the nonlinear temperaturedependent thermal property case. F or the whole domain problem, it may not
be possible to locally linearize the problem by evaluating the thermal properties
at the previous time step as it was for the sequential procedure o f Section 6.3.
Considering Eq. (6.2.7) as the basic model for the temperature field, the
This sensitivitycoefficient differential equation is linear in 111 if the temperature field is known because terms like k, pc, ok/oT, and O(pc)/oT can be treated
as known functions of position. However, the differential equation for 111
is no longer the same as that for T. Although the parameters (PI' P2, P3) d o not
explicitly appear in Eqs. (6.5.8)  (6.5.11), they appear implicitly in the predicted
temperature field T(P., P l' PJ). I f the heat flux model is nonlinear in the parameters, then the parameters appear explicitly in the sensitivity coefficient boundary condition. F or example, using Eq. (4.4.4) as the heat flux model,
(4.4.4)
then the boundary condition for
11, is
11
(111 oT Ox x=o  (k 0ox ')1 x=o =Plte fl ,'
ok T)I
o
(6.5.12)
which is nonlinear in the parameters.
The foregoing procedure can be readily extended to heat flux models with
a large number o f parameters. In such a case, it is beneficial to use matrix
notation; additional details on matrix calculus can be found in Chapter 6 of
Beck a nd Arnold. 4 Let Y be the vector of the discrete temperature measurements and T be the vector of the corresponding computed temperatures at the
given sensor location,
(6.5.13)
The subscripts on T and Y denote time and there are r discrete times. The
matrix form o f the least squares function, Eq. (6.5.2), is
(6.5.14)
236
C HAP 6.
O NEDIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM
T he vector o f p parameters,
I, is
SEC . 6.6
237
SINGLE T EMPERATURE SENSOR
where
(6.5.24)
(6.5.15)
(with p <r) a nd it is t o be determined so t hat S is a minimum. This procedure
will involve the matrix derivative o f Eq. (6.5.14); this derivative is defined as
a
a pi
VI=
a
ap2
(6.5.16)
6.6 SINGLE T EMPERATURE SENSOR, F UNCTION
S PECIFICATION ( q=C)" F UTURE T IME STEPS
a
app
Differentiating Eq. (6.5.14) with respect to
F or a linear problem, X (I) is i ndependent o f IJ a nd only one iteration o f Eq.
(6.5.22) is required. I f t he sensitivity coefficients are calculated using differences
as in E q. (6 .5.7), a p parameter problem will require p + 1 direct problem solutions for each iteration. F or p large and approaching r, as was the case for the
sequential method o f Section 6.3, the number o f calculations required for the
whole domain estimation becomes very large, even for the linear problem.
Consequently, the whole domain estimation procedure is n ot recommended
except when the functional form o f t he heat flux model is strongly suggested
by the physics o f t he problem a nd the number o f p arameters is reasonably small.
I,
VIS =  2(VI TT)(y  T)=O
(6.5.17)
T he sensitivity coefficient matrix is defined as
X (I) = (V, TT)T
(6.5.18)
or,
I

aTI a T
ap, ap2
a TI
app
aT, oT,
api ap2
aT,
ap p
(6.5.19)
X=
tr
Substituting Eq. (6.5.18) i nto Eq. (6.5.17),
X T(YT)=O
T he discussion in Section 1.6 p ointed o ut t hat the sensitivity coefficient associated with a heat flux pulse o f finite duration can be nonzero after the heat flux
returns to zero. Consequently, q(t)=q"" over the interval t"'_I~I~t", a nd
zero otherwise, will influence sensor measurements Y"" Y",+ I ' . ... Therefore,
powerful procedures for estimating q", s hould use future temperature measurements. BeckS was the first "to recognize the importance o f future temperature
information and apply it t o the Duhamel's theorem solution o f the I HCP.
Beck and Wolf' were the first t o c ombine difference methods a nd future temPeratures for the solution o f the I HCP. Beck et al. 6 applied sensitivity coefficient
concepts and substantially reduced the number o f c omputations required for
difference methods with future temperatures. This section considers a single
temperature sensor for the function specification method t hat assumes a
constant heat flux and an arbitrary number o f future time steps.
Suppose the I HCP h as been solved up to time t ",  I; the estimated heat
flux li",  I a nd the entire temperature field
 I a re known. Next. the time is
advanced one step to t ", a nd a n estimate o f q", is calculated. A single temperature sensor is located a t a d epth x . below the active surface. A temporary
a ssumption is m ade that the heat flux is c onstant over r future times. The I HCP
is shown schematically in Figure 6.3. An estimate is s ought o f t he value o f q""
c onstant over r future time steps. t hat minimizes the least squares e rror between
the computed and measured sensor temperatures. The least squares function is
(6.5.20)
T he matrix form o f t he Taylor series expansion o f the temperatures in terms
of parameter space is
(6.5.21)
E quation (6.5.21) is a matrix form o f Eq. (6.5.5). Retaining only the linear terms
in Eq. (6.5.21) a nd substituting into Eq. (6.5.20) gives the normal equations
(6.5.22)
which is a matrix form o fEq. (6.5.6). T he solution can be symbolically written as
,
s= iI [ Y(x., t "'+i d T(x•• t "'+iI. t "'_I' q "'_I' q",)]2
=l
Differentiating S with respect to q",. replacing q", by l i",. a nd setting as/aq",
equal to zero yields,
as
o=~=  2
u q",
(6.5.23)
(6.6.1)
r
I [ Y(x •• t "'+i1)
i =1
 T(x., t "'+iI' t "'_I' q "'_I, l i",)]Z(x., t "'+I_I, t "'_I' q "'_I) (6.6.2)
238
C HAP.6
O NE·DIMENSIONAL INVERSE HEAT C ONDUCTION P ROBLEM
S EC.6.6
239
SINGLE T EMPERATURE S ENSOR
Substituting Eq. (6.6 .5) into Eq. (6.6.2) and solving for qM gives.
q (t)
,
T (t)
L (¥t.M+I  1q M=qt, +1= I
f:'+I' ) Z•. I
(6.6.6)
,
LI (Z.Y
j=
E quation (6.6.6) can be written in terms of a gain coefficient as was Eq. (5.3.4a)
(6.6.7a)
I
i
M IM
M +2
T M +rl
M+ 1
M +r2
Time
11
I
I
where
M IM
M +2
M +rl
M+ 1
M +r2
Time
Z •. i
K •.1 =,
F IGURE 1 .3 Pictorial representation of I HCP using r future temperatures.
L
(6.6.7b)
Z l.j
j= I
where the sensitivity coefficient is defined by
Z
.
'.1 = Z(xto I M +.I. tM I, qM I )_aT(XtotM+il.tMI.qMI.qM)
aqM
(6 .6.3)
Since qM is assumed constant over r future time steps, the Z sensitivity coefficient
is used instead o f X ; see Figure 6.1 for a comparison of the two sensitivity
coefficients. Applying a Taylor series expansion similar t o Eq. (6.3.4) gives
T (Xl. t M+i  l • ' M I • q MI' q M)= T (x • • t M+/ I , t M  I , q MI. q *)
+ (qM  q*)Z(x• • tM + 1 1, t M  I , q M.) (6.6.4)
Figure 6.4 presents a sketch o f the Taylor series expansion. In terms of index
notation, Eq. (6.6.4) c an be written as
T M+i1 ~fM+I1 (A
* )Z
•
.
+ qMq
'.10
.1
1 = . 2. ... . r
(6.6.5)
I f thermal properties are independent of temperature. then the unit step function
solution ljJ used in Chapter 5 is identical to the sensitivity coefficient Z solution of
this chapter and the Z values need t o be calculated only once for a given problem.
If the thermal property variation with temperature is treated in a quasilinear
manner, then Z =1= ljJ a nd Z m ust be recalculated for each value o f qAl' T he subscript k in Eq. (6.6.4) is kept as a reminder that the sensitivity coefficients depend
o n t he depth x . o f the sensor below the active surface; in the next section where
the analysis is extended t o include multiple sensors, the subscript k is even more
important. I t is necessary to remember that the consta{lt heat ftux assumption
indicated in Figures 6.3 and 6.4 is a temporary assumption. The computed heat
ftux qM is retained only over the time interval tM  1 ~ t ~ t M . F or each subsequent
time interval, a new heat ftux is calculated.
The detailed computational aspects of the algorithm given by Eq. (6.6.7) for
the constant heat ftux assumption. single temperature sensor, and r future times
are summarized below:
1. Assume q (t)=q* over tM  1~t~tM+'1 and perform the forward
elimination portion of the TDMA. Starting with the last node N,
T *M+rl
•
1
w N=b '
N
eN=wNcN,
f j=wj(dj+aJj+.),
f N=wNdN
e j=wjcj.
j =NI.N2 . .... 2
f l = W.(dl + atf2). d l =dl(q*)
M IM t M+2 ~rl
M+ 1
M +r2
Time
FIGURE 1 .4
M IM t M+2
I tM+,l
M+ 1
M +r2
Time
Reference condition in the calculation of q ., by differences.
Note that all a j. bj , cj • ej ' a nd Wj are fixed during the calculation of both f 'j'+i  I
Z j . M + 1_ I . i = 0, 1, . .. • r  1 because all thermal properties are held fixed a t
~nd
T'j'.
240
C HAP.6
O NEOIMENSIONAL I NVERSE H EAT C ONDUCTION P ROBLEM
OM
°
'
.
..
2. C aIcu Iate T j ' TOM+I , . .• , T M+, _ I usmg the back substitutIOn portion
j
j
o ftheTDMA,
o
T M +iI_ f I.i
I

M+
t M+IIe) t ) 1 1 1+ J'),i'
I
)

j =2,3, . .. , N}
Repeat for
; =O,I, . . . , rl
T he ; subscript on h,i is a reminder that this term must be recalculated for
each value o f future time (;). However, the e/s are identical to those o f step 1
and do not change for ; = I, 2, . .. , r.
3. Calculate the step function sensitivity coefficients ( Z) using the difference
equations developed from Eq. (6.2.7). The TDMA coefficients a j' bj , Cj ' e j' and
W j have already been calculated in step 1. Therefore only dj a nd j j need be
recalculated and the back substitution must be performed.
Z I,M+II
} Repeat for
= f t.i
Z j,M+il=ejZjI,M+il+h,i>
j =2,3, . .. , N
S EC.6.7
6.7 M ULTIPLE T EMPERATURE S ENSORS, F UNCTION
S PECIFICATION ( q=C), , F UTURE T IME STEPS
Multiple temperature sensors were considered in Section 6.4 for a single future
time step. In this section, the analysis is extended to include an arbitrary number
of future time steps with the temporary assumption of constant heat flux. T he
least squares error function must be modified to include a summation over the
number of sensors. F or J sensors and r future times, define a sum of squares
function by
J
,
t= I
i= I
LL
S=
[ y",M+jIT(xt,t M+i  I ,IMI,QMI,qM)]2
oS
N ote that ! he dj a nd h i terms for the calculation of Z j,M+i  1 are different from
those for T f + i  I. Also, the boundary conditions o n Z are different from those
on
4. Calculate qM from (6.6.7a)
t.
,
(6.7.1)
The value of qM, constant over r future time steps, that minimizes S is sought,
J,
LL
, ,=0=2
uqM
t =1
; =O,I, . .. , r1
241
M ULTIPLE T EMPERATURE S ENSORS
[ T(Xl>tM+jl,tMI,qMI,qM)y",M+jI]Zt,j
(6.7.2)
j =1
where the step function sensitivity coefficient is defined by
 Z(
Zt,j
Xl>
t M+j  l , t MI'
)
q MI
o T(Xt,tM+jl,tMI,qMI,qM)
OqM
( 673
.. )
Expanding the temperature field in a Taylor series about an assumed heat flux
"
* '"
q M=qM+ I... (Y.t ,M+II TOM+iI)K t ,i
t
q*,
i= I
5. Knowing qM, reevaluate d l a nd I I a nd repeat the back substitution
portion o f T DMA for t (x, tM, tM  I ' q M  I , qM)
T :+jl=t:+JI+(qMq*)Zt,j,
j =I,2, . .. , r1
(6,7 .4)
substituting Eq. (6.7.4) into Eq. (6.7.2) a nd solving for the desired heat flux,
t~=fl
t f = e)Tf1 + f i'
j =2, 3, . . . , N
(6.7.5.)
EXAMPLE 6.2. Using the information in Table 1.1, calculate the K factors o f Eq. (6.6.7)
for r =2,&I+ = (J.!u/L2 =0.05, a nd a temperature sensorlocated a t x+ = x/L=0.5 . Repeat
the calculation for x+ = x/L= 1.0.
Equation (6.7.5) can be written more compactly as
J
So/ulion.
&1+ = 0.05, x+ = 0 .5
Z . = 0 .0153659
Z l =0.0593109
0.0153659
K  ~1  (0.01 5 365W + (0.0593109)2
= 4.0933
K
0.0593109
2 (0.01 5 365W +(0.0593109)2
=15.80
&1+ =0.05,
q M=q*+
x+ = 1.0
,
LL
jl
Kt,j(Y,..M+j  lt:+  )
(6.7.6a)
t= I j= I
Z . =0.0002693
Z2 =0.0078853
where
K_
0 .0002693
•  (0.0002693)2 + (0.0078853)2
(6.7.6b)
=4.3261
0.0078853
K  ~~~
2  (0 .0002693)2 + (0.00785W
EXAMPLE 6.3. Repeat Example 6.2 with the modification that two sensors ( J = 2)
a re located at x+ = 0.5 a nd x+ = 1.0. F rom Eq. (6 .7.6b),
= 126.7
The x+ = 1.0 a nd 1+ = 0.05 case is also included in Table 4.1.
o
Z •. i
K •. i = Z2
Z2
Z2
Z2
•.• + 2.• + •. 2 + 2.2
b Zt,j
SEC. 6 .9
CHAP. 6
242
Note that it is necessary to use the pulse sensitivity coefficient X .,I_ j + 1 =
because the heat flux is allowed to vary in steps over the r future
time steps. Substituting Eq. (6.8.1) i nto Eq. (6.8 .3) gives
Solution. F rom E xample 6.2.
A Z.,I _ j + I
Z I.I = 0.01537
Z2.1 = 0 .0002693
ZI.2=0.05931
b =262.0
Z 2.2 = 0.007885
K I.I = 4 .027
K 1 • 2 =15 .54
K 2. 1= 0.07057
T :'+ j I = t :,+ j I + q",I/I.,j qMI I /I.,jl  q*Z.,i
K 2 • 2 =2.066
L
i
i
1
o
T he presence o f a second sensor a t x + = 1 does not add significant information
as evidenced by the fact t hat K 1 .1 a nd K I . 2 for this problem are very close to
KI a nd K2 for Example 6.2 (for x + =0.5). T his is because o f t he relative size
o f the sensitivity coefficients for x + =0.5 a nd 1.0.
A Z• . I _j+I=Z•. i >
L
j =1
j AZ•. I _ j + l
=L
Z .,j=I/I.,j
(6.8.5)
j =1
Substituting Eq. (6.8.4) i nto Eq. (6.8.2) yields
s=
t
( t:'+II Yt."'+I1 + q",I/I.,IqM_II/I•. I _Iq*Z.,1)2
(6.8.6)
i =1
The heat flux component q", _ I is treated as a quantity known from the previous
time step. T he value o f qM t hat minimizes S is found t o be
Repeat Example 6.3 for t wo identical sensors. b oth a t x+ = 0.5.
Solution.
(6.8.4)
where the following relationships are used:
j= I
EXAMPLE 6 .4.
243
SECONDORDER SEQUENTIAL REGULARIZATION M ETHODS
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
b = 133 .2
K 1 • 1 = K 2 • 1 = 2.047
o
K 1 • 2 =K 2 • 2 =7.900
6 .8 SINGLE T EMPERATURE SENSOR. FUNCTION
SPECIFICATION, L INEAR HEAT FLUX
(CONNECTED S EGMENTS)
(6.8.1)
U p t o this point. all the function specification methods in Chapter 6 have utilized
the (temporary) assumption o f a c onstant heat flux. O ther assumed functional
forms were considered in Chapters 4 a nd 5. All o fthe h eat flux variations utilized
with the convolution integral methods in Chapter 5 c an also be utilized with the
difference methods o f this chapter. An example is the connected linear heat flux
functional form shown in Fig. 4.8. T his approach is based on a straight line
e~tra~olation in time through the two points q"'_1 a nd q",. F or equally spaced
time mcrements. the heat flux variation can be expressed as
q "'+II=(li)qMI+iq",.
i =O.I • ...• r
(6.8.1)
N ote t hat Eq. (6.8.1) c ontains only a single unknown. qM' Additional details can
be found in Section 4.4.3.2. T he following least squares e rror is minimized
,
s=
L
[ Y(x •• t M+1_I)T(x•• t "'+I_I. t "'I' q MI. q", •. ..• qM+,_1)]2 (6.8.2)
1 =1
A comparison of Eqs. (6.8.2) a nd (6.6.1) indicates t hat t he temperature field for
the linear heat flux case depends on the r heat flux values q", • .. . • q M+,I;
however. there is only one independent heat flux (qM) since they are related
through Eq. (6.8.1). Using the standard Taylor series expansion. the temperature
field can be written as
1
T M+1IT"M+II
•
.
+ " AZ. ,Ij+1 (q M+jIq *)
~a
i I
(6.8 .3)
A c omparison o f the constant heat flux result. Eq. (6.6.1). s hould be m ade with
the linear heat flux result Eq. (6.8.1). Both are o f t he same general form b ut
the coefficients are different a nd the latter depends explicitly on q",  I ' T he
linear heat flux variation yields a more complicated equation t han does the
constant heat flux result. Experience has shown that the linear heat flux
algorithm requires larger time steps for stability t han does the constant heat
flux a lgorithm for the same value o f r.
6.9 S ECONDORDER SEQUENTIAL REGULARIZATION
M ETHODS
The general aspects o f regularization methods were discussed in Section 4.5
a nd the convolution integral method was applied to firstorder sequential
regularization in Section 5.4.3. In this section. the equations for the secondorder sequential regularization will be developed for the case o f r = 3 a nd then
generalized by using matrix notation for an arbitrary r. T he matrix derivation
also applies t o o ther orders o f regularization.
The secondorder regularization simulates the continuous term
r'"''
j,.,1
d2
I (
i)2 dt
(6.9.1)
dt
This integral is used in the cubic interpolatory spline 7 a nd the minimization
o f it yields the socalled minimum curvature property.
244
CHAP. 6
O NEDIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM
Following the material in Section 4.5, the standard least squares function
is augmented as follows:
SEC . 6.9
2 45
SECONDORDER SEQUENTIAL REGULARIZATION M ETHODS
3
L
[ f:'+II Y"'+iI
1 =3
r 2
,
S=
L
( Y"'+II T:'+II)2+CX
1 =1
L
( q"'+I_I2q"'+I+q"'+I+d 2 (6.9.2)
1=1
i
+
L
X 1,Ij+l(q",+jlq·)]X1,i2+ CX (q",2q"'+1 + q"'+2)=O
(6.9.6c)
J= 1
with r~ 3. N ote that cx is dimensional and has units of T/q2. I n order to keep
the development simple, a value of r =3 is assumed; therefore, the second
summation in Eq. (6.9.2) contains only one term. Differentiating Eq. (6.9.2)
with respect t o the three unknown heat flux components q"" q", + 1 , a nd q", + 2
a nd setting the results equal t o zero gives,
After straightforward but lengthy algebra, Eq. (6.9.6) can be written as
C tl X t'iXl'i+ CX )
C tl X1,,+1 X t,i 2CX )
C tl Xt ,i+2 Xl,i+ CX)
C tl X 1,i X l,i+I2CX)
C tl X 1"Xt ,i+ 4CX )
C tl X t ,i+ I X 1,,2CX)
(.±
(.±
(.±
1 =1
X l,iXl ,i+2+CX)
p =l
X1,iXl,i+I2CX)
. .= )
X1"X1'i+ CX)
3
L
X 1,;(Y"'+I_If:'+il)
1= 1
2
(6.9.3c)
L
• M +·
X1,;(Y"'+i  T t
(6.9.7)
')
i =1
1
T he lower limit o n the summation does not always start with unity because a
sensitivity coefficient is identically zero for all times prior to the initiation of the
particular heat flux component. Sensitivity coefficients are defined through
i JT:'+I1
X1•1 j+ 1
iJq"'+j_1
(6.9.4)
The usual Taylor series expansion is written as
1
T I"'+II = T IO"'+II+'" X t .lj+1 (q "'+jIq . )
L...
(6.9.5)
j= 1
Substituting Eqs. (6.9.4) a nd (6.9.5) into Eq. (6.9.3), one obtains
3
' " [TO " '+/1
L...
t

iI
Y"'+/1
.
L
X t ,lj+l(q"'+jIq·)]X1 ,I+CX(q",2q"'+1 + q"'+2)=O
i =1
The symmetrical coefficient matrix affords some computational savings. The
regularization parameter cx improves the diagonal dominance o f the coefficient
matrix and reduces the sensitivity to measurement errors. Since the equations
are linear in q jq., the value o f q . is completely arbitrary. If q ·=O is used,
then f is replaced by Tlq=o. I n the sequential procedure, generally only q",
found from Eq. (6.9.7) is retained. The procedure is then repeated for the next
time step.
F or the secondorder regularization, the smallest allowable value o f future
times is r = 3. Even for this simplest case, the algebra is rather lengthy. In order
to extend this procedure to arbitrary r, matrix methods are used. The matrix
analog of Eq. (6.9.2) is
S =(Y  T)T(y  T)+cx(H 2 Q)T(H 2 q)
1
+
'L" Xt,i (Y." '+/+1 TO"'+i+l)
...
t
(6.9.6a)
j= 1
The Y and T matrices have already been defined in Eq. (6.5.13). Since the
procedure being developed is sequential, the heat flux vector contains only r
elements,
3
' " [TO I" '+/1
L...
1 =2

Y"'+/1
.
Q=
L1 X 1,IJ+I(q"'+jIq·)]X1,iI2cx(q",2q"'+1 + q"'+2)=O
j=

q",
q"'+1
:
[ q "'+r2
1
+
(6.9.6b)
(6.9.8)
q "'+r\
]
(6.9.9)
246
C HAP.6
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
H 2=
0
0
0
0
0
1
2
1
0
00
10
2
0
0
0
0
1
r 2
0
0
0
0
2 1
00
00
01
00
00
.
T =T+X(qq*)
(6.9.10)
1
(6.9.20)
(6.9 .11)
for T, and q* is a vector of r elements all equal to q*, and X is the sensitivity
coefficient matrix,
Xl,r
X l,l
0
0
0
X1,r  1
0
0
0
(6.9.12)
0
X1. 2 Xl,.
Note that X is lower triangular and each element is defined by Eq. (6.9.4) . The
matrix derivative of S with respect to the heat flux vector q is given by
VqS = 2[Vq(Y  Tn(y  T) + 2exVq(H 2q)T(H 2q) ==0
Equation (6.9.19) is a matrix form o fEq. (6.9.7) a nd similarities can be observed.
t is replaced by Tlq=o. With matrix manipulation subroutines
available, Eq. (6.9.19) is much preferred over a generalization of Eq. (6.9 .7) t o
arbitrary r. The solution to Eq. (6 ,9.19) can be written symbolically as
I f q*=O, then
.
where the elements of T are defined in a m anner completely analogous to that
0
0
247
(6.9.19)
F or r = 3, only the first row o f H2 is nonzero. The matrix equivalent of the
Taylor series expansion, Eq. (6.8.3), is
Xu 0
Xl.2 Xl,l
X = X u Xl ,2
SPACE M ARCHING TECHNIQUES
The matrix product H21 is the vector, the elements of which are the sum of
the elements of each corresponding row o f H 2' T he elements of H 21 indiv,idually
sum to zero, Equation (6.9.16) can be written as the matrix normal equation
By analogy with Eq. (4.5.16c), the H2 matrix is
2
SEC. 6.10
(6.9 . 13)
F rom Eq. (6 .9.11),
I t is c ommon in the sequential procedure t o retain only the qM c omponent o f
q and repeat the process for each successive time step. In the linear version o f
the sequential regularization, the thermal properties are maintained constant
at values corresponding to the temperature field at time t M _ I ' I f t he thermal
properties are independent of temperature over the range o f interest, the whole
domain regularization equations are identical to those o f Eq. (6.9.20) provided r
is the total number of heat flux components being estimated. Since the sequential
regularization method requires the solution of a system of r algebraic equations
for each value of qM, it is computationally more costly than the other sequential
methods presented in this chapter but for small values of r such as r~ 5 the
increase in computation may not be large.
The individual elements of X are calculated by means similar t o those
discussed in Sections 6.2 a nd 6 .3. Namely, X 1,jl+ 1 is determined from the
difference solution of a direct problem in which the initial profile is zero, the
heat flux is unity during the time interval tM+i  2~t~tM+il a nd zero thereafter, and the inactive boundary is perfectly insulated.
F or a discussion of the choice o fthe regularization parameter ex, see Chapter 5.
Although Eq. (6.9.19) is derived for a secondorder regularization, it also
applies to other orders. F or a zerothorder regularization, for example, H2
is replaced by Ho which is developed in Section 4.5 .
(6.9. 14)
Also,
Vq(H 2q)T = vq(qTHIl=HI
(6.9 . 15)
Substituting Eqs. (6 .9.11), (6.9 . 14), and (6.9.15) into Eq . (6.9 . 13) and simplifying:
XTX(q  q*)+exHIH 2q =X1"(y 
t)
(6 .9.16)
Equation (6.9.16) can be put in a different form by noting that
exH~ ' H2q*=0
(6.9.17)
I f I is a vector with all elements equal unity, then for q* = q*1,
H 2q*=q*H 21 = 0
(6.9. 18)
6.10 SPACE M ARCHING T ECHNIQUES FOR
O NEDIMENSIONAL P ROBLEMS
T he difference methods presented in the previous sections can be termed "time
marching" because starting with the initial temperature profile, the entire
temperature field is calculated at the next time step. O ne marches ahead in time.
F or the inverse problem, two boundary conditions are known at the sensor
location; namely, temperature, and heat flux. The knowledge of two conditions
at the same spatial location suggests that the I HCP is related to an initial value
problem instead of a boundary value problem. This section explores several
space marching techniques which are applicable to "initial" value problems.
248
CHAP. 6
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
I n the analysis that follows, a planar geometry with constant thermal properties
is assumed; these restrictions are for simplicity only.
F or the IHCP, the geometry can be divided into inverse and direct regions,
as shown in Figure 1.4. T he direct region is the classical boundary value
problem; the experimental temperature measurements are the boundary conditions for the interface between the direct and inverse regions. The boundary
conditions are known a t the remaining boundary. Within the direct region,
it is possible to solve this wellposed problem. Once the temperature field is
known for the direct region, the heat flux a t the sensor location can be determined by differentiating the temperature profile. In the discussions that follow,
all temperatures and heat fluxes are assumed known within the direct region;
hence, both temperature and heat flux are known at the inverse boundary.
6.10.1
A nalytical S olution
SEC.6.10
Equation (6.10.3) demonstrates the significance of the dimensionless time step
for inverse problems.
.'
I f an error is made in the measurement of Y , there wIll be errors m both
I
I
q£,' and qO , I' A lower bound o f this error can be . estimate~ by re~lacing Y by
Y, + r, YI a nd ignoring the error in qli,l' The resultmg error m qO,1 IS found t o be
t:.t£=~6t/E2
r,QI = _1_ sinh ( _1_)
b YI
(6.10.1)
I f the time derivative is replaced by a backwards difference, Eq. (6.10.1) can be
converted from a partial differential equation t o a n ordinary differential
equation
(6.10.2)
where the superscript i denotes time index. The boundary conditions for Eq.
(6.10.2) a re
I
x =1i
Equation (6.10.2) has the appearance of an initial value problem because the
function and its first derivative are both specified at x = E. I fthe initial temperature profile (To) is uniform throughout the body, a very simple solution of
Eq. (6.10.2) exists for i = l (t = 6t)
T '(z) To=(YI

Yo) cosh (
~)+qli",J;Xi sinh ( _ z _ )
~~6t
, J;Xi
where z =Ex. The unknown flux a t x =O, ( z=E), is found by differentiating
this equation with respect to x t o o btain
O,1 E qli" E
QI =qk=k cosh ( E ) + (Y,  Yo)
C A.
~~6t
~
+q~1 ~sinh(~)J
An analytical solution of this equation is possible; however, due to the complexity o f the inhomogeneous terms, the algebra is more c0'!lplicated an~
solutions a t subsequent times even more so. Consequently, this approach IS
not pursued further. The important points t o remember from this analysis are
that for a given value of time, it is possible to march in space, and that the
important length scale in the dimensionless time step 6 tE is the sensor distance
below the heated surface, E.
6.10.2
M ethod o f D 'Souza
D'Souza 8 ,9 used the pure implicit difference approximation t o Eq. (6.10.~);
the time derivative was replaced by a backward difference a nd the spatIal
T i(x=E)= Y;
q£.i=  k dTi(x)
dx
~
which is tabulated in Table 6.6; these results are consistent with earlier observations that small values of dimensionless time 6 tli can produce large errors in
surface flux.
The foregoing process can be repeated for t = 2 6t by solving
2
d T2(Z) _ _1_ T2(Z) = _ _1_ [ To+(Y,  Yo) cosh ( ~)
dz 2
~6t
~6t
~~6t
Consider the constant property heat conduction equation for a planar geometry
02T(x, t) 1 o T(x, t)
iJx2 =~  ot
2 49
SPACE M ARCHING TECHNIQUES
E
C A.
~~6t
.
smh ( E ) (6.10.3)
C A.
~~6t
T ABLE 6 .6 D imensionless
H eat F lux E rror f or
S ingle T emperature
E rror «SY,
d "tE
cSQd{)Y,
10
1
0.5
0.25
0.125
0.0625
0.03125
0.1017
1.175
2.737
7.254
23.843
109.160
809.618
250
C HAP. 6
O NEDIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM
S EC.6 .10
261
S PACE MARCHING TECHNIQUES
(6.10.4)
T he heat flux qh is unknown a nd is t o be estimated while q~ is the heat flux
leaving t he inverse region shown in Figure 1.4; this latter heat flux is e ntering
the direct region a nd is found from calculations for the direct region. I n Eq.
(6.1O.6b) t he temperatures T~' a nd T~ a re replaced by l j_1 a nd l j, respectively.
E quation (6.1O.6b) is solved for T~  I for i = 1, 2, . . . , n t o g et
(6.10.5)
(6.10.7)
where p =aAt/Ax 2 • U sing t he spacetime grid shown in Figure 6.5, t he calculation procedure is s tarted a t n ode j = K a nd sequentially steps t hrough time i = 1,
2, . .. , n. Because t he t emperature field is known along b oth t he s patial grid lines
K a nd K + 1, Eq. (6.10.5) is a n explicit relationship even t hough it was derived
from t he fully implicit equations for direct problems. After t he t emperature is
c alculated all along space line K  1, t he procedure is repeated for K  2, . .. , 1.
A finite difference e quation for node j = 1 yields a n e quation for calculating
t he 'surface heat flux.
T he D 'Souza p rocedure is intended t o b e implemented o n a c omputer as
just outlined. Insight c an b e gained, however, by obtaining some algebraic
expressions for a small n umber o f nodes, K , from the heated surface t o the
sensor. U sing t he p ure implicit procedure, equations a re given for each node.
F or n ode 1, t he h eated surface, the e quation a t time t i is
T i _ Ti
qi + k 2
I
(6.10.6a)
o
Ax
F or K > 2, Eq. (6.10.7) is i ntroduced into E q. (6.10.4) w ritten for j = K  1 t o
obtain a n expression for T~2' i = 1, 2, . . .. F or K = 2, Eq. (6.10.7) gives expressions for T~ a nd (by replacing i by i I) T~' which a re i ntroduced into Eq.
(6. 1O.6a) t o o btain
d erivative by a central difference
T~_,2T}+ T}+,
Ax 2
Solving for T~_, in Eq. (6.10.4) yields
i
T J I = (2 +
!)
T i._TI
J
p
J+ I
_!
P
Ti.  I { j=K.Kl • ...• 2
J
.
1 = 1• 2, ... , n
T he e quation for a n i nterior n ode j is given by Eq. (6.10.4) a nd t he e quation
for n ode K is
(6.10.6b)
(6.10.8a)
i _q~E
QE =T
(6.10.8b)
since Ax = E for this case.
T he s ymbol V is the backward difference o perator,
V lj= l jlj_I'
V llj= l j2lj_, + l j2
(6.10.8e)
I f t he foregoing procedure is used for K =3, a nd t hus A x=E/2, t he heat
flux expression is
i
V lj
3 V 2lj
1
V3lj
i
1 VQ~
1 V2Q~
Qo= AtE + 16 ( A't'd+ 128 (A't'E)3+QE+ 2 A't'E + 32 (V't'£l2
0
Direct region
Inverse region
0 1
•
D"Souzal8.91
molecule
•
r 
~
OJ
•
JI~" OJ)"
 "\
0 2
"
II
E
i=
.)
•
( i.
I
i 2
2
•
j l
r
I
I
I
\ !I
j
•
j +l
Known
conditions
•
•
I nitial conditions
j =O
1
•
•
2
3
X 2
X I
Space index
x
X +l
FIGURE 6 .6 Spacetime grid for space marching methods.
(6.10.8c,d)
(6.10.9)
Several c omments regarding Eqs. (6.10.8a) a nd (6.10.9) a re n ow made. First,
this is a procedure t hat utilizes exact matching o f t he calculated temperature
with the measured temperature a nd t hus is sensitive t o m easurement e rrors.
Second, these equations a re n umerical analogs o f t he e xact solution which is
given by t he B urggraf equation, Eq. (2.5.20). T he coefficients o f Eqs. (6.10.8a)
a nd (6.10.9) a pproach t hose o f Eq. (2.5.20) as the n umber o f n odes increases;
the first difference coefficients in b oth o f t he former equations a re the same as
each o ther a nd t hose in Eq. (2.5.20) whereas the V2lj s tarts a t 0.25 [given by
Eq. (6.10.8a)], decreases t o 0.1875 in Eq. (6.10.9), a nd a dditional nodes bring
the coefficient closer t o t he exact value o f 0.1 667. Third. only past temperatures
are used t o o btain t he p resent heat flux estimates. I t h as repeatedly been noted
that future temperatures a re needed t o o btain accurate estimates with small
time steps. F ourth, if a n a nalogous procedure were utilized starting with the
explicit equation [i.e., left side o f Eq. (6.10.4) e valuated a t i I], forward difference o perators w ould replace the backward difference operators o f Eqs. (6.10.8a)
a nd (6.10.9), e nsuring t hat o nly future temperatures a re used. Fifth, as m ore
I
r
262
CHAP. 6
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
SEC.6.10
nodes are used, higherorder differences are needed, increasing sensitivity to
temperature measurement errors.
6 .10.3
253
SPACE M ARCHING TECHNIQUES
:0~
n = 1 0·
I
9·
)(
Weber worked with the hyperbolic form o f the constant property heat conduction equation
molecule (R&Blll,12
7· )
6•
X
0
0
X
X
0
)(
0
X
)(
X
0
0
X
X
~
3
' ;'I
I
X
\
I
"5.I
~ _ __
~
X
r  _J
X
1
L(TiI_2Ti+I+Ti+I)+_(Ti.+I_TiI)=~(Ti _ 2Ti+Ti. )
.&t2 j
j
j
2.&t J
j
.&X2 j + I
J
rI
)(
X
)(
~ __ .0
I
2
X
(6.10.11a)
i =00
j= 1
Solving for T }_I, Eq. (6.10.11 a) yields,
Cl.&t
( 1=
y 2 (.&X)2
IX
.&t
K Future times
)(
where y is a normalized thermal wave speed. Central differences in both time
a nd space were used,
p = .&x 2 '
0
X
_~'1
5 •_
II
i=
(6.10.1lb)
0
X
)(
E n K+l=4(x _ __ X I
(6.10.10)
0
0
x
lO
I
o
Computation
8·
M ethod o f W eber
r 
!
r0='~
x
oo
3
2
S pace,j
0
Computation
'..£;: molecule (Weber10)
, __~
,
~
0
oooo
7
o
X
4
x
5
)(
6
K l
K
8
K +l
FIGURE 6.6 Spacetime grid for Weber'o and Raynaud and B ransier"·
(node points), n = 10 (temperature measurements), 0 known. x calculated.
11
(6.10.11c,d)
This algorithm is explicit; it can be started a t node j = K a nd time i = 2, 3, . .. ,
n  1. T he value i = n is not permissible because T~+ I is n ot defined. Consequently the maximum time index for the grid locationj is one less than that for the grid
location j + 1. This grid pattern is shown in Figure 6.6. I f the surface temperature
is desired a t time M, then the temperature must be measured out to time
i =M+K. T he algorithm given by Eq. (6.10.11b) has the appearance o f using
K future times where K is the number of spatial grid points. Refining the grid
has the effect of increasing the number offuture times. For the methods presented
in Sections 6.3 a nd 6.4, the number o f future times is dependent o n the physics
o f the problem and should be relatively independent o f the number o f spatial
grid points.
Weber indicated that the parameter (1 should be chosen small, without
indicating how small. In fact, the algorithm works with (1=0.
Algebraic expressions similar to those for the D'Souza procedure can also
be found for Weber's. Following the same procedure as in Section 6.10.2 one can
derive for two nodes; that is, K = 2, a nd with ( 1=0
(6.10.12b)
where Q~ is a set equal t o zero. These equations are quite similar to Eqs. (6.10.8a)
and (6.10.9), found for D'Souza's method, and comments one, two, and five
made in connection with the latter equations also apply for Eq. (6.10.12).
However, Eq. (6.10.12) has two advantages over Eqs. (6.10.8a) .and (6.10.9):
(I) both past and future temperatures are used and (2) central differences are
used that are much less sensitive t o measurement errors.
6.10.4
M ethod o f R aynaud a nd B ransier
The method o f Raynaud and Bransier ll ,I2 consists of a twoste~ marching
procedure. The first consists of a space marching step followed by a tIme marching step; the final result is an average of the two procedures. F or the space
march, the heat conduction equation is differenced as follows:
.&x  .
.
k . i
k
p c(T,.+I_Tj)=(TjTjd+ "A (Ti+t l  f i+l)
j+
j
(6.10.12a)
...
methods; K =7
.&t
J
.&x
(explicit)
aX
(implicit)
( 61013)
..
264
CHAP. 6
O NEDIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM
' Note t hat the conduction terms are mixed explicit/implicit. Solving for t i. I
from Eq. (6.10.13),
J
t~_1 = (1 !) t i+(1 +P t J.+ t~+
!) i J+
P
J
1:...
I
I
the same as the Weber algorithm; see Figure 6.6 for details.
The second part o f the algorithm approximates the heat conduction equation
as follows:
(implicit)
(explicit)
Again, the conduction terms are mixed explicit/implicit. Solving for
from Eq. (6.10.15),
...
T~:!:~=(I+~) T~+I+(I~) T~T~+I
of this approach, see References 13 a nd 14. Some o f the basic features of this
method can be understood by working with the combined form of Eq. (6.10.18).
t he socalled leapfrog difference scheme was used,
.1x
.+ I
.I
e x.
..
 (T'.  T'  )=(T"+12T'+T"_d
j
2M J
. 1xJ
J
J
TI+
I
j I
(6.10.16)
.
6.10.5
oq
ox
(6.10.18a)
oT
ox
(6.10.20)
~ ~
)(
9x
)(
x
)(
x
)(
8x
x
)(
)(
x
x
G>
G>
7 )(
X
)(
X
X
x
G>
G>
6 )(
)(
X
X
X
)(
G>
G>
5 )(
)(
)(
X
)(
x
I
.,
(6.10.21)
0 __ /i>, _~
_
)(
I
\@ G>
'C!)I
X
X
X
x
3 )(
x
x
x
x
)(
X
X
X
I
I
I
~,
r?
x
X
X
)(
G>
G>
0
G>
I
X
'@
X
~ G>
,, __ J
1 )(
G>
,
rx' G> '  0'
I
i=OG>
0
I
I
" ' _ _ "
~§._%
2
3
4
5
(6.10.18b)
G>
G>
0
K l
j= 1
and then converted them to finite difference equations individually; the two
unknowns were T(x, t) a nd q(x, t). F or a discussion of the purported advantages
.
 Tj+1
)(
and
q =k
' +1
)(
2 )(
equations
oT
ot
1
.
4 )(
Hens~l a nd Hills 1 3 wrote the heat conduction equation as a pair o f first order
p c= 
I
n = 1 0 )(
E
i=
M ethod o f H ills a nd H ensel
1.
2p T r + 2Tj+ 2p T j
.
I
ex
_1
1x (T.2  Tj)=(T jI + l  2T jI + TjI )
M
J
.1x
(6.10.17)
Although there is a time marching component to the Raynaud and Bransier
algorithm, the com~ined e~ect ~fth~ two steps has the same overall appearance
as the Weber algorithm given 1 0 Figure 6.6 but it is even less sensitive to the
~eas~rement e rrors than Weber's method. The apparent number of future
tImes IS the same as the number o f spatial grid points ( K).
Raynaud and Bransier also considered temperaturedependent properties.
=
See Figure 6.7 for the computational molecule for Eq. (6.10.20). (A "computational molecule" is a schematic representation o f a difference equation such as
Eq. (6.10.20).) T he preceding algorithm is explicit and is identical to the Weber
algorithm, Eq. (6.10.11 b), for u = 0. Although Eq. (6.10.20) is valid for j = 1,
Hensel and Hills l3 chose to replace the centered time difference in Eq. (6.10.19)
with a forward time difference
T he Ti:!: ~ al.gorithm is explicit and time marching. Starting with j = 1, the
~pace mdex} steps through K  1, K  2, . .. , 1. T he final step o f the algorithm
IS a n average o f t he foregoing two procedures.
T 'j =Z(tlj +Tj)
~
•.
(6.10.19)
Although Richtmyer and Morton lS indicate that the foregoing difference
technique is unconditionally unstable for direct problems, it can still be useful
for inverse problems. I f Eq. (6.10.19) is solved for T~_I'
Tj_1
(6.10.15)
266
SPACE M ARCHING TECHNIQUES
(6.10.14)
~he e~plicit algorithm for t~_1 s tarts with n odej = K a nd steps through the time
hnes 1 =2, 3, . .. , n  1. h i general, the maximum value o f j is n +j  K which is
. 1x...
•.
k.
•
k ..
pc  (T'.+ 1 _ T j)= (Ti.+ 1 _ T i+ 1 )+_ ( T i
T i)
.1t J
.1x J
j l.1 x j +l j
SEC.6.10
K
K +l
6
7
8
Space
F IGURE6.7 Spacetime grid for Hills and Henscl'3.'4 method; K =7 (node points), n =10
(temperature measurementtimcs), 0 's known, x 's to be calculated.
2 68
C HAP.6
O NEDIMENSIONAL I NVERSE HEAT C ONDUCTION P ROBLEM
Solving for 1'/_ 1 ,
1
1 ')1 = ( 2 
1)
P
1
1
1')  1')+ 1
1)
+p 1'2
(6.10.22)
The computational molecule for Eq. (6.10.22) is shown in Figure 6.7 .
Equation (6.10.20) is n ot applicable for i = n. Instead of using the approach
of Weber and generating the solution up to the diagonal line of Figure 6.6,
Hensel and Hills l3 replaced the time derivative with a backward difference
(6.10.23)
Solving for 1 'j1o
1 ' 1=
)
(2 +p1) ) p)
1
1 '
1 '1
 1'"
)+1
(6.10.24)
which is the same equation used by D'Souza. The computational molecule for
Eq. (6.10.24) is shown in Figure 6.7.
The depth of the discussion presented here does not do justice to the analysis
presented by Hensel and Hills l3 and Hills and Hensel. I4 They extended the
analysis to consider temperaturedependent thermal properties, nonuniform
spaceandtime grid, planar, cylindrical, and spherical geometries, and an
estimate ofthe variance of the calculated surface heat ftux provided the variances
of the temperature measurements and heat ftux a t the inverse boundary are
known. Also, their algorithm has the capability of considering a digital filter
of the measured temperatures.
6.10.6
C omparison w ith P rior M ethods
A major weakness of most of the methods in Section 6.10 is the lack of a statistical
basis. Another is that the methods do not readily extend to multiple interior
sensors and to two and threedimensional cases. The function specification
and regularization methods can have a statistical basis and the same formalism
applied to onedimensional problems can be directly extended to two and
threedimensional cases. Moreover, quite different problems such as those
involving integral equations o r ablation can be treated by the latter methods
whereas the space marching techniques m ayor may not be appropriate.
S EC.6.11
I t is appropriate to ask what difference method is the most accurate in the
calculation of sensitivity coefficients. Although an exhaustive study has not
been pedormed, Figure 6.8 provides some insight into this question; this figure
considers a sensor at the midplane of a planar slab. The sensitivity coefficient
is evaluated at a dimensionless time equal to the dimensionless time step of the
difference method; that is, t + = At+. The backward difference ( 8= 1) method is
superior for At+ near unity and larger whereas the central difference method is
superior for At + < 0.1.
As pointed out earlier, the sensitivity coefficient calculation time step A tx
can be smaller than the problem time step At. This should produce more
accurate sensitivity coefficients.
Chapter 5 presented several example calculations related to a triangular
heat ftux pulse. This same problem will be utilized here for a number of examples. The triangular heat ftuX pulse is maximum at t+ = 0.6 and then falls
linearly to q+ = 0 a t t + = 1.2. The analytical solution of this problem was
presented in Chapter 5 [see Eqs. (5.2.4)  (5.2.6). Table 5.1 tabulates the exact
dimensionless temperature response (to six decimal places) for a sensor located
1 .o1.__
0.1
0.01
o
6.11
2 67
N UMERICAL CALCULATIONS
N UMERICAL C ALCULATIONS
 Analytical
0 6 = I , Lumped capacitance
o 6 = 112, Lumped capacitance
[ )'I/[)'Ix
This section presents numerical calculations that bring out the salient features
of several difference methods presented in this chapter.
Most methods in this chapter require the calculation of various heat ftux
sensitivity coefficients (iJ1'jiJq). The most powerful of these use differences.
=1
o
0.OO~:::
. 0c:1~;:0";. 1~~t1.0
[ )'I+
= a [).tlL2
F IGURE 6 .8 Calculation o f p lanar slab sensitivity coefficient by difference methods for x /elL =0.5.
268
C HAP. 8
O NE  DIMENSIONAL INVERSE HEAT C ONDUCTION PROBLEM
T ABLE 6 .7 T emperature
R espon.e t o a T riangular
H eat F lux P ulse. A dditive
N ormal U ncorrelated
R andom E rror. w ith
S tandard D eviation o f
0 .0017
 0.24
 0.18
 0.12
 0.06
0.0
0.0600
0.1200
0.1800
0.2400
0.3000
0.3600
0.4200
0.4800
0.5400
0.6000
0.6600
0.7200
0.7800
0.8400
0.9000
0.9600
1.0200
1.0800
1.1400
1.2000
1.2600
1.3200
1.3800
1.4400
1.5000
0.00034
0.00281
0.00135
0.00090
 0.00020
 0.004241
0.000639
0.001307
0.007161
0.012874
0.021764
0.038954
0.054821
0.072952
0.098381
0.127722
0.155529
0.193275
0.223066
0.248541
0 .275499
0.293646
0.311805
0.330811
0.342215
0.348997
0.350227
0.355722
0.359318
0.361609
a t the insulated inactive surface Te
by using a random numbe
. mpe~ature measurement errors simulated
deviation o f 11+  00017 r ge~erat~r With normal distribution and standard
r .
are given I D Table 6 7Th '
approximately 5% o f the m '
d'
...
IS standard deviation is
.
0
aXlmum Imenslonless t
.
tnangular pulse. All difference I I .
emperature nse for the
ca cu atlons that follow are performed using
S EC.6.11
259
NUMERICAL CALCULATIONS
( /= I and lumped capacitance. Similar results are expected for (J = 1/2 and/or
other capacitance distribution schemes.
Figure 6.9 compares results for exact and inexact sensor data, both with
, = 3 for the function specification method, q = C. The results are quite good
for the exact «(1; < 10  6) temperature data. However, the inclusion o f temperature errors causes scatter in the computed heat flux. The results for (1; =0.0017
are similar to those using Duhamel's theorem that are given in Figure 5.13;
the (1; /11; ratio is 23.2 for Figure 6.9, which is close to the 25.4 value o f Figure
5.13.
Figure 6.10 compares the function specification method ( q=C) for r =4 and
5 for data containing errors. Both values of r yield satisfactory values of heat
flux considering that the d ata contain errors. As r increases, the ability to
follow sharp changes in q is impaired, as evidenced by the calculations near
1+ = 0 ,0.6, a nd 1.2. T he optimal choices o f r a nd at are discussed in Section 5.6.1.
The results shown in Figures 6.9 and 6.10 are nearly identical to those
obtained using the convolution approach, Figures 5.12 a nd 5.13. This is expected because the I HCP algorithms are the same; the method o f solving
the heat conduction equation and the evaluation of the sensitivity coefficients .
are different. The random errors are also different.
The same example problem was also solved using both first and secondorder sequential regularization with CXI =CXl = 0.001 ; these results are shown in
0.6
0.5
0.4
q+
0.3
0.2
o
,,~
< 106
'" ,,~ == 0.0017
0.1
'"
0 .0
 0.1
 0.3
" : == 23.2 for ,,~+ == 0.0017
"''''
"y
00
'" '" '"
'"
0 .0
0.3
0.6
0.9
1.2
'"
1.5
t+
FIGURE 1 .9 Comparison or exact (11; < 1 0 6 ) and inexact (11; =O .OOt 7) d ata ror runction
specification method, q =C, 20 clements, r =3,lu/Al x =2, A I+ = 0.06,6= 1,Iumped capacitance,
x /L= 1.0.
SEC. 8.11
0.6
0.5
0.4
q+
0.3
0.2
0.1
&6
&
0.0
0
0
6
 0.1
 0.3
0.0
0.3
0.6
0.9
1.2
1.5
t+
FI~URE 8.10 Function specification mcthod, q =C, 20 clemcnts, a ; =0.0017, A t/lux =2,
At = 0.06,8= 1,Iumped capacitance, x /L= 1.0.
281
N UMERICAL C ALCULATIONS
Figure 6.11. The firstorder regularization appears superior to the secondorder by a considerable margin. The choice of (X was based on the results in
Chapter 5 for the convolution integral approach. The analysis for the optimum
choice of (X2 has not been performed. However, Figure 6.12 compares results
for (X2 = 5 x 1 0 3 with (X2 = 5 X 1 0 4 and shows very little difference.
Numerous examples were presented in Chapter 5 illustrating how a single
temperature error affects the heat flux calculation at times both before and
after the temperature error. Similar results are shown in Figure 6.13 and were
determined by setting Y = 1.0 a nd all other temperature measurements identit
cally zero. The secondorder sequential regularization exhibits much greater
swings than either the function specification (q = C) o r firstorder sequential
regularization. This partially explains why the secondorder regularization
performed poorly in Figure 6.11. The function specification method has the
smallest swings.
I t has been demonstrated that comparable results can be obtained for both
function specification (q = C) and firstorder sequential regularization methods.
The presence of the regularization parameter (X requires that trial calculations
be performed for each new class of problems solved by the regularization
method. Both the function specification and firstorder regularization methods
have given reasonable results.
0.6
0.6
0.5
0.5
0.4
0.4
q+
0.3
q+
0.3
0.2
0.2
o 1st order,
0.1
6
6
a
0
0.0

0
III
6
0.1
2nd order, 0 2 = 0.00 1
o
8
&
0.0
112
= 5X 1 03
0
= 0.001
02
= 5X 1 04
•
8
III
fI
f!
~~.3~~0~.0nrO~.3'~0.~6~0~.9~~I~.2~~1.5
t+
FIGURE 8.11+ Comparison of first and secondordcr sequcntial rcgularization with «=0.001,
20 clemcnts, a r =0.0017, A t/Atx =2, At+ =0.06, r=5, 8 = 1,Iumped capacitance, x /L= 1.0.
280
0.3
0.6
0.9
1.2
1.5
t+
FIGURE 8.12 Influcnce of «2 on secondordcr sequcntial regularization, 20 clcmcnts,
0.0017, A t/Atx =2, At+ =0.06, r =5, 8 = I, lumped capacitance, x /L= 1.0.
a; =
2 62
CHAP. 6
ONEDIMENSIONAL INVERSE HEAT CONDUCTION PROBLEM
SEC. 6 .12
263
COMPUTER PROGRAMS
T ABLE 6 .8 L ist o f C omputer P rograms f or t he I nverse H eat
C onduction P roblem
O q=C
1st order sequential regularization
10
Program
Name
Developer/
Organization
C ONTA
6
o 2nd order sequential regularization
Beck; Michigan
State Univ.Sandia, Albuq.
Muzzy, Avila,
R oot;GESan
Jose
H CODE
HEATINV
v=[l.i= 1
. 10.
ti
=
all other i
INVERT 1.0
itJ.t
Bell, Wardlaw;
Naval Surf.
Weap. Center
Snider; E GG
I daho
O RMDIM
Bass, Drake,
O tt;ORNL
O RINC
O tt, Hedrick;
O RNL
S ODDIT
Blackwell ;
Sandia, Albuq.
SMICC
Hills, Hensel;
New Mexico
State
Drake; ORNL
Time index. i
FIGURE 6 .13 Dimensionless heat ftux error due to a single measurement error at t + = At + = 0.06.
r =S. 20 elements, At/At x =2.8= l,lumpedcapacitance, x /L= 1.0.
6.12
C OMPUTER P ROGRAMS
A great many computer programs have been written for the inverse heat conduction problem. The purpose of this section is t o mention some of these. A
review o f some computer programs is given by Beck. t 7
A list of some computer programs is given in Table 6.8. All o f the computer
programs, except the last one listed, utilize some difference procedure for
solving the transient heat conduction equation and thus can treat temperaturevariable properties. Each program has unique features and capabilities. Brief
descriptions o f the programs listed in Table 6.8 are given next.
C ONTA was developed by Beck 18 and was supported by both Michigan
State University and Sandia National Laboratories. I t uses the function specification procedure with the temporary assumption o f a constant q. T he sensitivity
coefficients are calculated in a very efficient manner.6 The program can treat
multiple internal sensors in composite plates.
H CODE was developed by General Electric 19 in connection with nuclear
reactor safety analyses. I t can treat composite cylindrical rods. A modification
o f the function specification method is used.
HEATlNV 20 was developed a t the Naval Surface Weapons Center and is
the only differencebased program listed that uses whole domain regularization.
...
ARIES
I HCP
Algorithm
Sequential
function
specification
Modified
sequentiai
function
specification
Oth + l storder
whole domain
regularization
Sequential
function
specification
Sequential
function
specification
Uses
temperatures only
a t time tAl
Sequential
function
specification
Space marching
Oth, 1st, and
2ndorder
whole domain
regularization
Difference
Method
Integral
Method
CrankNicolson
No
CrankNicolson
No
CrankNicolson
No
Generalized
CrankNicolson
No
Finite element
No
Backward
difference
No
Finite control
volume, backward
difference
Yes
No
No
Yes
A combination o f zeroth and firstorder regularization is employed. A sequential regularization procedure, instead of the whole domain method, would be
more computationally efficient.
INVERT 1.0 was written by Snider 21 a t EG&G, Idaho, and uses the function
specification method. I t can treat composite, cylindrical rods.
O RMDIM was developed at O ak Ridge National Laboratory by Bass
et at. 2 2 I t is a twodimensional, finite elementbased program.
O RINC 23 was also developed a t O ak Ridge National Laboratory. I t is for
composite cylinders. The fully implicit method o f solution of the difference
equations is employed. The temperatures only up to time t are used to calculate
the flux a t time t .
2 64
C HAP . 6
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
S ODDIT ~:s deve~oped a t Sand~a N ational .Laboratories, Albuquerque,
by Blackwell. The difference equations are denved using the finite control
v?lume. method. T he sequential function method is used a nd a rbitrary onedimenSional geometries are possible.
i SMI<?C is a ~rog~am de~eloped by Hilts a nd Hensel 1 4 t hat uses space marchng a nd IS d escnbed I n Section 6.10.5. I t also provides estimates for the resolving
power and accur~cy o f the calculated heat flux components. Hilts a nd Hensel
are a~ New. MeXICO S tate University and the research was funded in part b
Sandia NatIOnal Laboratories.
Y
ARIES 1 ' is t~e ~nlr progr~ listed in Table 6.8 t hat is based o n a n integral
m~e~ ~nd t hus IS h mlted t o h near problems. The integral model does give the
flex~bJllty t o tr~at a variety o f inverse problems, however. T he user has the
OptlO~ o f choos~ng ~roth, first o r second (or some combination thereof) whole
d om.un regulanzatlon.
Copies o f some o f t he reports may be obtained from:
National Technical Information Service
U.S. D epartment of Commerce
5285 P ort Royal Road
Springfield, VA 22161
a nd possibly from the authors. In some cases the reports contain listings of the
programs.
PROBLEMS
2 65
I I. Raynaud, M ., D etermination d u F lux Surfacique Traversant U ne P aroi Soumise a U n Incendie
a u Moyen D 'Une M ethods D'inversion, Laboratoire D'Aerothermique G roupe ' Echanges
Thermiques' Universite Pierre e t M arie C urie, Paris, F rance, August 1983.
12. Raynaud, M., a nd Bransier, J., A New Finite Difference Method for N on L inear Inverse Heat
Conduction P roblem, t o be published in Numerical Heat Trans/er.
13. Hensel, E. C. a nd Hills, R. G., A Space M arching Finite Difference Algorithm for the O ne
Dimensional Inverse Conduction H eat T ransfer P roblem, ASME P aper No. 84HT48,
presented a t 22nd N ational H eat T ransfer Conference, Niagara Falls, NY, August 6  8,1984.
14. Hills, R. G . a nd Hensel, E. c., S MICC, t he Space M arching Inverse Conduction Code, S AND
841563, Sandia N ational L aboratory, Albuquerque, N M, 1985.
IS . Richtmyer, R. D . a nd M orton, K.W., Difference M ethods/or Initial Value Problems, 2nd ed.,
Interscience Publishers, New Y ork, 1967.
16. Alifanov, O. M., Inverse Boundary Value P roblems o f H eat C onduction, J . Eng. P hys.l9(I),
8 21830 (1975).
17. Beck, J. V., Review o f Six Inverse H eat C onduction C omputer Codes, A NL/RAS/LWR 811,
Argonne National Laboratory, Argonne, Illinios, February 1981.
18. Beck, J. V., U ser's M anual for C ONTAProgram for Calculating Surface H eat F luxes F rom
T ransient Temperatures Inside Solids, S andia N ational Laboratory, C ontractor Report,
SAND837134, December, 1983.
19. Muzzy, R. J., Avila, J. H., a nd R oot, R. E., Topical Report: Determination o f T ransient H eat
T ransfer Coefficients a nd t he Resultant Surface H eat F lux from Internal Temperature
Measurements, General Electric, GEAP20731, 1975.
20 . Bell, J . B. a nd W ardlaw, A. B., Numerical S olution o f a n I llPosed Problem Arising in Wind
Tunnel H eat T ransfer D ata R eduction, N SWC T R 8232, N aval Surface W eapons Center,
Dahlgren, Virginia 22448, D ecember 1981.
21. Snider, D . M., Invert 1 .00A P rogram for Solving t he N onlinear Inverse H eat C onduction
P roblem for OneDimensional Solids, E GG2068, E G&G I daho, Inc., I daho Falls, I daho
83415, February 1981.
Bass, B. R., Drake, J. B., O tt, L . J., O RMDlN: A F inite Element P rogram for TwoDimensional
Nonlinear Inverse H eat C onduction Analysis, N UREG/CR1709, O RNL/NUREG/CSD/
T M17, U .S. N uclear Regulatory Commission, Washington, D .C., December 1980.
23 . O tt, L. J . a nd H edrick, R. A., O RINC: A O ne D imensional Implicit Approach t o t he Inverse
H eat C onduction P roblem, N TIS R eport O RNLjNUREG23, O ak Ridge National L abora
22.
REFERENCES
~k,
J. V. a nd Wolf, H., T he N onlinear Inverse H eat C onduction P roblem, ASME Paper
SHT40, p resented a t t he A SME/A IChE H eat T ransfer Conference a nd Exhibit, L os Angel'"
CA, August 8 11, 1965.
....,
2. ~OC 7600/CYB~R ~O M odel 76, C omputer Systems H ardward Reference Manual, Control
a ta C orp., P ubbcauon no. 60367200, Minneapolis, 1972.
3. Blackwell, B. F., Efficient Techni~ue for the Numerical Solution o f t he OneDimensional
Inverse P roblem o f H eat C onductIon, Numer. Heat Trans/er 4 ,229  238 (1981).
4. Beck, J. V. a nd Arnold, K. J., Parameter Estimation in Engineering and Science Wil~y New
'
,
York, 1977.
S. Beck, J. V., Surface H eat F lux Determination Using a n I ntegral M ethod N ucl Eng D 7
170 178(1968).
'
.
. es. ,
I.
6. Beck, J. V., L itkouhi, B., a nd St. Clair, C. R., Effective Sequential S olution o f t he Nonlinear
.........
~
...
Inverse H eat Cond~ction P roblem, Numer. H eat Trans/er 5 ,275  286 (1982).
7. d~Boor, C ., A Pracllca/ GUIde to Splines, SpringerVerlag, New Y ork, 1978.
8. ~ Souza, N., Inverse H eat ~onduction P roblem for P rediction o f Surface T emperatures a nd
e at Transfer from I ntenor T emperature Measurements, Report No. SRCR74 S ace
Research C orporation, Montreal, D ecember 1973
'P
D 'So.~ N: Numerical S olution o f O neDimensional Inverse Transient H eat C onduction
.
9. b
y FI~Jte DIfference Method, A SME P aper N o. 75WA/HT81, presented a t W inter Annual
M cetlOg, H ouston, TX, Nov. 3O Occ. 4, 1975.
10. Weber, C. F~ Analysis a nd S olution o f t he IllPosed Inverse H eat C onduction P roblem
I nt. J . H eal Mass 1/'ans/er, 24(11), 1783 1792 (1981).
'
tories, O ak Ridge, T N ., 1977.
Blackwell, B. F., User's M anual for the S andia O neDimensional Direct a nd Inverse Thermal
( SODDlT) C ode, S andia N ational Laboratories internal report o f Div. 7537, 1983.
25. Drake, J . B~ ARIES : A C omputer P rogram for t he S olution o f F irst Kind Integral Equations
with Noisy D ata, K /CSDfTM43, C omputer Sciences a t O ak Ridge Gaseous Diffusion P lant,
Post Office Box P , O ak Ridge, Tennessee 37830, O ctober 1983.
24.
P ROBLEMS
6 .1.
P rove that E q. (6.2.5) is valid for constant properties.
6 .2. Write the difference equations for the solution of Eq. (6.2.7) for the
step sensitivity coefficient Z for constant properties and onedimensional
planar geometry. Use t he time integration of your choice. Demonstrate
that the same difference equations can be obtained by differentiating
Eq. (3.3 .45) with respect to q,.,.
6 .3. I n the algorithm for a single temperature sensor that exactly matches
286
CHAP. 8
O NEDIMENSIONAL INVERSE H EAT C ONDUCTION PROBLEM
the temperature data, can the X sensitivity coefficient be replaced by
Z in Eq. (6.3.8)? Why?
6 .4. Write an efficient computer subroutine. t o solve systems o f linear
algebraic equations that are tridiagonal in form; see Eqs. (6.3.10) and
(6.3.11).
6 .5. Verify the algebra for &is. (6.3.12)  (6.3.14), the algorithm for calculating
the sensitivity coefficients.
6.6. Calculate the sensitivity coefficients Zl.1 a nd Z l.2 for x 1/L=0.5 and
1 for L1t+ = 0.05 a nd 0.1 respectively and for planar geometry with a
perfectly insulated surface. Use the pure implicit method. Compare your
results with the tabular values given in Table 1.1.
CHAPTER
7
M ULTIPLE H EAT
FLUX E STIMATION
N
6.7. Prove that
I
W1 = 1 in Eq. (6.4.11).
1 =1
6 .S. Repeat Example 6.1 for x 1/L= 1 a nd for M + =0.05 a nd 0.5.
6 .9 . I n t he analysis of Section 6.6, the Z (step) sensitivity coefficient was used
in the Taylor series expansion Eq. (6.6.4). Why Z and not X ?
7.1
I NTRODUCTION
6 .10. Repeat Example 6.2 for L1t+ =0.1.
6 .11. Extend the analysis o f Section 6.9 t o include multiple temperature
sensors.
6 .12. C ompare the computational molecules for all the space marching
methods o f Section 6.10.
6 .13. As a h eat transfer consultant, you have been asked t o comment on the
d ata quality from the following inverse problem : A thermocouple is
located 0.001 m below the surface o f a 0 .005 m thick copper slab. The
known back face b oundary condition is a specified heat flux of 50 kW/
m 2 • T he person describing the experiment to you " thinks" the unknown
heat flux you are trying to estimate is o f the order o f magnitude of
1 kW/m2. Will this experiment yield meaningful information?
In the previous chapters the case o f the single unknown surface heat flux history
was considered. I n this chapter the case o fthe multiple heat flux I HCP is treated.
A multiple heat flux case arises when both surfaces of a onedimensional
body are exposed t o unknown heat flux histories. See Figure 7.1a for a plate
heated o n b oth sides. The unknown heat flux histories are q1(t) a nd Q2(t).
Figure 7.1b depicts either a hollow cylinder o r sphere. Again there are two
unknown heat flux histories. In both Figures 7.1a a nd 7.1b there must be at
6 .14. a.
F or t he D'Souza procedure, derive a twonode expression for the
surface heat flux for a solid cylinder with the sensor a t the center.
Derive the difference equations from an energy balance.
b. Repeat p art (a) for three nodes.
c. C ompare the expressions with that obtained using the exact ex pression given by Burggraf.
(a)
6 .15. Solve Problems 6.14 for a solid sphere.
6 .16. F or the Weber procedure, derive a twonode expression for the surface
f iGURE 7 .1 Some onedimensional cases o f
bodies exposed t o two unknown heat lIux histories,
q I (t) and qa(/). (a), plate; (b), hollow cylinder o r sphere.
heat flux for a solid cylinder with the sensor a t the center. The wave
speed is infinite. Derive the difference equations from an energy balance.
Compare the equation with the exact expression given by Burggraf.
6 .17. Solve Problem 6.16 for a solid sphere.
(b)
287