I'HEBIS v— ;_.. -vii' "L_i I WW 3‘! .....EL. :35} ;& that, This is to certify that the thesis entitled €71 fl/yeéc/{AEU ( G + U/LJ-éj’r‘y )1: $7.“ Job gru‘wg ”@443‘11/5, presented by \A’é’l ‘ C/Iit YOH} has been accepted towards fulfillment of the requirements for "(”52 c1", lather/degree in Hrchmaf' 1737‘ 144% Major professor Date V [0 - 5’6 0-7639 MS U is an Waive Action/Equal Opportunity Institution ~‘- th— _._..__a_ RETURNING MATERIALS: IV1531_] Place in book drop to remove this checkout from w your record. FINES wiH be charged if book is returned after the date stamped below. An AFFLICATION 0F UNSRADY SURFACE EMT NE'IBOD '10 13mm. cavm Roam BY WEI-CHI YANG A THESIS Sbni tted to Michigan State University in partial fullfillnent of the requirements for the degree of MASTER OF SCIENCE Department of lechanical Engineering 1985 ABSTRACT An APPLICATION OF UNSTEADY SURFACE ELEMENT METHOD TO mm MAL CAV ITY PROBLEM BY W El- (3111 YANG The unsteady surface element method (USEM) is a new numerical. systematic method for solution of transient. linear heat conduction problems. The method is deveIOped from Duhamel's integral and the principle of superposition. It involves a set of Volterra integral equations. one for each element. From the development of USER. another method called the quasi-steady state method is obtained and is an unique contribution presented herein. It involves a set of linear algebraic equation. one from the energy conservation principle, others from interface condi- tions between connected geometries. The quasi-steady state thermal behavior can be obtained without the knowledge of the transient thermal behavior. To show the advantages of the USEM, it is applied to the thermal cavity problem which is described in Chapter 2. In this thesis. one-. three-, and five-node solutions of USEM are used to analyze the problem instead of several hundred elements required by finite difference. It is found that very high accuracy is attainable with a relatively small number of surface elements. ii TABLE OF CONTENTS LISTOF TBLE OOOOOOOOOOOOOO0....0....OOOOOOOOOOOOOOOOOOOOOO LISTOF FIGURE OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO LISTOF symms OOOOOCOOOOOCOOOOIII0.0...OOOIOOOOOOOIOOOOOOOO 1. INTRODUCTION OF USEM .................................. 1.1 INTRODUCTION ...................................... 1.1.1 What is USEM ................................ 1.1.2 Characteristics of USEM ..................... 1.1.3 Motivation .................................. 1.2 Discussion of USEM ................................ 1.2.1 Duhmmel's Theorem ........................... 1.2.2 How to Apply Duhamel’s Theorem to USEM ............................. 1.3 Literature Review ................................. 1.4 Outline of the Thesis ............................. INTRODUCTION OF PR0 (1-2b) 6(r,t) = 0 for t=0 in region R (1-2c) Then from Duhamel's theorem T(r,t) can be expressed in terms of 9(r.t) [6] as t T(r,t) = 51:0fi(1)(39(r.t-T)/at)dt (1'3) where fi(t) is the forcing function in (l-lb). If there is more than one time-dependent boundary condition. the principle of superposition can be used. The number of separate prob- lems is equal to the number of time-dependent boundary conditions. If there are n time-dependent boundary conditions. the temperature can be expressed as n t T(r.t) = zi=1jt:0fi(t)(061(r.t-t)/Ot)dt (1-4) For each 9i(r.t) there are a corresponding field equation. boundary and initial conditions (1-2) to be satisfied. The influence function of the corresponding heat conduction prob- lem is denoted as 9- 1f fi represents temperature. 9 is called the temperature based influence function. If fi represents hcgt fluxes. 9 is called the heat-flux based influence function. For many geometries the influence function are known or can be obtained simply by analyti- cal or numerical procedures. In this thesis the influence functions are all based on the heat flux. 1.2.2 How to apply Duhamel's theorem in USEM In applying the USER, the interface condition between two connect- ed geometries is used to combine the temperature integral equations for the point on the interface derived from Duhamel's theorem and the prin- ciple of superposition. In Figure 1.2. the thermal cavity problem. point 1 belongs to body 1 and body 2. Body 1 is a solid cylinder with radius R and height E. Body 2 is a semi-infinite body outside a cyl- indrical void with radius R. Body 3 is a semi-infinite cylindrical void with radius R. qlm Body 2 Figure 1.2 Possible discretization fot the thermal cavity problem For body 1 as shown in Figure 1.3. the temperature at point 1 can be expressed by Duhamel's theorem and the principle of superposition T1 t 1 t t 1 = quoaelo/atdl - Ibqlaeil/atdl - IbqZBOizlatdl (1-5) 'hOtO q0' q1 and q2 are heat flux histories and are assumed uniform in space. but may be function of time; qo(t) is known value. For body 2 as shown in Figure 1.4. the temperature at point 1 can be expressed as t t T§ = Ibq066%o/6td1 + Ibqlae§llatdi (1.6) where the radial surface of body 3 is assumed to be insulated. efi is the temperature rise at point j of body k due to unit heat flux qi and is called the heat-flux based influence function. Q0 1 1 q0 l) - ' a ‘ 2 F"‘ ____ 2, k4 ‘ r r—" 1 ‘____ . ___.. 1 i p l 91(t) 91(t) h r r——' q (t) *—*“ ———~ 1 I *—_V ‘ I t q2(t) I Figure 1.3 Geometry of body 1. Figure 1.4 Geometry of body 2. From the interface condition. equations (1.5) and (1.6) are the same because there is no physical boundary between body 1 and body 2. Then combine (1.5) and (1.6) to get ‘ 12 t 12 t1 quOia(910-9107/3t1dl=1091[6(911+911)/at]dx+I0q23912/3tdl (1.7) where in some simple cases influence functions can be found in [7] and [8]. For more complex boundary conditions Green's functions can be used to calculate the influence function. The use of Green's functions in the procedure of calculating influence functions are given in the Appendices. Equation (1.7) is an integral equation for heat flow for the thermal cavity problem. The heat fluxes q1(t) between bodies 1, 2 and q2(t) between bodies 1 and 3 are uniform in space. but are func- tions of time. Equation (1-7) is a one-node integral equation of heat flow. The left hand side of (1-7) is known and the right hand side contains q1(t) and q2(t). two unknown functions. If the thermocouple is insulated, qZ(t) is equal to 0 and (1.7) is an integral equation with one unknown function. 91(t). This case will be dicussed in chapters 4 and 5. If the thermocouple is not insulated, another inter- face condition between body 1 and body 3 is needed. By using a similiar procedure, the temperature at point 2 can be expressed in two different forms. After the consideration of interface condition between body 1 and body 3. there are two equations with two unknowns (ie. q1(t) and q2(t) ). The influence of the thermal properties of thermocouple will be considered in Chapter 5. Equation (1.7) can be simplifed by Laplace transform or Stolz's method. Both methods will be used in this thesis. Laplace transform will be introduced in Chapter 4 and Stolz's method will be given in Chapter 5. 1.3 Literature Review N. R. Keltner and J. V. Beck [1] published the first paper discussing the USER. The authors used USEM to solve the intrinsic thermocouple problem which is a semi-infinite cylinder (wire) attached perpendicular to the semi-infinite solid (subtrate) as shown in Figure 1.5. In this paper the transient thermal behavior of the interface between the thermocouple wire and the subtrate due to a unit tempera- ture change in the subtrate was investigated by using the one-node USER. This resulted in a single integral equation to solve. The authors used both the temperature based and heat-flux based USER to get the integral equation and solved it analytically by utilizing Laplace transforms. A comparison was made with existing analytical and numeri- cal solutions given by other investigators. The results show excellent agreement between the one-node USER and other solutions. An advantage of USER is that only one node was needed in the interface and the results can be easily obtained without the need of a computer. I. V. Beck and N. R. Keltner [2] discussed the problem of two semi-infinite bodies contacting over a circular area by using one-node USEM as shown in Figure 1.5. Both the temperature and heat-flux based USER were employed. In this paper the perfect as well as imperfect contacting conditions were considered. The authors used the Laplace transform to simplify the integral equation derived from USEM. 10 Figure 1.5 The geometry for the intrinsic thermocouple problem Figure 1.6 The geometry of thermal contact conductance problem between two semi-infinite body. 11 B. Lithouhi and J. V. Beck [3] published the first paper dis- cussing the mulitnode USEM for two arbitrary bodies contacting over part of their surface boundaries. To demonstrate the capability of this method. it was applied to the problem of two semi-infinite bodies initially at two different temperatures suddenly brought into perfect contact over a small circular region. The temperature at the interface was expressed in the form of an integral equation with space and time variables by utilizing Duhamel's theorem and the principle of superpo- sition. In order to solve numerically the integral equation given by USEM, both the space and time domains was subdivided into small ele- ments. Then. the integral equations were transformed into a set of algebraic equations by Stolz's method of constant heat flux in each time step. The heat flux and temperature for each element of the interface can be obtained step by step in time. B. Lithouhi and J. V. Deck [4] discussed the same problem as [1] by using the multinode USEL. In multinode U313»; the single numeri- cal solution contains the whole time domain. 1.4 Outline of the Thesis Chapter 2 describes the thermal cavity problem and gives the mathematical model and literature review. Both the integral equations derived by USER and algebraic equations derived by the quasi-steady state method for the thermal cavity problem are given in Chapter 3. The procedure and results of the one-node USER using Laplace tranSfOIUS 12 is included in chapter 4. Chapter 5 gives the procedure and results of one-node USEM with Stolz's method. The influence of thermal properties of thermocouple on the temperature measurements is considered in this chapter. Chapter 6 describes the procedure and results of three-node USEM with Stolz's method. The five-node quasi-steady state solution is given in Chapter 7. Chapter 8 gives the conclusions from the results of thermal cavity problem by the USEM. Appendices include the deriva- tion of influence functions needed in the thermal cavity problem. CHAPTER 2 INTRODUCTION OF THE PROBLEM 2.1 Introduction In many dynamic heat transfer situations. temperature measurements are used to determine surface heat flux. In this thesis the tempera- ture is considered to be measured by a thermocouple which is located in a cylindrical hole normal to the heated surface as shown in Figure 1.1. The analysis is concerned with the case in which the thermocouple possesses much lower thermal prOperties than its surrounding material. Due to the existence of the thermocouple. the temperature distribution is disturbed. In turn the temperature measurements obtained from the thermocouple are different from those would exist if the thermocouple were not present. By modeling the heat transfer associated with the thermocouple, it is possible to deve10p corrections for the temperature measurements. One purpose of this thesis is to present an analysis for a thermocouple located as shown in Figure 1.1. 2.2 Motivation One motivation is to improve and extend the transient solution given by the finite difference (FD) method [5]. Because the body is semi-infinite. the number of FD elements needed to analyze the problem becomes greater as time increases. The dimensionless time (t*=ot/E3). 13 14 where a is the thermal diffusivity of the material and E is the dis- tance between the top of thermocouple and the surface of material, was used only up to 30 for the FD method. Since in the USEL the time domain of the transient solution can be any value. this method can be used to contain a larger time domain for the transient temperature his- tories. 2.3 Objectives 1) The first objective is finding the transient temperature at the t0p of the thermocouple as well as at the hot spot. which is at the heated surface right above the thermocouple. of the thermal cavity problem as shown in Figure 1.1. 2) The second objective is to determine the "null point" [5]. In order to obtain an accurate surface temperature history. the thermocou- ple has to be placed at some distance below the surface. The temperature difference between the top of thermocouple and undisturbed heated surface (ie. without the thermocouple) is the function of dimensionless time (t+=ct/Rz) and g(=R/E). So. one objective of this thesis is to determine the value of g _that makes the temperature difference between the tOp of thermocouple and undisturbed heated sur- face zero at time infinity. The location at which the temperature difference is zero at time infinity is call the "null point". It is the location at which the thermocouple should be placed in order to measure directly the surface temperature with the highest accuracy. 3) The third objective is to find the quasi-steady state tempera- in it 15 ture disturbance at the tap of thermocouple and the hot spot directly by the quasi-steady state method. and then to compare with the results extrapolated from the transient solution by USEM. It is a unique con- tribution presented herein. The detailed description and procedure will be given in Chapter 3. 4) The fourth objective is to find the needed influence functions of USER. These function are derived in the Appendices. 5) The fifth objective is to find the influence of the thermal preperties of thermocouple on the temperature measurements. It will be discussed in Chapter 5. 2.4 Mathematical Model A thermocouple embedded in the heat sink. which has a much higher thermal conductivity than the thermocouple. disturbs the surrounding temperature. The maximum temperature disturbance would be caused by replacing the thermocouple with an insulated cylindrical void with radius R and the distance E beneath the surface and perpendicular to the surface as shown in Figure 2.1. The heat conduction equation in cylindrical coordinates is: a’T/az3 + 1/r0(roT/dr)/Or = 1/adT/dt (2—1) The boundary and initial conditions for the insulated cylindrical ther- mal cavity problem are: 16 de(0.r.t)/az = -q.(t) (2-2) 3T(E.r.t)/az = 0 r\L :(:§;; /’/’./r Figure 3.4 Geometry representation of Ogi. If there are n elements on the radial surface of cylinder. only one half of the influence functions. eji' must be calculated from the expression given by (A-9). Others can be obtained from the following relationship. I _ 1 9ji " 9n+1—j.n+1-i (3‘5) 31 which means that the temperature rise at node j due to unit step heat flux Qi is equal to the temperature rise at node n+1-j due to unit step heat flux qn+l-i° For example. Gil = 9&3 for n=3 as shown in Figure 3-5. _. F‘— / N q1=1 ‘b— / \ .__, / \ —7 r— / \ / \ j \ / \ / Q ~ - L. I‘— j : ‘___~‘3 arc—— / ""— —"~ \ —"-." I ‘_ K ~— (13:1 A \\ —/ \ ._ Figure 3.5 6&1 = 9&3 for n=3 The symbol 93'n+1 denotes the temperature rise at node j of body 1 due to the unit step heat flux on the bottom surface of soild cylinder (ie. the n+1th element) and no heat flux over the other surface. These influence functions. 9}'3+1. can be obtained from the knowledge 0f 9}0 because qn+1 passes through the bottom surface of while qo passes through the tOp surface of the cylinder. relationshop can be used in the calculation of 9} n+l° the cylinder The following 32 1 1 °J.n+1 = °n-j+1.o (3-6) For example 9&0 - 0&4 for n=3 as shown in Figure 3.6. \\\\\ U) m /‘——‘\ ‘\\5\t‘\§f\\ z'./’,z' /",/' r /////////7 :>;*~ef~s;5<\:~efi~:>x_‘~ss‘\eu {3‘ \/ The temperature at node j on the body 2 due to qi(t), i=0. 1, ....n. can be expressed as I ah t I§(.) = qu0(l)[89§0(t-l)/8tldl + 21:1qui‘1)[39§i=satlkR+R<11-6s‘—123’c>/<24kg>+R’<3g’+2g’c-1)/4kg’(ant)"’ -2 = c1 + c2(t*)"/’ <4—24) 51 where 1* 21(t) is the dimensionless temperature difference between that at the top of thermocouple and that at undisturbed heated surface and c1 = (533+123’c-15)/(24g) (4-25.) C2 = -(3g'+2g’c-1)/(4g’n‘/‘) (4-25b) t+ = et/R’ (4-25c) Equation (4-24) is an important equation which shows that the transient temperature disturbance caused by thermocouple is function of 8 and (t+)-1/' for large times and the steady state disturbance caused by the thermocouple is function of g only. Figure 4.2 contains the temperature history. equation (4-24). at the top of thermocouple versus the value of g. Table 4.2 contains the value 0f C1 and C2 versus the value of g. C1 is the quasi-steady state temperature difference between that at the tap of thermal cavity and that at the heated surface far away from the thermocouple. If C1 is greater than zero. it means that the temperature measured by the ther- mocouple is greater than that at the undisturbed heated surface for 13:80 time. If C1 is less than zero. it means that the temperature measured by the thermocouple is less than that at the undistured sur- face for 18183 time. If C1 is equal zero. then the temperature measured by the thermocouple is the same as that at the undisturbed heated surface for large time. The value of g corresponding to a zero 52 value 0f C1 can be obtained by drawing the value of 01 versus the value of g. The location corresponding to this value of g is called the ”null point". which is the location that the thermocouple can measure the undisturbed heated surface temperature with minimum error. C2 represents the trensient disturbance due to the existence of thermocou- ple. From (4-24) transient disturbance decays as the function of (t+)-1/3. 1.4 l 1.2 \ 1.0 o“ :94 _m —°—- g}; //C/ -.4 / '1-00 0.5 1.0 2.0 3.0 4.0 g = R/E Figure 4.2 The transient temperature histories at the tap of the thermocouple. 53 Table 4.2 The steady and transient disturbance coefficient (i.e.. C1 and C2) at the top of thermocouple. g .1 .25 .5 1.0 1.333 4.0 C1 -6.086 -2.251 -.7579 .14 .3714 1.7283 C2 18.890 1.2712 -.2732 -.5727 -.5910 -.5391 The temperature at the "hot spot" is obtained by setting n=1 in equation (3-12). t T; = qoego - qu1(k)69§1(t—A)/atdx (4-26) Taking the Laplace transform on both sides of equation (4-26) and using (4-13) gives 7% = goégo - $31511 (4-27) Substituting equations (4-14) and (4-19) into (4-27) yields T; = qoléfio - 5] (4-28) Taking the inverse Laplace transform of (4-28) and substituting (4-21) into (4-28) gives T%*=I243~<3g‘+2g‘c—1>/4g’”‘ 54 =D.+n,(t+)“/’ (4-29) where T3+“) is the dimensionless temperature difference between that at the hot spot and that at the undisturbed heated surface and 01 (63’ + 1zg‘c - 3)/24g (4-30a) 02 .(331 + zg’c — 1)/4g’n‘/' (4-30b) Both the value Of 01 and D2 are in Table 4.3. The symbol D1 represents the quasi-steady state temperature difference between that at "hot spot" and that at heated surface far away from the thermocou- ple. Due to the disturbance of the thermocouple. the value of D1 should always be greater than zero for any value of g. This charac- teristic can be used to check the results of one node USEM by using the Laplace transform method. The symbol D2 represents the transient dis- tubance due to the thermocouple. The transient disturbance also decays as a function of (t+)_3/3. 55 Table 4.3 The steady and tranient disturbance coeficient (i.e.. 01 and Dz) at the hot spot. 3 .1 .25 .5 1.0 1.333 4.0 D1 -1.0858 -.1883 .2422 .6400 .8238 1.8533 DZ 12.8964 1.2712 -.2732 -.5727 -.5910 -.5391 4.6 Observations 1) From the above calculations. it is shown that the one-node USEM with the Laplace transform method is a very easy and quick procedure for the solution of the thermal cavity problem compared with the pro- cedure needed in the finite difference method [5]. 2) From Figure 4.2 the transient temperature histories by the one-node USEM are greater than those obtained by the finite difference method from at/E’= 0 to 30. But the difference decreases with the increment of g. It means that the one-node assumption is suitable for large values of g. 3) From the inspection of Figure 4-2. the temperature histories at the top of the thermocouple from one-node USEM with the Laplace transform. the best location for the thermocouple is at g=.925. which is the location of "null point". At this location the temperature measurements have the minimun error with respect to the undisturbed surface temperature. From reference [5] the calculated value of g is equal to 1. 4) From (4—25b) and (4-30b). the expressions for C2 and D2 are the 56 same. It means the temperature difference between that at the hot spot and that at the tap of thermocouple is a function of the value of g only. This phenomenon is true only for large times. 5) From equation (4-24) and (4-29). the transient disturbance caused by the thermocouple is function of t“/‘. This characteristic will be used in the following chapters to extrapolate the quasi-steady state temperature from the transient solution by USEM. Then these values can be compared with those given from quasi-steady state method to show the accuracy of these two methods. 5) From Thble 4’3: the V4130 0f ”1 for g=.1 and .25 is negative. It shows that the quasi-steady state temperature at the hot spot is less than that at the heated surface far away from the thermocouple. This situation can not happen in the physical world. So. it means that the one-node assumption is not sufficient to analyze the problem if the value of g is smaller than .5. For the same thermocouple assembly ( R=constant ). a small value of g means that the distance between the tap of the thermocouple and heated surface increases. The assumption of uniform heat flux in space over the interface between body 1 and body 2 is not correct. Three and five nodes on the interface between body 1 and body 2 will be considered in analyzing the situation of small values of g in Chapters 6 and 7. CHAPTER 5 ONE-NODE TRANSIENT'SOLUTION BY THE STOLZ'S METHOD 5.1 Introduction Equations (3-14) and (3-15). integral equations of heat flow derived from the USEM. represent a set of n+1 Volterra equations of the first kind with the unknown heat fluxes. qi(t)'., appearing inside the integrals. These integral equations can be approximated by a system of linear algebraic equations by replacing the integral with suitable qua- drature formulas. . The first step is to divide the time region 0 to t into n finite equally small time intervals. At. so that t represents the value of t at the end point of the mth interval. t = tn = mAt (Sfl) where At+ (ie. eAt/R’ ) in this thesis is taken to be 1. 2 and 3 in the calculational procedure. Then equations (3-17) and (3-18) can be rewritten as. ‘n m ‘k 1 2 " tk 1 [i=1 kglItk_1qi(X)a[9j.i+ej.illatdx+2k=ljtk_1qn+1(1)39j.n+1/atdx =‘m ]tk 1 e2 ._ 2k=l tk-1q0(l)a[ej.o- j.0]/atd3. J-lsZs eees B (5-2) Li=12t=1[tk,19i(1)39n+1.i/atdA + 2k=1jtk_1qn+1(1)3[9n+1.n+1 + \ 'II t 92+1'n+1]/3td1 = 2k=1It:_1q.(1)66:+1.0/atdx (5’3) 57 58 where to = o (5-4) In the simplest form of approximation the heat flux histories qj(t) are assumed to be constant values in each time interval. which is Stolz's method. and is shown in Figure 5.1. 4 q(t) ,//// \ // / qn= q[(m-.5)At] ) / + t to ‘1 ‘2 t.-1 hm Figure 5.1 Geometry illustrating uniform heat assumption over each time interval. Then ‘1: tk ]tk_lqi(1)l30§,i(t-x>/8tldx qi(k)Itk_1[30§’i(t-A)/dt]dh qi(k)[9§.i‘tm’tk-a)'95.i(tm'tk)] = qikee§1.e-t (5-5) 59 where 91k = qi(k) = qi[(k-1/2)At] (5’5) 1‘ Aejim-k ‘ °§i.m-t+1 ’ 9§i.e-t (5-7) and eji.k is the temperature rise at node j of body r due to unit step heat flux qi at time tk=kAt. The initial temperature. 0§j.o. is 0. So. equations (5-2) and (5-3) can be simplified by using the assumption of constant heat flux over each time interval to obtain \3 m 1 2 m 1 zi=a2t=1qi.k[49j.i.e—k * Aej.i.m-k]+2t=19n+1,kAej.n+1.m-k ‘ m =2k=1q0klA93.0.m'k-Ae§.0.m-k] j=1' 2' .... n (5-8) \n ‘m A01 + ‘m [ 01 +Ae3 J [i=12k=lqik n+1.i.m-k 2k=1qn+1.k A n+1.n+1.m-k n+1.n+1.m-k .‘9 1 ’2t=190kA9n+1.0.e-k (5'9) From equations (5-8) and (5-9). the heat fluxes qim" ( for i=1.2. .... n+1) can be determined at different time intervals one after another by marching forward in time for. n=1. 2. 3. ... Thus for each time step. equations (5-8) and (5-9) represent a system of n+1 linear algebraic hquations with n+1 unknowns. q1.m' 92.m' .... qn+1'm. and the heat fluxes at previous times. qi.k for k=1. 2. '60 .... m-I. are all known. By using the same assumption of constant heat flux over each time interval the temperature integral equations (3-9) and (3-12) can be simplified as. 1 "n 1 n+1m 1 Ih+1("’ ‘ 2k=1q0kAen+1.0.m-k ‘ 2i=12k=1qikA9n+1.i.m-k (5‘10) n+2‘m) = 2k=1q0kA9n+1.0.m-k ' 2k=1 k=1qikA9n+2.i.m-k (5‘11) where (5-10) and (5-11) are the temperature histories at the top of thermocouple and the hot spot. respectively. for time t=tm, T(m) means that the temperature at time t=tm, Thus for each time step. equations (5-8) and (5-9) represent n+1 linear algebraic equations with n+1 unknowns qlm' 92m, ..., qn+1.m° After these n+1 equations are solved. these n+1 quantities. qim'g, .rg then substituted into equations (5-10) and (5-11) to get the tempera- ture at time tIn (for m=1. 2. 3. ...). The outline of this chapter is as follows: Section 5.2 includes the objectives of this chapter. Section 5.3 gives the problem descrip- tion. Section 5.4 lists the influence functions needed in the calculational procedure. Section 5.5 includes the procedure of tran- sient solution by Stolz's method. The quasi-steady state method will be presented in section 5.6. Section 5.7 gives the observations. 61 5.2 Objectives 1) The first objective is finding the influence of the thermal prOperties of thermocouple on the transient temperature history at the t0p of thermocouple by the one node USEM with Stolz's method of con- stant heat flux in each time step. 2) The second objective is finding the influence of the thermal properties of thermocouple on the location of null point. 3) The third objective is finding the quasi-steady state tempera- ture at the tap of thermocouple as well as at the hot spot directly by the quasi-steady state method with the assumption of insulated thermo- couple. 5.3 The Problem Description 1) The initial temperature distribution of the whole domain is assumed to be 0. 2) The combination of thermal properties of thermocouple. kpc. is assumed to be much lower than that of surrounding material. k is the thermal conductivity. p is the material density. and c is the thermal capacity. 3) The contact between the t0p of the thermocouple and the materi- al is perfect. 4) There is no heat flux passing through the radial surface of thermocouple. 5) The heat fluxes between bodies 1. 2 (ie. q1(t) ) and 1, 3 (ie. 92(t) ) are uniform in space. but function of time. So. the one node -62 USER is considered in this Chapter. 6) Point 1 is at r=R and z=E/2: Point 2 is at r=0 and z=E; Point 3 is at r=0 and x=0. 7) The situation is shown in Figure 5.2. //// // / '1 Figure 5.2 One-node USEM with the consideration of the thermal prOperties of thermocouple. 5.4 The Influence Functions. The following are the influence functions needed in the procedure of one node DSEM and quasi-steady state method. These influence func- tion. except (5-23) can be obtained from reference [7]. The symbol 93k means the temperature rise at node j of body i due to unit heat flux '63 0}o(t+)=R/k[gt+-1/24g-2/gn'2:=1(-l)nexp(-4n'n'g't+)l4n']0(5-12) e}1(t+) = R/kl 2t+ + .25 - 22:=oexp(-a;t+)/a; ] (5-13) 9}2(t*> = 0}0(t+) <5-14> 0%0(t+)=R/k[gt+-1/6g-2/gn’2:;1(-1)nexp(-4n’n‘g’t+)In’] (5-15) e§1(t*) = R/Kl 2t+ - .25 — 22:;lerp(-e;t+)/a31.(en) ] (5-16) e§2(t+)=R/k[gt++1lsg+2/gn’2n=1(-1>“erp(-4n’u’g’t*)In’] (5-17) 950(t+) = 922(t+’ (5-18) 9%1(t+) = 0%1(t+) (5-19) 952(t+) = 9%o(t+) (5-20) 932(t*) = 900(t+) = (zk/k)(t+/n)"'tT (5-21) e§o(t*)=R/k[2(t+/n)"’erp(—1/1eg‘t*)—erie(1/4gt*"’)I23] (5-22) 2 t . 911 = a/kjodtj‘eazocxzods (5-23) -64 where t+(=ut/Ra) is called the Fourier number and is a measure of the rate of heat conduction in comparison with the rate of heat storage in a given volume element. In equation: <5-13> and (5-16) on is satisfied by the relation of 11(un)=0 where 11 is the Bessel function of the first kind with order 1- JO is the Bessel function of the first kind with order 0. An important ratio of thermal properties is kT=(ktpt°t/k9°]1/a' c and p are the thermal capacity and density of the material. The sub- script t means the pr0perty of the thermocouple. In this chapter the values Of kT are taken to be 0. 0.01. 0.05 and .10. In equation (5-23) the influence function is derived from Green's functions. which is discussed in the Appendix A. In general. GRij meani that the Green's function in cylindrical form with the ith kind inner boundary condition and jth kind outer boundary condition. The 1st kind of boundary condition is that the temperature is prescribed along the boundary. The 2nd kind of boundary condition is that the normal derivative of temperature along the boundary is precribed. The 3rd kind of boundary condition is that the combination of temperature and its normal derivative is prescribed along the boundary. The symbol 0 means that there is no boundary. In this thermal cavity problem i=2 and j=0. The Green's function in rectangular coordinate is GXij with the ith kind of boundary condition on the left (or tap) boundary and jth kind of boundary condition on the right (or bottom) boundary. In this thermal cavity problem i=2 and j=0. In equation (5-23) ds'=2anz' and z' is from 0 to E. The value of '65 9&1 cab be obtained from Appendix B. equation (8-10) and (8-14) by set- ting n=1. In equations (5’13) and (5'15) “n is satisfied by the relation of 31(un)=0 where 11 is the Bessel function of the first kind with order 1. JO is the Bessel function of the first kind with order 0. An important ratio of thermal prOperties is kT=(ktptot/kpo)‘/’, o and p are the thermal capacity and density of the material. The sub- script t means the preperty of the thermocouple. In this chapter the 741005 0‘ k1 are taken to be 0. 0.01. 0.05. and 0.1 In equation (5-23) the influence function is derived from Green's functions. which is discussed in the Appendix A. In general. GRij means that the Green's function in cylindrical form with the ith kind inner boundary condition and jth kind outer boundary condition. The 1st kind of boundary condition is that the temperature is prescribed along the boundary. The 2nd kind of boundary condition is that the normal derivative of temperature along the boundary is precribed. The 3rd kind of boundary condition is that the combination of temperature and its normal derivative is prescribed along the boundary. The symbol 0 means that there is no boundary. In this thermal cavity problem i=2 and j=0. The Green's function in rectangular coordinate is GXij with the ith kind of boundary condition on the left (or top) boundary and jth kind of boundary condition on the right (or bottom) boundary. In this thermal cavity problem i=2 and j=O. In equation (5-23) ds'=2anz' and z' is from 0 to E. The value of 9&1 can be obtained from Appendix B. equations (8-10) and (8-14). by setting n=1. 66 5.5 Calculational Procedure In this chapter only one node is assumed in the interface between body 1 and body 2. So. n=1 in the heat flow algebraic equations. (5-8) and (5-9). Now there are two equations with two unknowns. q1(k) and 02(k). for each time step. These calculated heat fluxes then are sub- stituted into (5-10) and (5-11) to find the temperature histories at the tap of the thermocouple and the hot spot. The algebraic equation of heat flow from the interface condition between body 1 and body 2 is found by seeting n=1 in equation (5-8) m 1 2 1 2t=1Q1kl [A911.m-k + A911.m—k]+q2kA912.m-k ' qotlAelo.n-t ’ A0§0,m_k] 1 = 0 (5-24) The algebraic equation of heat flow from the interface condition between body 1 and body 3 is obtained by setting n=1 in equation (5-9) ‘ m 1 1 3 2k=ll glkAGZI’m_k + q2k(A922sm-k + ABZme—k) - quAe%0sm-k] =0 (5-25) Now there are two equations (5-24) and (5-25) with two unknowns. q1k '04 92k- 91k q1[(k-1/2)At] (5-263) 92: q2[(k-1/2)At] ' (5-26b) 67 Define A1 = 0&0 — 0&0 (5-27.) A3 = 0&2 (5-273) 81 = 0&0 (5-23.) 82 = 9&1 (5’28b) B = 1 + 3 - 3 922 922 (5 28c) Substituting the above expressions. (5-27) and (5-28). into equa- tions (5-24) and (5-25) and assuming that qo is constant in time as well as in space gives. .m qu1(m) = 2k=1l 91(k)AA2(m’k) + q2(k)AA3(m-k) ] (5'29) m qoBl(m) = Zkzll q1(k)ABz(m-k) + q2(k)AB3(m-k) ] (5'30) where n=132531e eO AA(m-k) = A(m-k+1) - A(m-k) (5-31a) AB(m-k) = B(m—k+1) - B(m-k) (5-31b) q1(k) = 91k q2(k) = 92k (5-31c) 2:=1AA(m-k) = AA(m-1) + AA(m-2) +...+ AA(0) [A(m)-A(m-1)] + [A(m—1)-A(m-2)] +...+ [A(1)-A(0)] A(m)-A(O) A(m) (5-32) '68 Define I“ 1 ' QOA1(m) - 2;: I q1(k)AA2(m-k) + q2(k)AA3(m-k) 1 = C(m) (5-33) m-l qoBl(m) - )k=1[ q1(k)ABz(m-k) + q2(k)AB3(m-k) 1 = D(m) 15-34) where C(m) and D(m) are known values at t=t Substituting (5-33) and (5-34) into equations (5-29) and (5-30) gives C(m) q1(m)A2(1) + qC2(m)A3(1) (5—35) D(m) q1(m)82(1) + q2(m)B3(1) (5-36) The left hand side of equations (5-35) and (5-36) are known. q1(m) and q2(m) are the unknowns to be calculated. They can be obtained step by step by marching forward in time for. m=lp2931ee These two unknowns. q1(m) and q2(m). can be obtained by the fol- lowing formulas. q1(n) A1/A (5-37t) qZ(m) Az/A (5‘37b) 69 where A2(1) A3(1) C(m) A3(l) A2(1) C(m) A 2 82(1) 33(1) .D(m) 83(1) 82(1) D(m) The temperature at the tap of thermocouple and the hot spot at tIn can be obtained by substituting n=1 and j=1 in equations (5-10) and (5-11) to get 1 -'.m 1 T2(m) - 2t=1l quA92O.m-k ’ qlkA921.m-k ‘ q2kA922.m-k 1 (5'33) 1 _ ‘9 1 1 1 T3(m’ - 2t=1l quA93O.m-k ’ qlkA931.m-k ‘ 92kA932.m-k 1 (5'39) where T%(n) and T%(m) are the temperature at the tap of thermocouple and the hot spot. respectively. at t=t”. The procedure of calculation is as follows: 1) Set m=1. 2) Substitute the value of m into equations (5-33) and (5-34) and calculate C(m) and D(m). 3) Substitute C(m) and D(m) into equations (5-35) and (5-36) to calculate the value 0f 91(m) and q2(m) by using equations (5-37a) and (5-37b). 4) Substitute 01(m) and q2(m) into equations (5-38) and (5-39) to calculate the temperature at the tap of thermocouple and at the hot spot. 5) Increase m by 1. then go to step 2. From the results of calculations. the transient temperature at the tap of thermocouple for large times can be expressed as T%+(t)= 7O (Ti-TbQ)I(QOR/k) = C1+Czt+-’/'. Tables 5.1 to 5.4 contain the values 7 Of C1 and C2 versus the values of g for kf=0. 0.01. 0.05. and 0.1. respectively. Figures 5.3 to 5.6 are the transient temperature histo- Z, ries at the tap of thermocouple versus the values of g for kT=°- 0,01, 0.05. and 0.1. respctively. Figure 5.7 contain the value of g cor- L, responding to the location of null point versus the values of kT, Table 5.1 C1 and C2 versus g Table 5.2 C1 and C2 versus g 7L 2. for kT=O. for kT=.01 g .50 .75 1.0 g .50 .75 1.0 C1 “.7474 “.1830 .1536 C1 “.8735 “.2982 .0384 C2 “.5475 “.6360 “.6440 _ C2 “.4470 “.5860 “.6180 Table 5.3 C1 and C2 versus g Table 5.4 C1 and C2 versus g a «Z for £;=.05 for kT=.1 g .75 1.0 1.5 g .75 1.0 1.0 C1 “.4480 “.1048 .3063 C1 “.5471 “.2112 .3063 C2 “.5040 “.5600 “.6050 C2 “.4290 “.5080 “.5830 '15 '72 .3 .2 .1 0 J) E -T°°/ //W 111/ eoR/k _ / Ox 0: 15 N N H X. N R f\\\ ' \ \\ ///// jig// ”1°12? .75 1.0 g = R/E Figure 5.4 Transient temperature at the tap of th ooooooo ple for k¥éo,01 \: qOR/k // // // / // -, ///// _.///// \I \ on N] 0‘ Us A u N H O H N W :5 110 ON g=R/E 2. Figure 5.5 Transient temperature at the tap of th ooooooo ple for kT=,05 Figure ///§ ///,//“* /./ // /%V/ //// // / // / / 2. 5.6 Transient temperature at the top of th ooooooo ple for kT=.1o 15 1.35 1.30 1.25 ——+ 1.20 /// 1.15 rm_-1 1.10 2/ . amt-me“ / -. 1 .95 ——- +~_—__nh- -WM__V, ”f .85 .75- 0 0.01 0.05 kT Figure 5.7 The null point location versus hf; 76 5.6 Quasi-steady State Solution In this section the quasi-steady state method derived in Chapter 3 is used to calculate some thermal behavior of interest at quasi-steady state. They are the temperatures at the hot spot and the tap of the thermocouple. For simple analysis the thermocouple is assumed to be insulated. The energy conservation equation is obtained by setting n=1 in equation (3-17) ql = RqO/ZE = gq0/2 (5“40) The large time temperature at node 1 of body 1 can be expressed by setting n=1 and j=1 in equation (3-18) T} = q09}0(t) — qle}1(t) + a (5-41) where the constant B is dependent on the boundary condition of body 2 and the interface condition between body 1 and body 2. The exponential terms in the influence functions 6&0 and 9&1 in (5-12) and (5-13) are drOpped because they are negligible for large times compared with other terms in the equations. Then substitute these modified functions into (5-41) to obtain I} = —q1R(2t+ + .25)/k + qOE(gt — 1/24g)/k + B (5-42) Substituting (5-40). the energy conservation equation. into (5“42). 77 then gives T} = q0R(“g/8 - 1/24g)/k + a (5-43) For g=.5 T} = qOR(-7/48)/k + a (5“44a) For g=1.0 T} = q0R(“1/6)/k + a (5“44b) For g=2.0 T} = q0R(-13/48)/k + a (5-44c) where equations (5-44) are meant that the temperature at node 1 of body 1 for large time is a function of g. the interface condition between body 1 and body 2 and boundary condition of body 2. The large time temperature at node 1 of body 2. can be expressed by setting n=1 and j=1 in equation (3-19) T? = q09%0(t) + q10§1(t) (5-45) Substituting (5-22). the influence function 9&0, and the steady state value of 9&1. which is from [13]. into (5-45) and replacing q1(t) by 90(t) by the energy conservation equation gives For g=.5 T} = q1R(1.515174)/k — qu/ng + Tom = qOR(-.6212065)/k + Toe. (5“46a) 78 For g=1.0 I} - q1a(1.05779)/t - qu/2gk + To. = qok(.028596)/k + 1b. (5-46b) For g=2.0 1} = qu(.703691)/k - eon/2.x + Tba . q0R(.453691)/X + 1b, ~ (5-46c) where equations (5-46) is meant that the temperature at node 1 of body 2 is determined by the value of g. The values of 1.515174. 1.05779 and .703691 are from [13] and The is the temeprature at the surface of semi-infinite body due to unit step heat flux over the surface. The interface condition between body 1 and body 2 is considered perfect. T} = 1%, For - s=-5 B = “qOR(.475373)/k + 1b, (5“47a) For 3:1.0 B = q0R(.195562)/k + Tb, (5-47b) For g=2.0 B = qu(.724524)/k + Toco (5“47c) The temperature at the top of thermocouple. T}, for 14183 times is Té = qu(gt - 1/6g)/k — q1R(2t — .25)/k + a (5-43) For g=.5 T%.= q0R(-.746206)/k + Ibo (5—49.) 79 For g=1. T} . qu(.153894)/k + 1b, (5-49b) For g=2. I} = q0R(.891191)/k + 1b” (5“49c) The temperature at hot spot. T}, for large times is I} = q0R(st + 1/3g)/K - q1R(2t — .25)/K + a (5-50) For s=o5 1% = q0R(.25379)/K e Ibo (5“51a) For g=l.0 Ti = qu(.653895)/K + Toon (5-51b) For 3:2.0 1} = q0R(1.141191)/K + The (5—51e) 5.7 Observations 1) The difference between temperature histories at the tap of the thermocouple by the one-node USEM. Figure 5.3. and the finite differ- ence method. Figure 2.5. with the assumption of insulated thermocouple decreases with increment of the value of g. The difference increases with the increment of t+. but the rate of difference approaches zero at large time. 2) Compareing the results by the finite difference method and one-node USEM with the assumption of insulated thermocouple. the value of g corresponding to the null point location is 1 by the finite difference method but it is only .85 by the one-node USEM. This is due 80 to the assumption of uniform heat flux in space passing through the interface between body 1 and body 2. Since this assumption increases the disturbance caused by the thermocouple. the thermocouple has to be placed a little far away from the heated surface. 3) From the Tables 5.1 to 5.4. the transient disturbance caused by the thermocouple (ie. the absolute value of 02 ) for the same value of g decreases with the increment of kT, The steady state disturbace caused by the thermocouple (ie. C1) also decreases with increment of kT for the same value of g. It means that the distubance caused by the thermocouple decreases with the increment of thermal proproties of the thermocouple. 4) From Figure 5.7. the value of g corresponding to the location of null point is increasing with the increment of kT- It means that the thermocouple should be placed near the heated surface with a larger value 0f k1 because the disturbance caused by the thermocouple decreases. Though one-node assumption is not accurate for the insulat- ed' thermocouple case. the error decreases as the kT value increases because more energy passes through the interface between body 1 and body 3. In this case one node in enough to model the interface thermal behavior between body 1 and body 2. 5) Comparing the dimensionless temperature at the tap of thermo— couple. 9%+=(T%-Tb¢)l(qu/k). the results from the quasi-steady state method and the transient solution of USEM are almost the same. It means that the both method are acceptible. CHAPTER 6 THREE-NUDE TRANSIENT SOLUTTONS OF USEM 6.1 Introduction In this chapter the interface between body 1 and body 2 is divided into 3 equally-spaced surface elements and the thermocouple is assumed to be insulated. Stolz's method of constant heat flux assumption in each time interval is used to simplify the integral equation of heat flow derived from USEM. Due to the use of 3 elements on the interface between body 1 and body 2. the results of the calculated thermal beha- vior is more accurate compared with those by the one-node USEM. especially for the small value of g. The outline of this chapter is as follows: Section 6.2 gives the objectives of this chapter. Section 6.3 describes the pgoblem. Section 6.4 includes the influence functions needed in the calculation- al procedure. Section 6.5 includes the calculational procedure. Section 6.6 discusses the quasi-steady state method. Section 6.7 includes the observations. 6.2 Objectives 1) The first objective is finding the transient temperature at the top of thermocouple by the three-node USEM with Stolz's assumption of constant heat flux over each time interval. 81 82 2) The second objective is to make an accurate prediction of null point location. 3) The third objective is to illustrate the quasi-steady state method with three nodes on the interface between body 1 and body 2. 6.3 The Problem Description 01(5) -——« \\\\\\\\ 1; 1.4 l 4, Figure 6.1 Geometry representation of 3-node USEM The following is a description of the problem. See figure 6.1 1) The thermocouple is assumed insulated. 2) The boundary between body 1 and body 2 is discretized into 3 83 equally-spaced surface elements. Each of which has radius R and height E/3. 3) Node 1 is at r=R and z=E/6. Node 2 is at r=R and z=E/2. Node 3 is at r=R and z=5E/6. Point 4. which is at the top of the thermocou- ple. is at r=0 and z=E. Point 5. which is the location of the hot spot. is at r=0 and z=0. 4) The transient temperature histories at the top of thermocouple is calculated by the three-node USEM. The quasi-steady state method will be used to find the quasi-steady state disturbance at the top of thermocouple as well as at the hot spot. 5) The heat fluxes 41(t). q2(t) and q3(t) are assumed uniform over each element in space and functions of time only. Stolz's method of constant heat flux over each time interval is used to simplify the cal- culational procedure. 6.3 The Influence Functions The following influence functions are needed in the calculational procedure. The temperature at node (or point) j of body m due to unit step heat flux qi is denoted as 0?}. 9‘}; = Ogi/(qiR/k) (6-1) The temperature at node j of body 1 due to a unit step heat flux qi which is over the radial region from a to b. can be obtained from Appendix A. equation (A-9) with n=3. and is denoted as 6}}, a4 91* - 2 */3 + c1 - A - - A (6-2 31 ' t ji 1 A2 3 ) where i. j = 1. 2. and 3 t+ = et/R’ (6-3) ' a: a s A} = 2 12exp(“Bm)/3Bm (6“4a) IF- “, _ A2 = 2 18cos(nnz+)sin(nn(a++b+)/2)sin(nn(b+-a+)/2) n: exp(“n'n'g't+)/(n'n’g') (6-4b) .c» ten A3 = 2 Z 8cos(nnz+)sin(nn(b++a+)/2)sin(nn(b+-a+)l2) n=1 m=1 erpt-(a;+n‘n‘g'>t*)Inn(a;+n’n‘g‘) (6-4c) 2+ = z/E z = (j-1/2)E/3 (6-5a) + _ a = a/E a - zi_1 (6“5b) h+ = b/E b = z. (6-5c) fin is satisfied by 11(Bm)=0. 11 is the Bessel function of first kind with order zero. Tables 6.1. 6.2 and 6.3 contain the constant c}i's for g=.5. 1.0. and 1.5. 85 Table 6.1 Values of C}i for g=.5 i 911 912 C13 921 C22 0.646841 “0.040476 “0.356447 “0.040469 0.330934 Table 6.2 Values of C}i for 3:1.0 1 1 1 1 1 C11 C12 C13 C21 C22 0.299399 0.027527 -0.076945 0.027534 0.194940 Table 6.3 Values of Cgi for g=l.5 1 1 1 1 1 C11 C12 C13 C21 C22 0.217573 0.047193 —0.014775 0.047197 0.155609 The principle of superposition can be used to check the accuracy 0f C}i"° For example. the temperature at node 2. z+=.5. due to unit step heat flux over the whole radial surface of body 1. 0}+, can be obtained from equation (6-2) as. 1+ _ 1+ 1+ 1+ 92 “ 921 + 922 + 923 (6’6) For the case of g=.50. the constant terms on the right hand side as of equation (6-6) can be obtained from Table 6.1 and the sum of these constants is .646841 “ .040469 “ .356447 = .249925. The theoretical value of constant term on the left hand side [7]. page 203. is .25. Other constant C}i's can be checked by using a similiar procedure. The temperature at point 4 due to a unit step heat flux qi can be found in Appendix C. equation (C-4) with n=3. and is denoted as 0}}, 0}: = 2t+/3 + Ci} + Bl + 82 + 83 (6-7) where .I on + Bl = 2F12exp(“B;t )IBB;J0(fim) (6-8a) .cm B2 = 42 ll sin(nnb+)“sin(nna+) ]exp(“n’n'g't+)/(n'n'gz) n: (6“8b) «>- on . + . + a a s a + B3 = 42 12 [ s1n(nnb )“s1n(nna ) ]exp[“(8m+n n g )t ] m= n=1 Inn(B;+n'n'g’)Jo(pm) (6-8c) where am is satisfied by J1(Bm)=0. The value of Ci} is in Tables 6.4. 6.5 and 6.6. 87 Table 6.4 Values of C}i for g=.5 1 1 1 C41 c42 C43 “.089575 -.033372 -.077055 Table 6.5 Values of C4i for 3:1,0 1 1 1 c41 c42 C43 -.121357 —.034432 -.043712 Table 6.6 Values of C31 for g=l.5 ' 1 1 1 C41 C42 C43 -.397066 -.121797 -.268917 The principle of superposition can be used to test the accuracy of C:i's. Fer example. the temperature at point 4 due to a unit step heat flux over the whole radial surface of body 1. 0}+; can be expresged as. 91* = H + 923’ + 923 (H) The constant terms on the right hand side of equation (6-9) for the case of g=.50 can be found in Table 6.4 and the sum is .268917 - .121797 “ .397066 = “.249946. The theoretical value of the constant term [7]. page 112. on the left hand side of equation is “.25. The temperature rise at node j of body 2 due to unit heat flux qh i=1. 2. 3. can be obtained from Appendix B. equation (8—10). and is denoted as 0%}, The steady state Pltt 0f 9%} is denoted as C3}. Tables 6.7. 6.8 and 6.9 contain the value 0f CE} for g=.5. 1.0 and 1.5. respectively. From these Table cg} = 0%}. c}3 = c3} and c§3 = cgz. Table 6.7 Table 6.8 Table 6.9 The value 0f Cg} The value of C3} The value of C3} for g=.5 for g=1.0 for g=l.5 C2 841121 02 543134 c2 413439 11 ' 11 o 11 . C12 -445051 Gig .324921 cfz .261967 0&3 .301813 ‘ c}3 .236434 c}3 .197672 C51 -444946 0%1 .324873 c§1 .261976 6%; .694849 ‘ 0&2 .454707 cfiz .349355 C33 -374903 c§3 .278038 c§3 .226652 031 .301315 c3} .236435 cg} .197672 0%; .374967 €32 .278038 c§2 .225579 c}3 .651746 c§3 .424298 c§3 .325793 89 The steady state dimensionless temperature at node 2. z+-1.5. due to a unit step heat flux over the boundary etween body 1 and body 2 for g=0.5 and 1.0 are [8] at g =0.5 c} = 1.515174 (6-10a) at g =1.0 C MN = 1.05779 (6“10b) From the principle of superposition the temperature at node 2 due to a unit step heat flux can be expressed as. C3 = 0%} + cg; + c§3 (6-11) where Cgi's can be found in Tables 6.7 and 6.8 for the case of g=0.5 and 1.0. at g=0.5 CE — .444946 + .694349 + .374903 1.514698 (6“12a) at n=1-0 Cg .324873 + .454707 + .278038 1.057618 (6-12b) The difference between (6“10) and (6-12) are very small.r So. this means that values in Tables 6.7. 6.8. and 6.9 are very accuracy. The temperature at node j of body 2 due to a unit step heat flux qo can be found in reference [7]. page 75. 90 036 - 2(t+/n)‘/‘exp(“z3/4g't+) “ z+erfc(z+/2gt‘/’ (6“13) The temperature at node (point) 1 of body 1 due to a unit step heat flux qO can be found in reference [7]. page 112. 0}+ = gt+ - (z+/2 “ 1/6) - 2/n’2 (“1)ncos(nnz+)exp(n’n’g't+)/n' (6-14) 6.4 Calculational Procedure In this section the 3-node USEM is used to calculate the transient temperature history at the tap of the thermocouple. Stolz's method of constant heat flux over each time interval is used to simplify the integral eqautions of heat flow derived from the USEM. The resulted algebraic equations are (5-8). with n=3 and no qn+1 term. because of the assumption of insulated thermocouple. The procedure of calculation is similiar as that in section 5.4. The algebraic equation of heat flow is ‘ 3' In 1 2 _ In 1 2 2i=12k=1qikt Aeji.m“k+A9ji.m-k J “ 2k=1q0kl 4930.n-rAeji.n—k 1 j=1.2.3 (6-15) Now there are three equations with three unknowns. qlk' 92k and q3k. qik = qi[(k-1/2)At] (6'16!) k k k Aejim-k = eji.m-k+l ' eji.m—k (6-l6b) 91 Set Aj(k) = e}o (6-17c) ”J“’ = 913(k) + 913(k) (6-17d) for j=1. 2. and 3 Substituting (6-17) into (6“15) and assuming that qo is constant in time as well as in space gives ‘ m quj(m) = 2k;1q1(k)ABj(m-k+1)+q2(k)ACj(m-k+1)+q3(k)ADj(m-k+1) (6-18) where ABj(m“k+1) = Bj(m“k+1) - Bj(m-k) (5-193) ACj(m-k+1) = Cj(m-k+1) “ Cj(m-k) (6“19b) ADj(m-k+1) = Dj(m“k+1) “ Dj(m“k) (6“19c) Define m-l Lj(m) = quj(m) — 2k=ll q1(k)ABj(m-k+1) + q2(k)ACJ(m—k+1) + Q3(k)ADj(m“k+1) (6-20) 92 Substituing (6“20) into (6-18) gives 11“") = 41(m)Bj(l) + q2(m)Cj(1) + q3(°)Dj(1) (6-21) where at each time step the left hand side of (6“17). Lj(m)'s. are known. The heat flux q1(m). q2(m) and q3(m) can be obtained by the following algorithm ql(m) = AIIA (6“22!) Q2(m) = Az/A (6“22b) 43(m) = A3/A (6“220) where 81(1) C1(1) 01(1) L1(m) C1(1) D1‘1) A: 82(1) C2(1) 02(1) A1: L2(m) C2(1) 02(1) 33(1) C3(1) D3(1) L3(m) C3(m) D3(m) 81(1) L1(m) 01(1) 31(1) C1(1).L1(m) 42= 82(1) L2(m) 02(1) A3= 32(1) c2(1) L2(m) 83(1) L3(m) D3(1) 83(1) C3(m) L3(m) The temperature at node 4. the tap of the thermocouple. can be obtained from (5-10) by setting n=3 . fl ' 3 ~ m - 1 1 T4‘fl’ - ZkFJQOkAe40.m-k’>i=12k=lqikAe4i.m-k (5‘23) 93 The procedure of calcuations is as follow: 1) Set m=1; 2) Substitute m into (6-20) and calculate Lj(m)'s; 3) Obtain the heat fluxes qi(m)" by (6-22); 4) Substitute the results obtained in step 3 into equation (6-23) to find the temperature at the tap of thermocouple; 5) Increase m by 1. then go back to step 2. Figure 6.2 shows the transient temperature historys at the tap of thermocouple versus the value of g. From this figure the location of null point is found to be at g=.95. From the results calculated by the USEM the temperature history at the tap of thermocouple can be expressed as Tfi+=c1+cz(t+)-1/' for the large times where C1 and C2 are given in Table 6-10 for the case of g=.5. 1.0 and 1.5. Table 6-10 The quasi-steady state and transient disturbance coefficients at the top of thermocouple by 3-node USER g 1.5 1.0 0.5 C1 0.5340 0.0612 “1.0304 C2 “0.6550 “0.6125 “0.2387 95 6.5 Quasi-steady State Solution From the principle of energy conservation. (3“17). the quasi-steady state equation is given by 01(t) + q2(t) + q3(t) = 3gqo(t)/2 (6-24) where (6-24) is suitable only when the rate of the temperature change within body 1 is almost zero. The temperature at node j of body 1 due to go, q}, 92 and Q3 for large times can be expressed as 1 _ 1 ‘ 3 1 TJ- - q0(t)0j0(t) +Zi=1qi(t)0ji(t) + a (15—25) where B is a constant which is dependent on the boundary condition of body 2 and the interface condition between body 1 and body 2. Substituting the related influence functions into (6-25) and deleting the exponential terms with “t as parameter gives 1 3 1 T1 = [ 13q0/7Zg - } lcliqim ]R/k + a (6—26a) 1: T2 = [ “Q0/24g “ 2 . 1C2}qi(t) ]R/k + B (6—26b) ‘= T1-[—11/72—"31 - 3 - go 3 2. 1C3iqi(t) ]R/k + a (6 26¢) 1: 1. where C i's are in Talbes 6.1. 6.2 and 6.3. J 9.6 The temperature at node j. j=1. 2. and 3. of body 2 due to unit heat flux qi"- i=0. 1. 2 and 3. for large times can be expressed as T2- = (“91- + 2 3 92-- o(t) (6“27) J “0 10 i=1 1191 Substituting the related influence functions into (6-18) and deleting the smaller terms gives « 3 T1 = [ "90/68 + 21=1C%}qi(t) ]R/k + Ibo (6“28a) T2 ‘ [ [2 + ‘ 3 2 2 - ’90 s zflcmut) 1m + To. (cs-231.) T2 - I -5 I6 ' 3 2 3 ' ‘10 8 + >i=lc3iqi(t) ]R/k + Tbm (6“28c) where the value 0f Cg} are in Tables 6.7. 6.8 and 6.9 and Tb, is the temperature rise at the heated surface of semi-infinite body due to unit heat flux over its surface. From the perfect interface conditions between body 1 and body 2 one obtains I} = I} j=1,2,3 (6-29) Now. there are four equations; one is from energy conservation equation. (6-24) and the other three are from the interface conditions. (6-29). There aslo are four unkowns. which are q1, q2' Q3 and B. After solving these algebraic equations one obtains 97 If g=.5 B = “.583536q0R/k + 10w If g=1. B = .125900q0R/k + T0” If g=l.5 B = .454823qOR/k + 10w (6-30a) (6-30b) (6-30c) The dimensionless quasi-steady state temperature at the top of thermo- couple. '12::(1‘4-1‘0‘”) /(QOR/k) s is If g=.5. 1: = -1.031338 If g=1-0. Ti = .057749 If g=l.5. It = .526313 The dimensionless quasi-steady state temperature at I; = (TS-Tbm)/(qOR/k). is If g=.5. I: = .306696 If g=1.0. 1* .610165 .869435 11 g=l.5. I: (6-31a) (6-31b) (6“31c) hot spot. (6-31a) (6-31b) (6-31c) 98 6.7 Observations a) Comparing the transient temperature histories at the t0p of the thermocouple. the difference between the USEM and finite difference method decreases with the increment of the value of g. b) Comparing the quasi“steady state temperature at the top of ther- mocouple. the difference between the quasi-steady state method and that extrapolated from the transient temperature derived from USEM increases with the increment of the value of g. The relative error at g = 0.5 is about 0.1% and 1% at g = 1.5. c) The null point location is found to be at g=.95. CHAPTER 7 FIVE-NODE QUASI“STEADY STATE SOLUTION 7.1 Introduction In this chapter the interface between body 1 and body 2 is subdi- vided into 5 surface elements. each of which has radius R and height E/S. For simple analysis the thermocouple is assumed to be insulated. The quasi-steady state method will be used to obtained the temperature at the tap of thermocouple as well as at the hot spot. The outline of this chapter is as follows: Section 7.2 gives the objectives of this chapter. Section 7.3 describes the influence func- tions needed in the calculations. Section 7.4 includes the procedure of calculation of the quasi-steady state method. Section 7.5 gives observations. 7.2 Objectives 1) The first objective is checking the results from the three-node USEN for the case of g=0.5. 2) The second objective is to find an accurate quasi-steady state disturbance at the t0p of thermocouple and at the hot spot for small value of g. 99 100 7.3 Problem description K\s:>\\\\\\x if L )5"1 [9 1.] 1.1 ////// 4 I {new} 1.1 \ Figure 7.1 Geometry representation of 5-node assumption The following is a description of the problem. See figure 7.1. 1) The thermocouple is assumed to be insulated. 2) The boundary between body 1 and body 2 is subdivided into 5 equally-spaced surface elements. each of which has height E/5 and radius R. 3) The heat flux qi(t) passing through element i is uniform in space. but function of time. 4) The quasi-steady state method for the thermal cavity problem 101 derived in chapter 3 is used to obtain the steady state error caused by the thermocouple at the tap of thermocouple and the hot spot. 7.4 Influence functions The influence function 0? J} means the temperature rise at node j of body m due to a unit step heat flux qi, 9?; = 9?i/(qu/k) (7-1) The temperature at node j of body 1 due to a unit step heat flux qi can be obtained from Appendix A. equation (A-9) with n=5. as where a}; = 26/5 * C11 ' A1 ‘ A2 ‘ A3 (7-2) i.j=1. 2. 3. 4 and 5 O a: A1 = 2F12exp(“8;)/SB; (7-3.) m 42 = 2 8cos(nnz+)sin[nn(a+ +b+)]sin[nn(b+ -a+)] n=1 - exP('(B; + n'fl’g’)/(n'n'g') (7-3b) _ 0° '00 + + + +_ + A3 —2 2 8cos(nnz )sinlnn(a +b )/2]sin[nn(b a )/2] n=1 n=1 exp[“(B;+n'n’g')t+]/nn(fi;+n’n'g') (7“3c) ' 102 r+ . (j-1/2)/5 a+ = zi_1/E b+ = zi/E; zi=iEl5 The value 0f Cji" are given in Tables 7.1 and 7.2 for Table 7.1 The value od Cji for g=.50 0}} .556266 c}2 .173013 C}3 -.049160 c}4 “.182574 c}5 —.247643 c}1 .173018 c}2 .334102 c}3 .039623 c}4 -.114232 c}5 “.182574 c}1 -.049159 c}2 .039624 c}3 .269042 g=0.5 Table 7.2 ’ The value of Cg} for g=0.5 c}1 1.641059 c}2 .560032 c}3 -.223333 c}. —.736041 c}5 -.992068 c}1 .560032 c}2 .857702 c}3 .047440 c}. —.479377 0}, -.746044 cg} “.223228 cg. .047440 0&3 .601691 and 103 The temperature at node 6 of body 1 due to a unit step heat flux q} can be obtained from Appendix C. equation (C-4) with n=5. as 0}§‘= 2t*/5 + 0}} - a1 - 32 - 33 (7-4) where i=1. 2. 3. 4. and 5 on B1 = 2mF12exp(“B;t+)ISB;J0(fim) (7-5a, hm + + 8 1+ 0 ’ 3 B2 = 42 [ sin(nnb )“sin(nna ) ]exp(n’n g t )/(n n g ) (7-5b) n=1 ' °° ‘°° + + + B3 = 42 12 ll sin(nnb )“sin(nna ) ]expl “(B;+n'n’g')t ] m: n: /[(a;+n‘n’g‘)nn10 0.5. 2) From Figure 7.3. dimensionless quasi-steady error at the hot 3P0t (=( Th'Tbo )/(qu0/k) ). the results between the 3-node USEM and finite difference method are very close for the case of g ) 1. But the .109 difference becomes greater as g is smaller than 0.5. So. the results obtained from the finite difference are accurate only when g is greater than 1. Furthermore the number of elements used in this problem is greater than 200 by the finite difference method. Finite difference needs more elements to analyze the problem for g ( 0.5. It can be con“ cluded that the USEM is superior to the finite difference in this problem. 3) From Figure 7.2. dimensionless quasi-steady state error at the tap of thermocouple ( =( To — To... )/(Rq0/k) ). the difference between quasi-steady state method and the finite difference method are negligi- ble. It proves the superiority of quasi-steady state method to the finite difference method. 12;“. qu/k “.4 -08 “1.2 '1.6 '2.0 “2.4 =3.2 =3.6 110 -25 -5 1.0 1.5 g = RIE Figure 7.2 Quasi-steady state temperature at the tap of thermocouple 111 z=E/E Figure 7.3 Quasi-steady state temperature at the hot spot 1 Us W 112 3- ode l-node 5 “node /1‘;6;:'(Lapl ce transform) / \\ l0 0.5 1.0 1.5 2.0 g = R/E Figutr 7.4 Quasi-steady state temperature at the hot spot CHAPTER 8 CONCLUSIONS The multinode unsteady surface element method (MUSEM) and the quasi-steady state method for solution of the thermal cavity problem with linear boundary conditions and perfect contacting interface condi“ tions have been presented. Both methods are developed from Duhamel's integral theroem and the principle of superposition. MDSEM involves a set of Volterra integral equation. one for each surface element. The quasi-steady state method involves a set of linear algebraic equations. one from the energy conservation principle and others from the inter“ face conditions. In this thesis the radial surface of the thermocouple is assumed insulated and the thermal properties of the thermocouple are very small compared with those of surrounding material. The numbers of elements taken in this thesis are 1. 3 and 5 instead of several hundred that are required by the finite difference method. From the results of calculations. the USEM and quasi-steady state methods can be concluded to have the following advantages compared with finite difference. finite element and boundary integral method. 1) They are more suitable for problems having infinite or semi-infinite space domain. 2) They are approriate for problems having a finite contacting area between connected geometries. The number of elements needed to analyse the problem is fixed for USEM but must increase with time for 113 114 other related methods if the space domain of problem is infinite or semi-infinite and very large time periods are needed 3) They are suitable for problems if the number of location of interest is not large. The number of influence functions needed in the calculational procedure is pr0portional to the numbers of location of interest. In FD and FE methods the solution at all internal nodes is unavoidable generated. whether or not this information is needed. 4) They are suitable for problems having a small but important region. such as body 2 in the thermal cavity problem. 5) Extremely accurate solution can be obtained with very few ele- ments compared with other related methods. 6) Simple analytical results can be obtained by using small number of nodes with the Laplace transform method. The disadvantages of USEM are as follows: 1) A difficult part in the USEM is determining the influence func- tions in the integral equations of heat flow. If the derivation of the influence functions is restricted to use of Green's functions. it is not suitable for problems having irregular boundaries. But if the same problem is needed to be analyzed again and again. finite difference and finite element methods can be used to find the influence functions. Then this information can be used later to analyze the same or related problems. 2) It is not suitable for problems having temperature dependent thermal pr0perties. 3) Since USEM is develOped form Duhamel's integral. it is suitable for linear heat conduction equation. But the boundary conditions can 115 be linear or nonlinear. In particular. the radiative boundary condi- tion can be investigated. APPENDIX A 6134an EXPRESSION FOR 911 The radial surface of solid cylinder with radius R and height E is divided into n equally“spaced surface element. each of which has height E/n and radius R. and is shown in Figure A.1. The center point in each element is called 8 node. 9}} means the temperature rise at node j due to a unit step heat flux over surface element i and insulated over the rest of the area. It will be derived by using the Green's function. (i-l)aln \l l f: T Tux/“:2, Figure A-l Geometry representation of 9}i 116 117 For the linear, transient heat conduction problems the Green's function method [6]. [7]. [14] and [15] is a simplified, systematic procedure for obtaining the solution. Consider the following three-dimensional. nonhomogeneous boundary value problem of heat conduction. v2T(r.t) + g(r.t)/k = (l/a)6T(r.t)/dt in region R t)0 (A-la) kiaT/ani + hiT = fi(r,t) on boundary Si t)0 (A-lb) T(r.t) = F(r) in region R t=0 (A-lc) where k is the thermal conductivity. a is thermal diffusivity. g(r,t) is the internal heat generation, F(r) is the initial temperature distribution and fi(r.t) is the non-homogeneous term in the boundary conditions. The solution of heat conduction problem (A-l) in terms of the Green's function G(r,tlr',t) is given [6] as T(r.t) = JkG(r;tlr'.t).t=oF(r')dv' +(a/k)Jo [IRG(r.tlr'.t)g(r'.t)dv']dr t m w] OIZFIJSin-nlrn)[ro=rifi(r.t)/kid3i]dt mm 1 where the first term of (A-2) is due to the initial temperature distri- bution, the second term of (A-Z) is from the internal heat generation and the third term comes from the nonhomogeneous boundary conditions. 118 The Green's function in (A-Z) is satisfied by the following equa- tions with the corrseponding homOgeneous boundary conditions. vaG(r.t1r'.t) + (1/a)6(r-r')6(t‘t) = (l/a)aG/dt t>t (A-3a) kidG/ani + hiG = 0 on si t>t (A-3b) G(r.tlr'.t) = o if t 119 If there is no internal heat generator and initial temperature distribution is zero. The solution in terms of Green's function can be rewritten as t n m T(r.t) = a] 0 dtZFlJ‘SG(r.t|r'.t)iro=rifi(r.t)/kid8i (A-4) i Then eji' the temperature rise at node j due to unit heat flux 9i with no heat generation and zero initial temperature distribution, can be expressed as 1 t 9.11 = a/kJo dtL’G(r.t|r'.t).rn=Rds' (A-S) where the unit heat flux qi is over the radial surface of the solid cylinder from a=zi_1 to bzzi and r=R and ds'=2anz' Th0- and three- dimensional Green's functions can be found for many cases by simple multiplication of one-dimensional Green's function [14]. This is true for almost all of the boundary condition of rectan- gular system. provided the problem is linear, the body is homogeneous and the geometry is ”simple". By simple geometry is meant that any boundary must be located only by one coordinate such as x=0 or x=L. but not at x+y=0 ,for example. For the cylindrical coordinate system one can use the multiplication method for the two-dimensional problem only for r, z coordinate and 9, z coordinate with the restrction of almost 120 constant value of r. For the Spherical coordinate system no multipli- cation relation is possible for the one-dimensional Green's function to get two- or three-dimensional Green's function. The one-dimensional Green's function of rectangular. cylindrical and spherical coordinate system can be found in [16]. [17] and [18]. By the multiplication method. eji can be written as t b 1 — ' . eji " u/LIO dTJ'.Gfloszz2(2flRdZ ) (A-6) where G " [ 1 >00 (- ?z( - le) ]/ P2 (A-7) [:02 - + Lmlexp GLUE t T) \ n» O >4 [‘0 N I oo [ 1 + 2) cos(nnz/E)cos(nnz'/E)exp(-n3n2a(t-t)/Ez) ]/E ‘-m=1 (A-S) GROZ and 6X22 are from [16] and [17]. respectively. pm is satis- fied by the relation of 11(bm)=0. J1 is the Bessel function of the first kind with order one. In general. GRij means the Green function in cylindrical coordinate with the ith kind of inner boundary condition and jth kind of outer boundary condition. GXij means the Green's func— tion in rectangular coordinate with the ith kind of boundary condition on its left (or tap) boundary and jth kind of boundary condition on its right (or bottom) boundary. 121 Substituting (A-7) and (A-8) into (A-6) gives 9}: = e}i/(qiR/K) = cji + 2t+/n - A1 — A2 - A3 (A—9) where A _. )(I 2 3 +)/ 3 1 - 1.31:1 exp(-[imt [3m (A—IO) M II 2:418 cos(nnz+) sin[nn (b+-a+)/2] sin[nn(b++ 3+ ) /2] exp(-n1nzgzt+)/(n’n’gz) (A-ll) A3 = Z:=12:18cos(nnz+) sin[nn(b+-a+)/2] sin[nn(b++a+)/2] epr-(B;+nznzgz)t+]/[nn(p;+n3nzgz)] (A'IZ) t+ = at/Rz g = lI/E zj_1 = (j-1)E/n 2+ = z/E; z = zj_1+E/2n b+ = b/E; b = 21; 2i = iE/n a+ = a/E: 8 = Zi-1 122 n is the number of equally-spaced surface elements C}. is the constant to be found If a(t-t)/Rz (.07. GROZ can be approximated [16] by GROZ = I (nx>"/’ + .s + .75(i/n)‘/‘ + sixa + 21A”’/(32n‘/’) J /(2nR2) (A-13) where l=a(t-t)/R3 (A-14) An alternate form for 6X22 [17] is given by CC 6X22 = ZnFiexp(-(2nE+x-x')z/4a(t-t))+exp(-(2nE+x+x')z/4a(t-t))] /<4na(t-r>’/’) (A-lS) Comparing (A-15) and (A-8). different expressions for GXZZ- the t variable in the exponential terms is at the numerator for (A-8) while is the denominator for (A-lS). So. (A—lS) is suitable for the small time domain while (A‘B) is appropriate for the large time domain. In the procedure of calculating the influence function. 931- a(t-t)/R3) is smaller than 0.07. (A—13) and (A-lS) are used for GROZ if l (1.3.. and GXZZ' respectively. If k is greater than 0.07. (A-7) and (A-8) are used for GROZ and GX2Z' respectively. 931 can be rewritten as 1 _ t b . eji - ZNRCL/kLzod‘tszojfixzzdz (A—16) i23 where GRZO' equation (A-7) or (A-13). is not function of 2. For the case Of a(t"t)/R3<.07. I 6X22 can be calculated by substi- tuing equation (A-15) and using the following formulas. b j zezdz' ~ [ erf((z++b+)/th+1/3) - erf((z++a+)/2gt+1/z) + l err((2+-.+)/2gt+"’) — ert<b+ (A-17a) b ] GXszz' = I err((z++b*>/2gt+"’) - erf< - erf< oexp(-p’x)[p(J;(p>+y;(p>>1"du (8-2) where k=a(t-‘c)/Rz (8-3) The integral form 0f GRZO is not easy to manipulate. For simple analYSiS the time domain 0f GRZO is subdivided into two parts. which are small time domain and large time domain. For each time domain there is a simple expression that can approximate (8-2). In this case the time domain is divided at dimensionless time i=6. For 1(6. by the least square error method with the polynomial model. the best one is given as GR20=(2nu’)”[ (nL)-1/z — .5 + .41343411/z — .299877i + .1544831’l’ -.o45263i’ +.oos484i‘/’] (e-4) The values of (8-4) and (B-2) can be compared to check the accura- cy of the approximate expression of GRZO for the small time domain. For example. at 1:1 the value of (8-2) is 0.2926330 the value of (8-4) is 0.2924500 the relative error is .0625% at k=6 the value of (8-2) is .0680578 127 the value of (8-4) is .0682995 the relative error is .355% For the case of k>6. GRZO can be expressed [16] by Gazo = I 1 - L(1+3(1-L)/4i)/2x - (n’+4)/32i 1/(4nn’i) (D-S) where L=10g(4k)-7 7:.57722. which is called Euler number. 6X20 can be found [17] to be 5x20 = [ exp<(z-z')’/4a(t-c)) + exp<(z+z')‘/4a(t-r)) 1 /(4na(t-c))1/’ (8-6) Since GRZO is not function of 2. (8-1) can be rewritten as . _ b. . . Is szocxzods - ZHRGazoj‘zeodZ (L‘?) where Jbezodz' can be obtained by (A-l7) a . b Define J.GX20dz'=F (8-8) Substituting (8-8) into (8-7) and changing the time variable obtains .128 t+ 9%. = 2nR/kIO FGRzodA (8-9) The time domain 1 in (8-9) can be subdivided into two region. One is from 0 to 6; The other is from 6 to t+. (B-9) can be rwwritten as t+ GJ.i = 2nR/kJZFGR20di + ZnR/kIgFGRZOdh (8-10) For the case of A(6. only the first term of the right hand side (8-10) is needed to calculate. This term contains the integration of expressions of the combination of polynomial and complemetery error function. which can be solved by the following relations. b J.U(zn-’)lzerfc(aU-1/3)dU = 2t(3n-1)/1[ erfc(at-1/z)- “En(2n+1)(2n-1)... (2k-1)n n=1.2.3.... (8-12) where an ie the exponential integral function with order n En(z) = IZpr(-zt)/tndt (8-13) The value and formulas of (8-13) can be found in [19] 129 For the case of A)6. it is difficult to calculate the second term of the right hand side of (8-10) analytically. Since the value is very small compared with that of the first term. (8-10) can be approximated 5 m 9. . = ‘, , _ J1 an/kJ‘OI‘CROZdA "l’ Zh=12flR/R[FGR02]A1 (B 14) where FGROZ(h) = FGROZI6+(h-1/2)At] (8-15) APPENDIX C GENERAL EXPRESSION OF 91 . n+1.i 9h+1.i is the temperature rise at the n+1th node. the center point of the bottom surface. of the solid cylinder due to unit heat flux over the element i. which is from zi_1=a to zi=b' and is shown in Figure C-l. -../fl“ \~ \ . ‘ l l 1 E _fl~ I \\\ \\ W4 1 r I I///// NO x .1 // €91 ,WH%/”-/:> Firgue C-l Geometry representation of 9i+1 i From the multiplication method of Green's function 9h+1.i can be expressed as t 1 ' 9n+1.i=“/KJ =OdTISGR026x22dS (C-1) 130 131 where ds'=2anz'; o and k are the thermal conductivity and thermal dif- fusivity of the material. GROZ = [ 1+Z:1exp( a(t- tlbm/Ra )/Jo(pm) ]lnR (C-2) where 11(pm)=0. II is the Bessel function of the first kind with order one and JO is the Bessel function of the first kind with order zero. (I) 1 2 3 6X22 = [ 1+22n=1cos(nnz'/E)exp(n n a(t-1:)/E ) ]IE (c-3) Substituting (C-2) and (C-3) into (C-I). one obtains 1+ _ + 1 en+1. i - 2tIn + Cn+1.i — A1 - A2 - A3 (C-4) where n is the number of equally-spaced surface elements on the radial surface of the solid cylinder ”('1‘ A1=Zm=12exp(—p;t+) lnp;10(pm) (c—5) A2=4L:1[ sin(nnb/E)-sin(nna/E) ]exp(- -nzn 3 g3t+)/(n n 3’) (C-6) A3= 42n 12: [ sin(mnb/E)-sin(mna/E) ]exp[-(d;+mzn 3 gzlt+ ] = n=1 /(mn