On the iterative regularization of inverse heat conduction problems by conjugate gradient method

Int. Comm. Heat Mass Transfer, Vol. 25, No. 7, pp. 999--1008, 1998

Copyright © 1998 Elsevier Science Ltd Printed in the USA. All rights reserved 0735-1933/98 $19.00+.00

Pergamon PII S0735-1933(98)00091-8

ON THE I T E R A T I V E R E G U L A R I Z A T I O N O F INVERSE H E A T C O N D U C T I O N P R O B L E M S BY C O N J U G A T E G R A D I E N T M E T H O D

M. Prud'homme and T. Hung Nguyen Department of Mechanical Engineering t~cole polytechnique de Montrdal C.P. 6079, Succ. Centre-ville Montrdal, Qu6bec H3C 3A7 Canada

(Communicated by J.P. Hartnett and W.J. Minkowycz)

ABSTRACT The convergence and regularization properties of the conjugate gradient algorithm applied to the inverse heat conduction problem are considered for a time-dependent boundary heat flux. An analysis based on both numerical and analytical results clearly shows that the convergence process of the algorithm is strongly frequency-dependent and provides in this way a very efficient regularization mechanism against the destabilizing effect of random errors in the input data. © 1998 Elsevier ScienceLtd

Introduction The first attempt to solve an inverse heat conduction problem (IHCP) was made by Stefan who obtained an infinite series solution in 1890. Although the formulation and solution of this problem were presented over a century ago, its development has rapidly grown only during the last twenty years, due to a combination of the advent of high technologies, new mathematical achievements and modem computational facilities. An excellent review of litterature and comprehensive bibliography can be found in the books of Beck et al. and Alifanov [1,2]. One of the most stable algorithms for IHCP is the iterative method of optimization by conjugate gradient technique, even without a Tikhonov regularization term in the objective fimctional [3,4,5]. It is therefore of prime importance to show how such a stabilizing effect is i~hherently built in this iterative process of minimization.

999

1000

M. Prud'homme and T.H. Nguyen

Vol. 25, No. 7

To this purpose, we will first solve the 2-dimensional IHCP by conjugate gradient technique using an implicit control volume approach. Results obtained at each step of the algorithm will be presented to illustrate the convergence process. These results will be analyzed next via the analytical solutions of the direct, adjoint and sensitivity equations to show that. due to the diffusive nature of these equations which are the building blocks of the algoritlma, the low-frequency structure of the heat flux is revealed after a few iterations, while high-frequency components can be recovered only at a much later stage. Satisthctory predictions of the unknown heat flux may therefore be obtained from noisy data by using the iteration number as a stopping criterion.

Problem Definition The IHCP under consideration is depicted in Fig.1. The problem is to determine the unknown heat flux

qr(t)

at

x=H

for

O<_t<_tf, given temperature

measurements T~, a flux qr(Y) on the left wall and

q (y,t)

T (t)

q (y,t)

adiabatic conditions on the top and bottom walls at y = 0 , H . The cavity is initially at temperature To when heating begins.

)X

Based

on the

usual

dimensionless

variables, the temperature T must satisfy the heat equation

H FIG. 1 IHCP geometzy and coordinate system c3T = V 2 T

(1)

Ot

Method of Optimization by ConiuRate Gradient The IHCP stated above is treated as an optimization problem : given a set of" N temperature measurements T m over the interval 0 _< t _< t / , find the boundary heat flux q(t) on surface S (the right wall) that minimizes the object functional

E(q)--½1Iv-,112 ---½(T-VmlT-Minimization is achieved by' constructing a sequence of approximations qO,q~ ...,qk

(2) k+~

by' the conjugate gradient method such that qk~ = qk + ockpk where a k is the step size and pk the conjugate search direction [6]. When no a priori information is available on q ( t ) , the minimization of E must be performed in an infinite-dimensional space.

Vol. 25, No. 7

INVERSE CONDUCTION BY CONJUGATE GRADIENT METHOD

1001

The Sensitivity Problem Let us define the temperature sensitivity T as the directional derivative of T at q in the direction Aq, i.e

= limT(q + gAq) - T(q) ~-~0

(3)

E:

It is straightforward to show that 0 ~ = V2 ~

(4)

Ot 0~__~x=~ = Aq(y,t )

(5)

with Tit=0 = 0 and adiabatic conditions on all the remaining walls. These equations constitute the socalled sensitivity problem, whose driving force is Aq.

The Adj0int Problem The gradient VE(q) is formally defined from the directional derivative by

D~E(q)= (T- Tm[:~)= (VEIAq)

(6)

It may be shown based on Eqs.(4), (5) and (6) that the gradient of the object functional is equal to the so-called adjoint temperature T at the surface where the unknown flux is being sought, that is VE = T(S,t). The corresponding adjoint equation may be written as [2,3] N

a_~T = - V 2T + " ~ ' ( T - T~)f(F - F~)

L~

Ot

(7)

i=1

where ~ are the sensors positions. The adjoint temperature must vanish at t = t l , and satisfy adiabatic conditions on all the boundaries. When the sensors are located say, on the left wall, the source terms are replaced by the boundaly condition

O~-xT

= T - Tm

x~0

Note that the driving force of the adj oint problem is the error between the given data and the calculated temperature field.

(8)

1002

M. Prud'homme and T.H. Nguyen

Vol. 25, No. 7

The Coniugate Gradient Method The minimization algorithm by conjugate gradient may be summarized as follows [5].

1. Set initial conditions and choose initial guess q0. Set iteration counter k = 0. 2. Solve the direct problem with q~ to obtain T k . 3. Evaluate the error T k - Tm at the sensor. 4. Solve the adjoint problem backward in time to obtain T k . 5. Evaluate the gradient VE k = T k ( S , t ) . 6. Calculate the search direction 7 k =(VE ~-VE'

~IVE~)/JlVEk

pk.

If k = O , p k = - V E k , o t h e r w i s e p k = - V E k + y k p k

l

vdth

'[I2

7. Solve the sensitivity problem with kq = pk on surface S to obtain ~k at the sensor. 8. Calculate the step size c~k = - ( V E k ]pk)/tl~kll= 9. Update qk+t =qk +ctkpk. 10. Set k = k + 1, go back to step 2, repeat until criterion E k < s , or k > km~x, is satisfied.

It is worth mentioning that since the adjoint temperature vanishes at t = t r , the conjugate gradient method will always give qk ( t f ) = q°(t r ), thereby convergence may be veD' slow if q°(t 7 ) has an unrealistic value at If.

Results and Discussion The direct, sensitivity and adjoint problems are discretized by a first-order implicit, control volume approach with a power-law scheme. The discrete equations are then solved by alternating line and colurml sweeps by the Thomas algoritban at each tilne step. For the adjoint equation, we first make the change of variable r = t j - t to solve for t < t j with the condition at t = t f , The adjoint problem is thus solved as an initial value problem in r with a positive diffusivity coefficient, which ensures

Vol. 25, No. 7

INVERSE CONDUCTION BY CONJUGATE GRADIENT METHOD

1003

numerical stability [7]. All the results presented here are obtained for a 21 by 21 mesh and a time step A t = A x = 1 0 -2. We first consider a 1-dimensional IHCP, with a zero heat flux on the left wall, and try to recover a two-component heat flux

q(t) = - sin(nO - 0.2 sin(107tt) on the right wall, using one sensor at the center

o f the cavity. The results for the first few iterations are shown in Figs 2-4. Although the corresponding T O- T,~ profile reveals the presence o f the high frequency component in the input data, the IT° profile as well as the heat flux qJ do not show-any significant influence of this high frequency component. It should be noted that this first update is far from being zero at t = 0, as the exact flux does. However, the second iteration yields a much more realistic heat flux q2 with a rather recognizable sinusoidal shape in Fig. 4. The lower frequency component o f the flux is thus reasonably well recovered after only two iterations.

0.0

S

-0 2

4o(

-04

3o!

i

!

~°(.l,t)

2oi

-0.6

i[

,ol

-08

A

-1.0 00

02

06

0.4

08

o

o ~'-~-~ . . . .

10

/

t

FIG. 2a

FIG. 3a

Error at tile sensor, first iteration

Adjoint temperature at x = 1, first iteration

i/ O2

-5

IL'I m

/-

01

' T'(1,t)

J

/ l

l

O0

-01

-02 -

-03

O0

02

04

06

08

J

j'0.6

8

~0

o2

---o-0_/o.o

t

FIG. 2b

FIG. 3b

Error at the sensor, second iteration

Adjoint temperature at x = 1, second iteration

1004

M, Prud'homme and T.H. Nguyen

Vol. 25, No. 7

The third iteration starts to uncover the presence o f the high frequency term. Figs.2c and 3c for the T 2 - Tmprofile and adjoint temperature clearly reveal the high frequency component. However, many more iterations would be necessary to fully recover the heat flux component 0.2 sin(lOnt).

T2-Tm Ot)[

0

-0 I hi

O0

oO -001

'~/~

4) 02 (} 0

02

04

06

08

I (I

/~0.6

08"-

(3"~"-" - _ 04

t

0 4 -

- ~ 02

/

02

~-¢ o o

t

FIG. 2c

FIG. 3c

E~oratthe sensor, thirditeration

Adjoint temperature at x = 1, third iteration

Let us now consider a 2-dimensional IHCP by applying a non-tmiform flux on the left wall. We know' that if the unknown flux remains only time-dependent, it can be recovered with a single sensor. Figure 5a shows the result after 20 iterations, using exact data Tm at the center o f the domain, for the botmdary conditions qt = c o s ( n y ) a n d qr = - s i n ( n O .

00

//

-o2 \\ \ -o.4 -

qL(t)

.-•

\

...

\

" ~

X

"0,6

..........

/

.

/

l '

\..," +o8

/

/

q2(t)

/

-°°i

\

\

-to

/ /

\

.... ',

"/

/'0.8

4, o ! . /

-I 2 00

~

04

0.6

(18

' /o 6

0,~ . . .0 . 4." . - . ---._ . . .... . . . t

02

/:

"O 2

--0.5/00

FIG. 4

FIG. 5a

Predicted fluxes q', q2

F l u x after 20 iterations. 2-d case without noise

The regularization property o f the inverse algorithm is studied next by introducing random terms in Tm to simulate the errors encountered in real measurements. Computations done with T~,(l+.02o), where -0.5 <_o -< 0.5 is a random number, are shown in Figs 5b, 5c : the predicted flux q,. appears much smoother after 6 iterations than after 20. We thus recover the true heat flux long before the undesirable

Vol. 25, No. 7

INVERSE CONDUCTION BY CONJUGATE GRADIENT METHOD

1005

higher frequency components of the noise are taken into account. The following section provides some insight into this sequential regularization.

~

q(y,t,

I

o~--~o.-6_._..a:~ . . . . . /o.~ o.~ --;:o_/o.o

~

'

"

"

q()7,t)

°z' o.-~

t

~,~,

/o.~ o.-~-~--~~.~/o. o

t

FIG. 5b

FIG. 5c

Flux after 6 iterations, 2-d case with noise

Flux after 20 iterations, 2-d case with noise

Analysis of the Regularization Mechanism

The regularization effect of the conjugate gradient method is best understood by following the inverse algorithm step by step to recover a heat flux of the basic formq(t) =ei~, assuming adiabatic conditions on the left wall. Using Laplace transforms, the direct problem under the above boundary condition can be solved exactly to obtain n

T(x,t)=

q(t)cosh(xf~x) t ~ sinh(xff~)

-n2~2l

/ _2y(-21~)f2e co ~ g ~ l n r c + i t o

cos(nnx)

(9)

where the denominator xff~ sinh(/~) may also be expressed in short hand notation as De -i# with ~b---~-n+ tan-t(~!tanhxf~)~tanx/co/2 Jl' and

As D~ ~ f ~ e ~

D= xf~[sinh2 x / ~ +sin 2 xf~-] vz

(10)

(at high frequencies), the magnitude of the temperature field will decrease

exponentially as co increases. Neglecting the rapidly decaying transient terms, the temperature measured by a sensor on the left wall is therefore

Tm(O,t) = -q(t) i~

i ~o

(11)

1006

M. Prud'homme and T.H. Nguyen

Vol. 25, No. 7

Let us now start the iteration process to recover q(t) with an initial guess qO = 0. The boundary condition for the first adjoint temperature field is then simply T o - T = - T m . The adjoint temperature at x = 1 will behave asymptotically as ~o(1,t ) = _ q(t)_ + ie '(* .... ) i +--f~ D2 coD co where f~ = t

-

(12)

t f + 1/6. The table below displays the values o f upper bound },-0(1,t), IT

1

1

The values o f If~[ are shown in Table 1 at t = t f / 2

1

(13)

for t / = 1. It is clear that the term involving

f l is the leading term as o is increased. TABLE 1. Values of terms of Eq. ( 13) at two frequencies

0)

l/ O =

1/o)0

If, l~0)

rr

9.1 x 10 -2

9.6 x 10 -2

1.1 x 10 -I

10~

4.6 x 10 -s

2.2 x 10 -4

1.0 x 10 -2

For the temperature sensitivity, we readily obtain from Eq.(4) the asymptotic (for large t) solution

T°(O,t)

=

e'* - q ( t ) D3

i

coD2

ie ~(~+'% ) +--~z~(t-1/6)+if2

co

(14)

where

f2 = t2 - ( t f + 1/6)t + 1/6(tf - 1/6) + 7/180

(15)

is once again the leading term, even more dominant in Eq.(14) than the fz term in Eq.( l 2).

Figure 6 shows a comparison o f the full numerical solution, including transients, for T ° at x ----1 with the asymptotic solution Eq.(12), when co = n . The agreement is satisfactory overall, with greater discrepancy near t ---t r where the transients o f T o are present and to a lesser degree near t = 0 where the transients o f T~, are also present, as we might have expected. Figure 7 shows the same quantities, but for co = 10n. The solution is practically a straight line as the fl term is now dominant. The amplitude o f both solutions (analytic and numerical) is one order o f magnitude smaller than for the lower frequency co = n .

Vol. 25, No. 7

INVERSE CONDUCTION BY CONJUGATE GRADIENT METHOD

025

0 030

020

0.025

015

0.020 - -

010

Nulnerlcal

""

1007

TU(l,t) "

""'

.

0015

0 05

0.010

000

0.005

-005

0 000

-0 l0

-0 005

-0.15

~).010 0.0

02

04

0.6

0.8

1.0

00

02

04

t

0.6

0.8

1.0

t

FIG. 6

FIG. 7

Adjoint temperatureat x = 1, co= n

Adjoint temperatureat x = 1, co= 10~

On the basis of the above fundamental solutions, let us now investigate what happens when two components of different frequencies are present in the unknown flux, for example

q(t)

= - sin(o~t) -

sin(co/et)

(16)

where e is a small number, say e <0.1. The expression for the first update qJ is obtained by superposition from Eq.(12). Using the leading term approximation for the adjoint temperature associated with the higher frequency component of the flux, we readily obtain the asymptotic solution

q' = ct°[T°(1,co,t) - c/¢o f,]

(17)

in which

is the adjoint temperature associated with the lower frequency component. This clearly shows that the main contribution to the flux update is due to the lower frequency component. Furthermore, the step size can be approximated as

~0

2

-IYit=L

II °li

'

IIT°ll 2

+o2]

(19)

For c << 1 the value of ct thus remains virtually unchanged from what would be found if the flux contains only the lower frequency component sin(cot). Proceeding into the next iteration, it is found that T ~ - T,, consists of transients, a quadratic term in t, a third negative power o f D, which in turn lead to a q2 involving a cubic polynomial in t, up to the fourth negative power of D, and so on...We thus find that the leading contribution to the second iterated function q2 is still due to the low frequency component of

1008

M. Prud'homme and T.H. Nguyen

Vol. 25, No. 7

the given data Tin. In fact, the effects of the high-frequency component build up very slowly and can be completely revealed only at a much higher iteration number. For example, while it takes only a few iterations to recover reasonably well the low tYequency component ~ = ~, convergence for the highfrequency component co / e = 107t is not yet attained after more than 30 iterations.

Conclusion A study was made of the regularization mechanism of optimization by conjugate gradient. An analysis of the iteration process based on the analytical solutions of the direct, adjoint and sensitivity equations agrees well with the numerical solutions to show that, due to the diffusive nature of these equations, the first few iterations only reveal the main structure of the boundary heat flux with the lowest

frequency. Higher fi"equency components are recovered at a much later stage. The conjugate gradient approach was thus proved to have a sequential filtering mechanism that can effectively eliminate the random noises in the given data by appropriately using the iteration index as a regularizing parameter.

Acknowledgements This research was supported by the Natural Sciences and Engineering Research Council of Canada.

References 1. J.V. Beck, B. Blackwell, and C.R. St-Clair Jr, Inverse Heat Conduction. Ill-Posed Problems, WileyInterscience, New-York (1985). 2. O.M. Alifanov, Inverse Heat Transfer Problems, Springer-Verlag, Berlin Heidelberg (1994). 3. T.H. Nguyen, Optimization Approach to the Inverse Convection Problem, Proc. Intl. Workshop Inv. Probl., Ho-Chi Minh City,, pp. 83-90 (1995). 4.

T.H. Nguyen, An Optimization Approach to Some Inverse Convection Problems, Proc. Intl. Conf. Analysis andMech, of Continuous Media., Ho-Chi Minh City, pp. 188-195 (1995).

5. M. Prud'hounne et T.H. Nguyen, Num. Heat Transfer, Part A, 32, 169 (1997). 6. R. Fletcher. Practical Methods of Optimization, John Wiley & Sons. New York (1987). 7. M. Prud'homme et T.H. Nguyen, M6thode de solution par volumes de contr61e du probl6me inverse en convection naturelle, Rept. EPM/RT-96/10 (1996).

Received April 30, 1998