Two dimensional steady-state inverse heat conduction problems

Int. Comm. HeatMass Transfer, Vol. 25, No. 4, pp. 541-550, 1998

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

Pergamon P I I S0735-1933(98)00041-4

TWO DIMENSIONAL STEADY-STATE INVERSE H E A T C O N D U C T I O N P R O B L E M S

N. M. AL-Najem, A. M. Osman, M. M. EI-Refaee and K. M. Khanafer Mechanical and Industrial Engineering Department Kuwait University, P. O. Box 5969 Safat 13060 Kuwait

(Communicated by J.P. Hartnett and W.J. M i n k o w y c z )

ABSTRACT In this work we estimate the surface temperature in two dimensional steady-state in a rectangular region by two different methods, the singular value decomposition (SVD) with boundary element method (BEM) and the least-squares approach with integral transform method (ITM). The BEM method is efficient for solving inverse heat conduction problems (IHCP) because only the boundary of the region needs to be discretized. Furthermore, both temperature and heat flux at the unknown boundary are estimated at the same time. The least-squares technique involves solving the equations constructed from the measured temperature and the exact solution. The measured data are simulated by adding random errors to the exact solution of the direct problem. The effects of random errors on the accuracy of the predictions are examined. The sensitivity coefficients are also presented to illustrate the effect of sensor location on the estimated surface conditions. Numerical experiments are given to demonstrate the accuracy of the present approaches. © 1998 Elsevier Science Ltd

Introduction The direct heat conduction problems are concerned with the determination of temperature distribution inside the solids when the boundary and initial conditions, heat generation rates and thermophysical properties are specified. In contrast, the inverse heat conduction problems (IHCP) are concerned with the estimation of surface conditions, heat source rates, and physical properties by utilizing measured temperatures within the solids. Inverse problems are encountered in various branches of science and engineering. In heat transfer field the inverse heat conduction analysis is used to estimate the surface conditions such as temperature and/or heat flux from measured temperatures. Furthermore, the inverse analysis can be

extended to include the estimation of thermal properties of solids, heat transfer

coefficients, energy generation, and so forth. In addition, one can also envision inverse problems of radiation, phase-change, and other parameters associated with heat transfer studies. 541

542

N.M. AI-Najem et al.

Vol. 25, No. 4

A literature review reveals that a very limited number of studies [1-3] considered the BEM in solving IHCP. Kurpisz and Norwak [1] used regularized BEM with future temperature to obtain stable solution for the estimation of the surface conditions. Lesnic et al. [2] introduced a regularized least-squares and minimal energy method into the BEM formulation to analyze one-dimensional, linear, unsteady IHCP in a slab geometry. The oscillations in the inverse solution resulted from the measurements error are substantially reduced by the minimal energy technique. In contrary to the other numerical methods, BEM does not require any discretization of the internal, only the external boundaries are need to be discretized. Furthermore, both temperature and heat flux at the boundary can be estimated at the same time. In this study the solution of two dimensional steady-state inverse problem with constant physical properties is analyzed by two different approaches. For the sake of simplicity constant and linear elements are considered in this investigation. The least-squares error is used to recover the surface temperature variations from the knowledge of temperature measurements. Numerous numerical experiments are carried out to examine the stability and accuracy of both techniques.

Mathematical Formulation

Consider a steady-state problem confined in a rectangular region 0 _
in 0 < ~ < 1 ,

0(Gn) = F(n)

at ~ = 0

(lb)

~---=0

at ~ =1

(lc)

O(Grl)=O

at r l = O

(ld)

--=0

at r l = B

(le)

at.

0
(la)

The measured temperatures at N discrete locations are 0(~',rl)=Y i

at ~" =1 and rl = rlj,

j=l,2 ...... N

(if)

Vol. 25, No. 4

INVERSE HEAT CONDUCTION

543

_Bo.undaryElement Method (BEM) Formulation

The governing equation for steady state heat conduction in a general domain D given in Eqs. (la)-(le) is transformed to a boundary integral form by the application of the weighted residual as explained in Ref. 7. Both temperature and heat flux are approximated as I(V20)0"dD = ~(q - ~)0'dD - ~(0 - O)q'dD D

aD 2

aO

,

(2)

OD t

3o"

where q = - - and q = . Here 6 is the applied temperature at the boundary 0D t and ~ is the On On applied heat flux at the other boundary aD 2 and $D 1 and aD 2 are the boundaries at which temperature and heat flux

are prescribed, respectively. The weighting function 0" must have continues first

derivatives, and later this function will also be required to satisfy the governing equations. Integrating the first term of Eq. (2) twice by parts, we obtain J'0(V20")dD= - J'~0*dD- ~q0"dO+ J'0q'dD+ I 0 q ' d D D

OD 2

dD t

(3)

aDl

,9D2

The function 0* is now taken as the fundamental solution to Laplace's equation, which is the solution corresponding to a unit point source, i.e. V20 *+Ai = 0

(4)

where A~ is the Dirac delta function and i the point where the source acts. Eq. (4) has the following property: I0(V20 * + Ai)dO = ~0V20*dD+ [0AidD = [0V20"dD +0 i D

D

D

(5)

D

Then E q (3) becomes 0'+ J0q'dD+ Jgq'dD= OD 2

0Di

I0"~dD+ Iq0"dD 0D2

(6)

ODt

where 0 i is the value of the unknown function, at the point where the source is being applied. The fundamental solution to Laplace's equation in two-dimensional case is known [3] to be: 0* = 1 in(l) 2rr r

(7)

where r is the distance between the point where the source is applied and the point under consideration. Eq. (7) is valid for any point in the domain, but it has to be taken to the boundary in order to be used in the boundary solution. Therefor, the boundary equation can be written as follows [3] 1 0 . + i 0 q ' d D = Iq0"dD 2 D O

(8)

This equation will now be applied on the unknown boundary of the domain under consideration.

544

N.M. AI-Najem et al.

Vol. 25, No. 4

Dis.cretization of the Boundary

The subdivision of the boundary can be done using different types of elements. In the present investigation, constant element and linear element are considered. The resulted discretized boundary integral equation at node k has the following form N¢ ~. H ke0 k e=l

m kt k = Y'.G q

(9)

e=l

where m .= N e for constant element and m = 2N e for linear element.

Evaluation of Integrals

For constant element, the components of the coefficient matrices I t and G are given by Hkt=Jq*dD

when k e g

(10a)

D

Hke=lq'dD+ / D 2

when k = g

Gke=10"dD

(10b) (10c)

D

A four-point Gaussian quadrature rule is employed for the numerical integration of the components of the coefficient matrices H and G. However, for calculating the diagonal components, special care must be taken. The diagonal terms H kk vanish due to the orthogonality of the element coordinate and the normal. The terms G ~' can be calculated analytically and it is given by the following equation

G kk= L o tin(2__) + 1] 2n

(l 1)

L~

where L c is the element length. The values of the off-diagonal coefficients of H and G can be written as

[31 Hk e

4

1 Lws ff ~(Ri - R2)2 + (S~ - $ 2 ) 2

4 1 4(RI - R 2 ) 2 +(S I - $ 2 ) 2 G ke = Y~ln(. - - ) w s s=~ d s 2

(12a)

(12b)

where Rt, Rz, St, and $2 are the coordinates of the extreme points of each element, w~ are the weighting for each point, L is the distance from the collocation point to the line element tangent to the element and d is the distance from the collocation point to the Gauss integration points on the boundary element. For linear elements, the main difference is that the G kk elements are assembled in a (N~ x 2 N , ) matrix

Vol. 25, No. 4

INVERSE HEAT CONDUCTION

545

instead of (N c x No). This is because two possible values of flux are considered at each node. The diagonal terms in I t are computed implicitly given as H kk = - ~ H ke

for k ~ g

(13)

t=l

Singular Value Decomposition (SVD) Technique The linear algebraic system of Eq. (9) is ill-conditioned and an appropriate method is needed to treat this system. In this investigation, the singular value decomposition SVD technique is used to solve this system of equations. In the SVD method, the matrix A can be written as the product of a columnorthogonal matrix U, a diagonal matrix W with positive or zero elements, and a transpose of an orthogonal matrix V [4]. The conditioning of a matrix is defined as the ratio of the largest of the the smallest of the

wj's where wj's are the elements of matrix W.

If the condition number is too large, the

matrix is ill-conditioned and its solution by SVD can be obtained by zeroing the small (i. e. replace 1/wj

wj's to

wj's of matrix W

"s by zero if wj ~ 0) and using the following modified equation

X = V. [diag.(l / Wy)]' (UV" F)

(14)

where X is the solution vector.

Inverse problem analysis Our objective in this inverse analysis is to estimate the applied surface temperature F(1] ) from the knowledge of temperature measurements Yj taken along the rlj (j=I,2,....N) direction at location ~ =~*. To estimate the unknown surface temperature F( rl ), we represent this function without loss of generality by the following quadratic function 2

F(~)= Z ~pn p

05)

p=O

where ~bp's are unknown parameters to be estimated by the inverse solution.

The solution of the direct problem defined by Eqs. (la)-(le) is obtained by the application of the integral transform method ITM [5] as

546

N.M. AI-Najem et al.

i=l

J-Li c o s h t J ~ i )

)

=

Vol. 25, No. 4

~1.~ c o s h 0.1.i)

(16a)

+ { 4 i=1]~[(-l)i+l B~ti - 1] e d° s hcosh(Ix ~[ l ~ i2( l - ~i)] sin(l~irl)} t ~ where I1i's are the eigenvalues given by (2i - l)n rt i = - , 2B

(16b)

i=1,2,3 .......

The inverse problem is mathematically ill-posed in a sense that the existence, uniqueness and/or stability is not yet insured. A successful solution of an inverse problem generally involves the transformation of the inverse problem into a well posed-approximation solution [6].

Least-Squares Approach (LSA)

To solve the inverse problem we require thatthe estimated temperature 0 ( ~ ' , r l j , ~ p ) a t t h e sensor location computed

from the solution of the direct problem by using the estimated values of the

coefficients d~p'S should match the measured temperatures Yj over the measurements domain in the least-squares sense. One way to realize such matching is to require that the standard least-squares norm is minimized with respect to each of the unknown parameter. The traditional least-squares criterion S(ff )is defined by N

S(d~p )= ~[Yj - Oj(q~p)]2 .

(17)

The next step in this analysis is the minimization of S(~p) by differentiating it with respect to the unknown coefficients and set the result equal to zero. We find c~S

8~ o

-

N o~j

= 0

p = 0,1,2

(18)

J=lO9 o

Now the inverse problem is reduced to that of solving the above normal equations resulted from the ^

minimization process, Eq. (18) to estimate the unknown coefficients d~p'S associated with the inverse problem. Knowing these coefficients, the applied surface temperature F(T1 ) is calculated from Eq. (I 5).The sensitivity coefficient ZI, is the first derivative of the temperature with respect to the unknown parameter (i.e., Z~p = .~-~-). It represents the changes in 0(~',q,d~p) with respect to the changes in the c~ o

Vol. 25, No. 4

INVERSE HEAT CONDUCTION

547

unknown parameters dpp's. It is preferable to have large and uncorrelated values of the sensitivity coefficients Z ¢ . For optimal estimation the, sensors should be placed at such locations where the readings are most sensitive to changes in the unknown parameters. A small value of Z¢, indicates insensitivity of the dependent variable to changes in the value of the unknown parameter; for such cases the inverse analysis become very sensitive to measurement errors and the estimation process become difficult. To establish the best sensor locations, measurements intervals, and so on, it is desirable to examine the effects of measurement locations, measurement intervals, and similar parameters on the relative values of the sensitivity coefficients. If the sensitivity coefficients are functions of the parameters to be estimated, then the resulting estimation problem is nonlinear; conversely, if they are not functions of the parameters, the estimation problem is linear. The sensitivity coefficients Zcp depends only on the geometry and

the boundary conditions for the problem considered. The sensitivity

coefficients for two different locations are plotted in Fig. 1. Clearly, the sensitivities at ~ = 0.5 are larger than the sensitivities at ~ = I as expected. Also, the sensitivity Z¢, reaches a constant value at about rl = 1.4 while the others continue to increase until about 11 = 1.8. Notice that the sensitivity functions cross each other at about r1=0.6 and 0.8. Furthermore, this figure indicates that the sensitivity functions are linearly independent (having different shapes) in the region between 11= 0.4 to 11 = 1.8. Therefore, for simultaneous estimation of ¢0, ¢~ and ¢2 most of the measurements should be concentrated in the interval 0.4 < rl < 1.8. 3.07 --

2.5-

¢=0.5

.... ¢=1.0 2.0Z,~1.5 1.00.50.0 0.0

I 0.2

[ 0.4

I 0.6

I 0.8

I 1.0

I 1.2

I 1.4

FIG. 1 Sensitivity coefficients at ~ =0.5 & 1.0

I 1.6

I 1.8

2.0

548

N.M. AI-Najem et al.

Vol. 25, No. 4

R e s u l t s and D i s c u s s i o n s

Next

we present the results of the steady state problems in rectangular region with B=2. The

temperature measurements Yj are generated from the solution (16) by choosing proper coefficients ~b's at ~ = 1 with Arl = 0.2. These locations being the farthest from the active surface (~=0). The results obtained under such condition provide the strictest test for the validity and stability of any approach. To simulate an experimental data we introduced random error to the exact data. The random errors in the measured temperature are additive, uncorrelated and have zero mean with standard deviation of 2.5E-02. Before we present our results, it is very useful to define some statistical measures such as the standard deviation c~ and the root mean squared error RaMS to evaluate the accuracy of the predictions by the two different approaches. The standard deviation a of parameter 13 is define by

while the root-mean-squared-error of 13 RMS(13) is given by

where 13, ~ are the true and the estimated values, respectively and ~ is the mean value given by l N^

The bias (D) can be calculated from the knowledge of both RMS and standard deviation ( ¢s ) as D = x/RMS 2 - ¢s2 It should be

(22)

pointed out that when the bias is equal to zero, the estimator is very sensitive to

measurement errors in inverse problems [5]. Figure 2 shows the estimated function 0 with two different methods along with the true value of for the case 0 =I.

The least-squares approach (LSA)

gives better results than the singular value

decomposition (SVD) results except at the region near rl = 0 where there is a maximum difference of 4.5%. The RMS of the LSA is 4.889E-02 and that of the SVD is 7.024E-02. On the average, the LSA estimates the function better than the SVD. Similar calculations are made to compare the true function and the estimated functions for the cases 0 = 11and 0 = l + rl + 112 using the LSA and the SVD. Results show that both methods deviate from the true values near 1] = 0. The RIMS values for LSA and SVD are 3.55E-02 and 3.77E-02, respectively. Also, it is noticed that as rl values increase, the LSA estimates the function better than the SVD. The residuals (ej = Yj - ~'j ) for the case 0 = rl are very small and the

Vol. 25, No. 4

INVERSE HEAT CONDUCTION

549

calculated standard deviation for the residuals are ~(e)=1.271E-02 and 1.478E-02 for LSA and SVD, respectively. The residuals for the other cases 0 = 1 and 0 = 1 + 1"1+.q2 show similar results. The statistical calculations show that the SVD has smaller bias but larger variance than the LSA which means that the SVD is more sensitive to measurement errors than the LSA.

2.0 Applied . . . . . BEM(Constant Element) . . . . . I'IT

1.81.61.41.201.0 0.8-

0.60.4-0.2. 0.0 0.0

I

1

0.2

0.4

0.6

I

I

I

I

I

I

0.8

1.0

1.2

1.4

1.6

1.8

2.0

FIG. 2 Estimated surface temperature for 0 = 1

Conclusions

In this work, singular value decomposition SVD with the boundary element method (BEM) and the least-squares approach LSA with integral transform method (ITM) are presented for estimating the surface temperature in rectangular conducting material. In the LSA, the variation of the unknown surface temperature is represented by a quadratic function with unknown coefficients. Then, the least-squares error is used to recover the surface temperature over the measurements domain. Constant and linear elements are used in the BEM analysis. In this investigation, the singular value decomposition SVD technique is used to solve the resulted ill-conditioned system of equations. Some statistical measures are applied to examine the stability and accuracy of

both techniques. The numerical experiments

demonstrates that both approaches considered in this study are capable to estimate the unknown surface temperature with high degree of accuracy from measurements contains random errors.

550

N.M. AI-Najem et al.

Vol. 25, No. 4

Nomenclature

B

Dimensionless variable, B=b/a.

F(rl)

Dimensionless unknown surface temperature, F(rl) =

Yj

Measured temperatures, Yj - yj - T,o

~'j

Estimated measured temperature.

~1

Dimensionless variable, rl =y/a. dimensionless variable, ~ =x/a

0

Dimensionless temperature, 0

f(y) - T.~

T~

T~

T(x, y) - Too

T~

References

1. K. Kurpisz and A.J. Nowak, Anal Bound Elem. 10, 291 (1992). 2. D. Lesnic, L. Elliot and D.B. Ingham, Int J. Heat Mass Transfer 39, 1503 (1996). 3. C.A. Brebbia and J. Dominguez, Boundary elements- An introductory course, Comp. Mech. Publications, McGraw-Hill (1989). 4.

W.H Press, B.P. Flannery, S.A. Teukolsky and W.T. Vettering, Numerical recipes, Cambridge University Press (1986).

5. M.N. O zisik, Heat conduction, New York, John Wiley and Sons (1993). 6. J.V. Beck, B. Blackwell and C.R.St. Clair, Jr., Inverse Heat Conduction, New York, John Wiley and Sons 1985. Received November 5, 1997