WS.rme- und Stoffffbertragung 26, 1- 5 (1990)

W~rmeu n d Stoffiibertragung © Springer-Verlag 1990

Determination of heat transfer coefficient using transient temperature response chart R. C. Mehta and T. Jayachandran Trivandrum, India

Abstract. An inverse finite element computer code is developed to predict surface heat flux and surface temperature in conjunction with the measured thermocouple temperature history. Determination of convective heat-transfer coefficient and combustion gas temperature is carried out employing transient temperature response chart. Examples are illustrated which are typical of those that arise in thermal design of rocket nozzle. The results demonstrate that the method is remarkable in its ability to estimate unknown boundary conditions.

Superscripts M i k T

present time step iteration step index time step index transpose of matrix

1 Introduction Bestimmung des W/irmeiibertragungskoeffizienten mR Verwendung yon Aufzeichnungen der transienten Ansprechtemperatur Zusammenfassung. Es wurde ein inverses Finite-Elemente-Computer-Programm entwickelt, urn den Oberfl/ichenw/irmestrom und die Oberflfichentemperatur aus dem mit Thermoelementen gemessenen Temperaturverlauf zu bestimmen. Die Bestimmung des konvektiven W/irme/ibertragungskoeffizienten und der Verbrennungsgastemperatur ist mit Verwendung yon Aufzeichnungen der transienten Ansprechtemperatur ausgeffihrt worden. Fiir die thermische Auslegung yon Raketendfisen werden einige Beispiele dargestellt. Die Ergebnisse beweisen, dab dieses Verfahren bemerkenswerte Ffihigkeiten besitzt, um unbekannte Randbedingungen ann/ihernd zu berechnen. NomencLature Bi c Fo h k L qw r S To Tg T~ T t x Y 0 q5 At c~ Q

Biot number, h L/k heat capacity of the material Fourier number, e t/L 2 convective heat-transfer coefficient thermal conductivity thickness of material surface heat flux number of future temperature least square function initial temperature combustion gas temperature surface temperature nodal temperature time abscissa thermocouple data dimensionless temperature, (Tg- T~)/(Tg-- To) finite element function time step thermal diffusivity density

Advanced technological applications such as thermal design of rocket nozzle, heat shield of launch vehicle, gas turbine blades, analysis of quenching processes and development of composites for high temperature applications require transient convective heat-transfer studies. To determine the necessary heat transfer data, one has to experimentally measure transient temperatures in an instrumented prototype or scale model. For bodies with high temperatures and subjected to an abrupt change in the convective boundary conditions, it is difficult to measure directly the surface heat flux or the surface temperature. The unknown surface conditions can be conveniently predicted from transient thermocouple measurements at anappropriate point inside or outside surface of a heat conducting material. This type of problem is referred as a mathematically ill-posed problem, because the boundary conditions on one of surfaces are not specified in the boundary val~ae problem. Inverse heat conduction problem may be solved using an exact solution [1], an integral method [2], a finite difference or a finite element methods [3-5]. In the previous analysis [6], an analytical solution was developed to determine convective heat-transfer coefficient for constant heat flux condition. A sufficiently powerful method to analyse the general problem appears to be numerical technique, although the basic concepts can be applied to the integral method. Solution of the linear inverse heat conduction problem can be obtained without iteration if smaller time step is used in numerical analysis and also thermocouple probe is placed at an appropriate location in order to avoid time delay in thermal response. Beck [7] makes this possible by using a procedure that involves minimization of the sum of squared

2

W/irme- und Stoff~bertragung 26 (1990)

difference between the actual and calculated temperatures at the known thermocouple location. The Newton-Raphson method [8] is found suitable when a higher initial time interval for starting the solution is employed in computation of unknown boundary conditions. The purpose of the present paper is to develop two different kinds of algorithms for predicting wall heat flux/surface temperature (unified thermal problem). This unified thermal problem is coupled with a transient temperature response chart in order to determine convective heat-transfer coefficient and combustion gas temperature. As well be shown, the graphical approach is straightforward and economical.

where the r,, are the N nodal values of the temperature. Following the usual Galerkin method, one obtains the heat conduction equation in the form of a system of semi-discretized ordinary differential equations in terms of nodal temperature, collected in the gobal vector T as C T+K r=o

(6)

where T is the time derivative of the nodal temperature vector, and C, K and g are the conventional heat capacitance matrix, thermal conductive matrix and heat flux load vector, respectively. Typical expressions for the matrices in Eq. (6) are L

C=~eE~*~OmdX m

2 Problem description

0

L

A temperature sensor is used at the insulated boundary at X = L and a transient heat flux is imposed on the other surface at X = 0 as depicted in Fig. 1. The temperature distribution inside the solid body of finite slab, assuming that the material is homogeneous and isotropic, is described by the heat conduction equation ~2T aX 2

i aT ~ ~t '

0
t>0

(1)

aT(0, t) = qw(t), aX

~T(L, t)=0,

t> 0

(2)

t>O

m

0

g~ =qw(t),

gin--0,

m--2, 3 . . . . . N

(7)

where the summation is over all the elements, q~,~is the local interpolation matrix for each element, and Bm defines the temperature gradient within each element as a function of nodal values. Using a modified Galerkin integration scheme for time integration procedure for the time derivative in Eq. (6) yields (C + 2At. K) ~ +1 = (C - ½At- K) T, + g At

subjected to the boundary conditions

k

K=k E I BT BmdX

(8)

where the subscripts denote the corresponding time station. This algorithm possesses unconditional stability and accuracy, and is also free from oscillatory behaviour.

(3)

aX

3 The Newton-Raphson method

and the initial condition T (X, 0) = To

for all X.

(4)

The transient heat flux q~(t) in Eq. (2) is the unknown function to be estimated by utilizing the temperature measurement Y (L, 0- Equation (I) may be solved using standard finite element procedure [9] where the geometric domain and the attendant temperature are divided into number of nodes such that

T(X, t)= T.,(t) ~bm(X, t),

qw

~

'1 2 3 . . . .

m = 1, 2 . . . . . n

/-U . . . . . . . . . L Thermocouple z j~_i'fN '' ~,-[ocation

(5)

The solution of the complete problem cannot be readily obtained because the q~.(t) is now known. The transient wall heat flux can be computed using an iterative technique. A special procedure based on the Newton-Raphson method and measured thermocouple data are employed to predict unknown wall conditions. During each time step k=n + i, initial guess is made for the nodal heat flux qw in Eq. (8). The temperature distribution is evaluated in the solid domain using this boundary condition for the solution of Eq. (8). The temperature calculated at nodal points corresponding to the location of thermocouple (denoted by F) is then compared to the temperature measured experimentally (denoted by D). The Newton-Raphson method is used to update the nodal vector qk to minimize the difference between F and D using the updating formula. (qk)i+ 1 = (qk)~ ((j- 1)k)i(((F)k)i_ D)

..f--m~

X=O

X=k

I --,.x

Fig. 1. Geometry of heated surface

(9)

where the superscript k stands for the time step and i for the nodal iteration number. The Jacobian is in numerical form

((F)k)i_ ((F)k)i- 1 (sk)i_ ((q)k), ((q)k)i- 1 "

(10)

R. C. Mehta and T. Jayachandran: Determination of heat transfer coefficient using transient temperature response chart The new heat flux (qk)~+~ is used to determine the new temperature distribution in the solid domain and iteration procedure is continued until the difference between the calculated and measured temperature data at a given thermocouple location converges to a predetermined tolerance; 1. E-4 say. The above mentioned minimization scheme is found suitable when a single thermocouple is used to measure the temperature history and a large time step in the computation is employed in order to satisfy the larger time interval of measurement. One disadvantage of this scheme is that it needs iterations to predict the unknown surface conditions.

4 The sequential future-information method The future surface heat flux qM can be estimated without iteration using additional future temperature measurements. The maximum amount of information about the surface heat flux is found using small time steps. The Beck procedure [7] is to predict at each time step the heat flux that minimizes the least squares error between the measured and calculated temperatures. The least squares function r

S= ~ (YM+J-- TM+J)2

(11)

j=l

is minimized with respect to unknown average heat flux qM between the time step At. yM+j and T M+j are respectively measured and calculated temperatures at time (M + j ) A t . The number of future temperatures used in the minimization is r. A solution for the unknown heat flux is given by Beck et al. [11] as 1

q~U=qM-t+AM ii=i (YM+~/-*-- TM'J -*)

8T ~ + j - 1

~qM

0.8

0.6

o

.5 0.4

(12a)

3

where

J:*\

aqM

7"

(12b)

The difference between qM_ qM- 1 represents a correction in the heat flux. The above notation T ~t'; means that T is evaluated by using qM-i for time t ~ - 1. The partial derivative in Eq. (12) is called a sensitivity coefficient and it is approximated by using a forward difference formula. For the linear inverse heat conduction problem, the sensitivity coefficients do not have to be updated for each qM because they are invariant with M. Thus, they are required to be computed at the starting of the computations and stored in the memory. The Newton-Raphson iteration method and the sequential future-information method described above are efficient for solving linear inverse heat conduction problems and can be implemented on microcomputers with small memories.

5 Determination of heat-transfer coefficient The transient temperature response chart as depicted in Fig. 2 is utilized here to determine the convective heat-transfer coefficient and combustion gas temperature in conjunction with the above mentioned unified thermal problem. The following simple procedure can be adopted for determination of Newtonian heating conditions: Step 1. Compute Fourier number, Fo=c~ t/L2 Step 2. Assume combustion gas temperature, To Step 3. Using previously predicted value of T~, compute 0 = ( T o - T~)/(7~-- To) Step 4. Find in Fig. 2 Biot number, Bi = h L/k corresponding to the values of F o and 0 Step 5. Compute %(t)=(BiL/k) (Tg-T~) and compare it with the previously predicted value of qw(t) Step 6. If it is not agree, change the value of Tg and go to Step 3 Step 7. If it coincides, repeat the above steps for next time level.

I

6 Examples

b...~

0.2

--01 0

1 2 Fig.2. Transienttemperatureresponse

3

1. Finite element method, iteration technique and graphical approach discussed in the above sections have been utilized for estimation of wall heat flux, surface temperature, heattransfer coefficient and combustion gas temperature for a typical divergent rocket nozzle made of mild steel in conjunction with the experimentally measured outer surface temperature data in a static test. The nozzle conditions and material properties taken are: L=0.0211 m, To=300K, c~= 8.1291 x 10- 6 mE/s, k = 35 W/m K (average) and burning time = 16 s. The heat input consists of a convective heat flux at the inner surface X = 0. The X = L surface is insulated. An

4

W/irme- und Stoffiibertragung 26 (1990) Thermocoupie location

Table 1. Estimated boundary conditions t

Fo

[s] 6 7 8 9 10 1t 12 13 14 15 t6

Y

0

Bi

Tg x 103 [K]

qw x 106 [W/m 2]

0.68 0.65 0.65 0.64 0.62 0.67 0.65 0.64 0.62 0.64 0.64

1.58 1.09 1.12 1.14 1.18 1.23 0.89 0.91 0.92 0.91 0.93

3.35 3.11 3.08 3.11 3.18 3.06 2.98 2.96 2.92 2.94 2.93

5.25 3.29 3.29 3.29 3.28 3.27 2.38 2.41 2.36 2.36 2.36

[K] 0,109 0.128 0.146 0.t64 0.182 0.201 0.219 0.237 0.255 0.274 0.292

326 342 356 380 402 425 440 460 479 507 528

~\\\\\\\\\\\\'%%%%%N~ .~

r oozzle Secondary nozzle

Fig. 3. Static test Rocket motor assembly Table 2. Estimated boundary conditions t

Fo

Is] 0.5 1.0 1.5 2.0

Y

0

Bi

0.1 t 3 0.226 0.338 0.451

315 343 384 428

To x 103

[K]

[K] 0.90 0.87 0.83 0.82

0.32 0.31 0.38 0.33

650 K

2.80 1.92 2.00 1.98

insulated chromel/alumel thermocouple of 30 gauge was used for measuring the temperature history at the outer surface of the nozzle. A time delay of 6 s was noticed in the measurement of temperature history. Due to time lag in thermal response, the Newton-Raphson method is employed to estimate unknown boundary conditions. Twenty one nodes and a time interval of I s are taken for computational purposes. Table I depicts predicted values of wall heat flux, surface temperature, convective heat-transfer coefficient and combustion gas temperature. The values of heat-transfer coefficient and combustion gas temperature determined using the graphical method coincides with the previously computed results [12]. Thus, it reveals that the graphical approach is accurate, simple and economical to obtain Newtonian heating conditions. 2. The sequential-future information method is now used in conjunction with the finite element method to solve unified thermal problem. The transient temperature response chart is employed to determine convective heating condition. A schematic sketch of system under consideration is illustrated in Fig. 3 along with a thermocouple location. The nozzle conditions and material properties taken are L =0.006 m, To=303 K, ~=8.129t x 10 -6 mZ/s, k = 3 5 W / m - K (average), burning time = 3.0 s. Eleven nodes, a time step At = 0.1 s and the future temperature r = 6 are used for computational purposes. This linear inverse heat conduction problem, in turn, utilized these transient measured temperature data to predict unknown surface conditions. Measured outer thermocouple temperature and estimated inner surface temperature histories, and predicted transient surface heat flux are shown in Figs. 4 and 5 respectively. Using the transient tem-

Beck's method

58O

~, 510

8.

E

440

370

~ / / ~ l e m p e m t u r e

300

0.6

1.2 Time

,

T 1.8

,

s 2.4

Fig. 4. Temperature history

5.0 W/m 2 4.Q

3.0 2.0 2E

1.0

00

0.4

0.8

1.2 Time

Fig, 5. Surface heat flux history

1.6

2,0 s

2.4

R. C. Mehta and T. Jayachandran: Determination of heat transfer coefficient using transient temperature response chart perature response chart, values of convective heat-transfer coefficient and combustion gas temperature are tabulated in Table 2. It may be observed that the present procedure is able to estimate sudden and large variations in heating rate. It would be worthwhile to mention here that one can obtain the Newtonian heating conditions using transient temperature response chart at any selected time interval without solving again the complete heat conduction equation.

7 Conclusions Two different types of parameter estimation and a graphical procedure are presented for estimation of transient heating conditions when temperature history at the outer surface of the heat conducting material is made available. The examples of the linear inverse heat conduction problem show that the Newton-Raphson method or the sequential future-information method and the graphical approach are having ability to handle abrupt and steep large changes in heating rate. The methods are efficient, convenient and economical and appears to be particular attractive with microcomputer. It has the advantage that convective heat-transfer coefficient can be directly determined at any specified time employing the transient temperature response chart.

References 1. Burggraf, O.: An exact solution of the inverse problem in heat conduction theory and applications. J. Heat Transfer 82C (1964) 378-382

5

2. Beck, J. V.: Surface heat flux determination using an integral method. Nucl. Eng. Des. 7 (1968) 170-178 3. Beck, J. V.; Litkouhi, B.; Clair Jr, C. R. St.: Efficient sequential solution of the nonlinear inverse heat conduction problem. Num. Heat Transfer 5 (1982) 275 ~286 4. Osman, A. M.; Beck, J. V.: Nonlinear inverse problem for the estimation of time-and-space dependent heat transfer coefficients. AIAA paper 87-0150 (1987) 5. Bass, B. R.: Applications of the finite element method to the nonlinear inverse heat conduction problem using Beck's second method. J. Eng. Ind. 102 (1980) 168-176 6. Mehta, R. C.: Estimation of heating rate using calorimetric probe. Rev. Sci. Instrum. 52 (1981) 1782-1784 7. Beck, J. V.: Nonlinear estimation applied to the nonlinear inverse heat conduction problem. Int. J. Heat Mass Transfer 13 (1970) 703-716 8. Mehta, R. C.: Solution of the inverse conduction problem. AIAAJ. 9, (1977) 1355-1356 9. Zienkiewicz, O. C.; Morgan, K.: Finite elements and approximations. New York: Wiley 1982 10. Owen, D. R. J.; Damjanic, F.: The stability of numerical time integration technique for transient thermal problems with special reference to reduced integration effects. Nmnerical Methods in Thermal Problems, Proc. Second Int. Conf. held in Venice, Italy, (1981) 487-505 11. Beck, J. V.; Blackwell, B.; Clair, S. R. St.: Inverse Heat Conduction, Ill-Posed Problems. New York: Wiley 1985 12. Mehta, R. C.: Estimation of heat-transfer coefficient in a rocket nozzle. AIAAJ. 19 (1981) 1085-1086 R. C. Mehta T. Jayachandran Vikram Sarabhai Space Centre Trivandrum 695022, India Received September 18, 1989