. . v v to 13' ¢ , I L. f |I . . J I -l - I. t - . . , a h. .1 oil I... .. . if I!” . $307.. I .014....P.i tilt}, £0: ‘. . J 1)., 1 .0! V V 3.41.11. 7". . ll: it). all I . w )l 5 YO. I . {Willtuuflh 0 Mid. w- ¢.Ic i'. roof—{IHHH i l l’ I JIII'II' a ”l I." a. j I | o . t It . . | I lilt‘ o g: 4 lel. .34. I [or x . u. . I h I! .u . J T]. .clfllufllnfltut- I. C. .l Illllucnnfllnv if litaxl‘ilniv l‘vtlr I - . . . D. It! .1. £3! .lll' :1 la- .{Illln‘cfl «illo! . -l -IugluocD ' 2.. I! z \ - .v.-. . - Ln .1... - .. .o- 1". é“! . 1117‘! . DA \l u) \‘ . . o.rol.-nlllllln‘|m lg“! .‘I 001 .‘tfitfillnlli‘ «I‘lafllkv .95? II...) All... I 1 1.1.":4914"; [911 .l . u 1 LI}..- oi. Ilia-cellulrvlnajl. .. . .nl... $51!.Lnbflurmb . . Ion-sill: {vol}: .‘0‘. {Illllill (all!!! fi’l‘ut.; It"! 10': cg; 3"qi 10.50.. 4"!!. l‘.,..vuvu|IA-I1 \) (II-II .\I ll.‘ i o .4’- . ‘fl'ol: I‘l'ls 6|! u 1...»!ng lull! 1.31.413!” ”fill... fill; I..v «clv I”. “unliqullili‘ 4 luv. h‘lprlllill‘j‘.‘ a. 1.. It, \Ill I! “h! . v '1: w rt {Ila}. .fiigifit. . .a -.-|o 1 - . rol- .900....!u1)..u.xthtflh.l tr- fifkp}.}ll -. E J u an... . . I! - Ill}! c (.Illlv ~ ‘I \ uuuuu 1 o ‘0 l .1! ‘iofl’. ‘!"|. .1 c - -t . it II ' .PI- 1 i. ' II n '1 u - ’Io!‘ ‘9. ill: I «I 1 l 1.? I {Elli . u I a ., "11. v {n .‘Iu'lul... it. . v l‘yglww“. "IIan 117. I. {I} ‘gi. 1%1’51. ulllud.‘ . ‘llull' 0011.. \-1Ib)l.' I ll Lit.” :1]: . IL.--’ ’5'! \ ‘flvl ldflflblllll: .nvtli‘llel flI. II - HJ’I!Iv-|[|LJI Ill fih‘f .n 1} infill A Efl‘filfilllliréjl(~nl1.I - 'ul‘llltlr. fl I‘lr (III-V: fi if? 0' f]; ’.I.AI|$1D.. liirt!\-hll IF\“} I llIlvl'IlII. I'll! I‘l lllll 'IJ‘II “I ‘I’ I.» 1‘ o . ‘II II I. '1'," 1i % ( Ito 1' J . THESIS i‘fl";;‘7§ .ifif'q'.___‘ ‘ .L fl. 2‘; , a «Ls-"u ‘x’ \:u‘~‘w ;..--_.r ‘4 rim; .' 4""! V I" :5: Mat-=93] This is to certify that the dissertation entitled SURFACE ELEMENT METHOD IN TRANSIENT HEAT CONDUCTION PROBLEMS presented by BAHMAN L l TKOUH I has been accepted towards fulfillment of the requirements for PH.D. degreein MECHANICAL ENGINEERING 7% Major pvrofess [ \ Date 8-11-82 MS U is an Affirmative Action/Equal Opportunity Institution 0- 12771 tam ~ RETURNING MATERIAggz 1V153‘_] Fléce in book drOp to LIBRARIES remove this checkout from III-KSIIIL. your record. FINES will be charged if book is returned after the date stamped below. !___ SURFACE ELEMENT METHOD 1N TRANSIENT HEAT CONDUCTION PROBLEMS By Bahman Litkouhi A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 1982 ABSTRACT SURFACE ELEMENT METHOD IN TRANSIENT HEAT CONDUCTION PROBLEMS By Bahman Litkouhi The heat transfer between two (or more) bodies with perfect or imperfect contact at the interface is of fundamental importance in heat transfer studies and it has accordingly received considerable attention over the last two or three decades. It is important in the problems involving electric contact, electronic cooling, welding, fins, contact conductance and many other applications for which two (or more) similar or dissimilar bodies are attached one to the other over a small part of their surface boundaries. In general, it is dif- ficult to obtain analytical solutions to such problems. The transient surface element method (SEM) is a new numerical method for solution of linear transient two- and three-dimensional heat transfer problems. The method is well-suited for the above mentioned problems compared with the other numerical procedures such as finite difference (FD), finite element (FE) or boundary integral equation (BIE) methods. In the SEM only the interface between the two geometries requires discretization as Opposed to discretization of the whole do- main needed in the FD and FE or discretization of the whole boundary in the BIE method. This in turns reduces the size of numerical calculation and computer time. Bahman Litkouhi In this dissertation a multinode transient surface element method for two-dimensional heat conduction problems with linear boundary condi- tions has been deveTOped and formulated. The method uses Duhamel's in- tegral and involves the inversion of a set of Volterra integral equations, one for each surface element. Computer programs were written and the following three different problems were solved: i) two semi-infinite bodies initially at two different temperatures suddenly brought together over a small circular area and insulated elsewhere, ii) the intrinsic thermocouple problem, and iii) a semi-infinite body with mixed boundary condition if a step change of the surface temperature over an infinite strip. For each of these problems the multinode SEM perfomed well. The results showed excellent agreement with those obtained by other investi- gators. It was found that very high accuracy is attainable with a re- latively small number of surface elements. This feature makes the method superior to the alternative numerical procedures such as FD, FE, or BIE methods for the type of problems considered. AKNOWLEDGEMENT The author wishes to express his deepest appreciation and gratitude to his major advisor, Chairman of his Guidance Committee, Professor James V. Beck for his contributions, guidance and encouragement during the course of this research and also for his friendship and painstaking re- view of this dissertation. The author is also grateful to the other members of the Committee, Professor Mahlon C. Smith, Professor Larry J. Segerlind, and Professor K. Jayaraman for their guidance and valuable discussions. Financial support for this research was provided by the National Science Foundation under grant number CME 79-20l03 which is greatly ap- preciated. Gratitude is alSo extended to the Department of Mechanical Engineering and the Division of Engineering Research for their financial support and cooperations. To his wife Mahshad, the author dedicates this dissertation for her understanding and moral support during the graduate study and research. ii w LIST OF TABLES ............................................... LIST OF FIGURES .............................................. TABLE OF CONTENTS LIST OF SYMBOLS .............................................. CHAPTER I. CHAPTER 2. CHAPTER 3. INTRODUCTION AND DEVELOPMENT OF INTEGRAL EQUATIONS ........................................ Introduction ......... - ....................... Mathematical Description .................... Duhamel's Integral .......................... Derivation 0f Duhamel's Integral Equation for Arbitrary Time and Space Variable Heat Flux Boundary Condition in a Two-Dimensional Region ...................................... 1.5 Integral Equation Formulation ............... —l—l—l—l o o o o th—I SURFACE ELEMENT METHOD FOR TWO BODIES IN CONTACT ........ . ................................. Introductory Remarks ........................ Problems to be Considered ................... Discretization Over Space ................... Surface Element Formulation for Two-Bodies in Contact .................................. Solution of Simultaneous Integral Equations ................................... N NNNN O I O O C 01 Ade TNO SEMI- INFINITE BODIES IN CONTACT OVER A CIRCULAR AREA .... ................................ Introduction ................................ Statement of the Problem .................... Surface Element Solution .................... Temperature at the Surface of a Semi- Infinite Body with a Disk Heat Source ....... Thermal Constriction Resistance of the Two Solids .................................. (A) wwww O O O O 0 U1 «thd Page vii xi l2 T3 18 22 22 23 23 26 32 4o 40 43 44 46 52 Page 3.6 Results and Discussion ....... . ............... 53 CHAPTER 4. INTRINSIC THERMOCOUPLE PROBLEM ................... 75 4.1 Introduction ................................ 75 4.2 Statement of the Problem .................... 83 4. 3 Solution .................................... 83 4.4 Results and Discussion ...................... 86 CHAPTER 5. SEMI-INFINITE BODY WITH MIXED BOUNDARY CONDITIONS ....................................... 98 5. l Introduction ................................ 98 5.2 Statement of the Problem .................... 100 5.3 Solution ............... ..................... l07 5.4 Solution to the Interior Region ............. 107 5.5 Results and Discussion ...................... 107 CHAPTER 6. TEMPERATURES IN SEMI-INFINITE BODY HEATED BY CONSTANT FLUX q0 OVER HALF SPACE ................. 117 6.1 Introduction ................................ 117 6.2 Mathematical Description .................... 118 6.3 Derivation of the Solution .................. 118 6.4 Other cases of Boundary Condition ..... ...... 132 CHAPTER 7. SUMMARY AND CONCLUSIONS .......................... 144 APPENDIX A DERIVATION OF EQUATIONS (2.5.11a) and (2.5.11b) .. 148 APPENDIX B EVALUATION OF I(r+,b+) GIVEN BY (3.4.5b) ......... 150 APPENDIX C EVALUATION OF THE INTEGRAL IN EQUATION (6.3.11) .. 154 REFERENCES ................................................... 156 iv -4. . aha—r ' Table 3.1 LIST OF TABLES Page Values of the influence function, aij, for semi-infinite body at dimensionless time t; = .001. ...... 54 Values of the influence function, ¢ij, for semi-infinite body at dimensionless time tM = l. ......... 54 Values of the influence function, aij, for semi-infinite body at dimensionless time t; = 1000. ...... 55 Results for elemental heat flux qj, for an isothermal disk region on the surface of a semi-infinite body (t+ = at/a2). ......................... 57 Normalized spatial variation of the surface heat flux at different times for the geometry of Case a, see Fig. 3.1. (q+ = q/qCL vs r+ = r/a). ............... 58 Normalized area averaged interface heat flux for the geometry of Case a, see Fig. 3.1. 6+ = a'/ 5;. ...... 58 Results for elemental heat flux, qj, for a disk-Shaped interface of two semi-infinite bodies, one of C0pper and the other of glass (t+=a1t/a2). ........ 65 Results for dimensionless constriction resistance, for an isothermal disk region on the surface of a semi-infinite body R: = Rc-a-k. ....... ................... 70 Results for dimensionless constriction resistance for a disk-shaped interface of two semi-infinite bodies, one of copper and the other of glaSS’ 4. RC - RC/Rc(w). ........................................... 7l Normalized values of elemental interface temperature, TI, as a function of dimensionless J time t+ = azt/az. ...................... . ................. 87 V Table 4.2 5.1 6.1 6.2 6.3 6.4 6.5 Normalized area averaged temperature histories for chromel substrate and alumel wire. .... ....... Comparison between the results obtained from the local-¢ solution and the averaged-a solution. Values of Function H(X,p). ....................... Values of the functions H(X,p), and ERFC(pX) for various X values. ................................ Number of terms in series of function H(X,p) given by equations (6.3.20a) and (6.3.23a) to obtain 8 decimal places accuracy. ......................... Values of dimensionless temperature for the solution given by equation (6.3.18), for different values of X and Z. ..................... Values of dimensionless temperature for various values of x+ and 2+ due to a strip source at time t = .25. ........................................ vi ... ..... 128 ........ 128 LIST OF FIGURES Some geometries for intrinsic thermocouple problem .... Some geometries for contact conductance problem ....... Some geometries illustrating mixed boundary conditions. ........................................... Geometry showing a two-dimensional heat conduction problem. ..................... ......................... Geometries illustrating various types of boundary conditions in a two-dimensional region ................ Geometry showing a two-dimensional region heated by an arbitrary heat flux only over a portion of its boundary. ............................ . ................ Geometry showing a discretization over the heated portion of the surface boundary. ...................... Uniform heat flux assumption over each element. ....... Geometry of semi-infinite slab attached to a semi-infinite body. ........................ . .......... Possible distribution of surface elements for connected semi-infinite body and slab. ................ Basic building blocks for the geometries of semi-infinite body and semi-infinite slab. ............ Geometry illustrating uniform heat flux assumption over each time interval. ......... ................. .... Geometry of a semi-infinite body with a step change of the surface temperature over a circular area. ...... Distribution of surface elements for two semi-infinite bodies in contact over a circular area. . Semi—infinite body heated over an annular-shaped region aj_]l. ............................................ 125 6.3 Function H(X,p) versus X for different values of p. ... 129 6.4 Dimensionless temperature T(X,Z) versus X for different values of Z in semi-infinite body heated by a uniform heat flux over half space x<0, z=O. ...... 134 6.5 Lines of dimensionless isotherms in semi-infinite body heated by a uniform heat flux over half space x < 0. ................................................ 135 6.6 Various possible cases that can be treated using the solution given for figure 6.1 as a building . block. ................................................ 137 6.7 Geometry of a semi-infinite body heated by a uniform heat flux, qo, over an infinite strip with width 2a. ..................... . ................... 138 . 6.8 Dimensionless temperatures T(x+,z+,t+) versus x+, for various values of 2* at t+=.25. ..... 7 .............. 140 . . + + 6.9 DimenSTonless surface temperature T(x ,O,t ) versus dimensionless x+, for various values of dimensionless times. 00000000000000000 0......000.0000 00000000 0000.00. 14] 6.10 Dimensionless isotherms in semi-infinite body heated by a uniform heat flux over an infinite strip at z+=0, for t+=1. . ............. . ........................ 142 LIST OF SYMBOLS a Raduis of the Contact area AC Contact area Aj Area of the surface element j b Raduis of cylinder b+ Dimensionless value of b given by (3.4.5c) B: Matrix defined by (2.5.12) Cp Specific heat C: Matrix defined by (2.5.8a) Ci(°) Cosine integral 0’ Vector defined by (2.5.8b) en(-) Truncated exponential function erf(-) Error function erfc(-) Complementary error function E' Vector defined by (2.5.8C) E(-) Complete elliplic integral of the 2nd kind defined by (3.4.4b) E](-) Exponential integral F Vector defined by (2.5.8d) F(t) Time-varying input function hs Heat transfer coefficient Rh Conductance matrix at time tM = M-At H(X,P) Function defined by (6.3.17) ierfc(-) Iterated error function xi Function defined by (3.4.5b) Bessel function of the kth order Thermal conductivity Complete elliptic integral of the lst kind defined by (3.4.4a) Time index ‘ Outward normal direction Number of surface elements Variable defined by (6.3.10c) Heat flux Constant surface heat flux Center line heat flux Area averaged interface heat flux Heat flux vector at time tM = M-At Total heat flow through the contact area as a function of time position vector Radial coordinate Dimensionless radial coordinate Interior region Thermal constriction resistance Normalized thermal constriction resistance Coordinate along the boundary Surface boundary Sine integral Time Dimensionless time ‘Temperature Surface temperature xii Ambient temperature Initial temperatures Area average interface temperature Initial temperature vector Average temperature of the contact area as a function of time Normalized Contact area temperature Instantaneous temperature Dimensionless temperature defined by (6.3.10d) Variable defined by (6.3.8) Width of the slab (see Fig. 2.3) Cartesian coordinate Dimensionless x defined by (6.4.2a) Dimensionless x defined by (6.3.10a) Cartesian coordinate Axial and Cartesian coordinates Dimensionless 2 defined (6.4.2b) Dimensionless 2 defined by (6.3.10b) Thermal diffusivity Variable defined by (4.4.2) Euler's Constant Incomplete gamma function A point along the surface boundary Temperature rise due to unit heat flux defined by (3.3.4) Dummy variable Density Temperature rise for unit heat flux (influence function) xiii Instantaneous temperature for unit heat flux Temperature rise at element k due to unit heat flux over element j at time ti = i-A for body m. Temperature rise at interior position (x,y) due to unit heat flux at element 3 at time tM = M-At Influence Matrix at time ti = i-At Fundamental solution Flux based fundamental solution Temperature based fundamental solution xiv CHAPTER 1 INTRODUCTION AND DEVELOPMENT OF INTEGRAL EQUATIONS 1.1 Introduction The transient surface element method (SEM) is a new numerical method for solution of linear transient two- and three-dimensional heat transfer problems. Its development originated with the need of solving certain transient composite body problems. There are numerous cases for which similar or dissimilar bodies are attached one to the other over a part of their surface boundaries. Some examples in conduction are problems associated with the intrinsic thermocouple, contact con- ductance and fins when the transient temperature distribution of the fin including the base is of interest. Other related problems are those with mixed boundary conditions. See Figures 1.1-1.3. In general it is difficult to obtain analytical solutions for such problems. The SEM described and developed in this dissertation is particularly suited for such problems compared with other numerical procedures. The most widely used numerical method in heat transfer is the finite difference method (FDM) [1,2]. It involves approximating the partial differential equations by simpler, localized algebraic ones valid at a series of nodes within the region. The method is extremely valuable for problems involving composite bodies, nonhomogeneous boun- dary conditions, and nonlinearity in the differential equations or boundary conditions. The finite element method (FEM) is similar to the FDM. The main difference between the FEM and the FDM is in the way of constructing the algebraic equations. The FDM involves approximating derivatives in differential equations, while in the FEM the approximating equations wire wire nsulated insulated various boundary 1 conditions (a) (b) Figure 1.1 Some geometries for intrinsic thermocouple problem. ’ r 5 F insulated fl "£22222:r ' 3 b2 ::g . o g /———JK (a) (b) Figure 1.2 Some geometries for contact conductance problem. prescribed heat flux history T81 . convective insulated flux boundary f.- r condition x E? constant temperature or insulated (a) (b) X Figure 1.3 Some geometries illustrating mixed boundary conditions. are cast in an integral form [3]. This approximation results from a minimization process based upon the theories of variational calculus. In general, the FEM is best designed to handle complex boundaries, while the FDM is superior for complex equations. For the problems mentioned above, the use of the FDM or the FEM is not entirely satisfactory. This is partly due to the necessity of setting up extremely fine grids near the interface, and many large grids further from the interface. Further, both methods involve whole-body discreti- zation schemes which require, finally, the solution of very large system of algebraic equations, especially for two- and three-dimensional problems. These methods unavoidably generate the solution at all internal nodes, whether or not this information is needed. This causes significant economic disadvantages for many applications where only interface results are of interest. Closely related to the SEM is the boundary integral equation method (BIEM), also referred as boundary element method (BEM) which.has become very popular in recent years. It has been used in a variety of engineering problems such as solid mechanics, fluid flow, seepage, soil mechanics, water waves, heat conduction, electrical problems and a broad range of other applications [4-8]. The method utilizes Green's theorem to formulate the problem described by a partial differential equation in a given region with some specific boundary conditions as an integral equation which applies only to the boundary of the region. Basic building blocks used in the BIEM are source solutions (Green's functions) for infinite homogeneous bodies. The main advantage of the method over alternative numerical methods such as the FDM or FEM is the reduction by one of the dimensionality of the problem under con- sideration, i.e. two-dimensional problems are reduced to one-dimensional problems and similarly three-dimensional to two-dimensional ones. This results in considerable savings in the input data and the computer time required to run the problem. A disadvantage of the BIEM is the presence of singularities at the boundaries. The BIEM is well-suited for solving steady state problems with infinite domain and irregular shaped boundaries. A number of papers have been written for steady state heat conduction problems [9-12]. Schneider [9] examined the constriction resistance problem using the BIEM for isothermal rectangular, hollow rectangular, and circular- annular contacts located on a semi-infinite body. He further (in a joint paper with LeDain [10]), introduced the method to solution of steady state heat conduction problems with special attention given to the corner problem+. Khader [11] and Khader and Hanna [12] employed BIEM in conjunction with Kirchhoff's transformation to solve nonlinear steady state heat condution problems. Application of the BIEM to transient problems have received less attention compared to the steady state problems. This is due to the complexity resulted by having the independent variable of time. There are two basic ways of handling the effects of time. One is to temporarily eliminate time as an independent variable by utilizing the Laplace transform and then solving the problem in the transform space by using BIEM. The time solution is then obtained by numerical trans- form inversion. This was the approach taken by Rizzo and Shippy [13] to solve the problem of heat conduction in an infinite cylinder of I This situation arises in the problems with mixed boundary conditions and at geometric corners [10]. an isotropic medium. The other approach (similar to the one used in this dissertation) is to treat the time directly in the same manner as the spatial coordinates are treated, integrating numerically over the time as well as over the boundary of the body [14-16]. Shaw [14] uti- lized the direct approach to investigate the heat conduction in a cir- cular sector of an isotropic medium. A similar approach was taken by Chang, et a1. [15] to treat anisotropic heat conduction in the transient case with heat generation. As a test they applied their method of solu— tion to three specific problems of i) a long square prism, ii) a long circular cylinder, and iii) a hollow eccentric cylinder. Recently, Nrobel and Brebbia [16] employed BIEM to solve three-dimensional axisym- metric transient heat conduction problems of i) a solid circular cylin- der with convection, ii) a prolate spheroidal solid with zero initial temperature subjected to a unit surface temperature at time zero, and iii) a solid sphere with time dependent boundary conditions. The time integrations were performed directly by dividing the entire time domain into small intervals. Temperature and heat flux were assumed to be constant over each interval. This assumption made possible the analy- tical evaluations of the time integrals by using series expansions. For the problems shown in Figures 1.1 through 1.3 (bodies connected over relatively small area) the SEM is superior to the FDM, FEM or BIEM. In the SEM only the interface between the two geometries required dis- cretization as opposed to the discretization of the whole domain required in the FDM and FEM or discretization of the whole boundary in the BIEM. This in turn reduces the size of numerical calculations and computer time. Further, the SEM does not require any modifications or special handling of points near the domain boundaries unlike the above mentioned alternative methods. The SEM uses Duhamel's integral and involves the inversion of a set of Volterra integral equations, one for each surface element. Though the method is limited to linear regionsit can be used for nonlinear boundary conditions. _ Two types of kernels (building blocks or fundamental solutions) can be employed in the SEM; temperature-based and heat flux-based. The SEM requires that these "building blocks" or kernels be known for the basic geometries under consideration. For example, for the con- tact conductance problem (see Fig. 1.1a) solved in Chapter 3, solution for a constant heat flux over a disk-shaped area on a semi- infinite body is needed. For many geometries the kernels are known or can be obtained Simply by analytical or numerical procedures. Yovanovich [17] suggested the name "surface element" and did early work on a steady state form of SEM. Keltner and Beck [18] and later Beck and Keltner [19], were the first to employ the SEM for transient problems. They have considered only one element along the interface and utilized the Laplace transform technique to obtain "early" and “late" time analytical solutions for certain cases [18,19]. Both types of kernels have been used in their solutions. In this dissertation a multinode transient surface element method for two-dimensional heat conduction problems with linear boundary con- ditions is developed and formulated. Only heat flux-based kernels are considered. three different problems are solved and the results are compared with those obtained by the other investigators. In Chapter 1 mathematical descriptions for two-dimensional heat conduction problems with various types of boundary conditions are presented first. Next, the Duhamel's theorem is introduced and its application to arbitrary time and space variable heat flux boundary conditions in a two-dimensional region is developed. Finally, integral equations with temperature-based kernel and heat flux-based kernel are discussed. In Chapter 2, the multinode surface element method formulations for two arbitrary geometries in perfect (or imperfect) contact over part of their boundaries are developed and described. In Chapters 3 and 4, the SEM presented in Chapter 2, is employed to solve the problem of two semi-infinite bodies initially at two different temperatures suddenly brought together over a small Circular region and insulated elsewhere (contact conductance problem), and the intrinsic thermo- couple problem, respectively. In Chapter 5, the SEM is further utilized to solve the problem of a semi-infinite body with mixed boundary condition of a step change of the surface temperature over an infinite strip and insulated elsewhere. The kernels required in this chapter are obtained from the exact solution for the problem of a semi-infinite body heated by a con- stant heat flux over half the surface which is presented in Chapter 6. Finally, the closure and conclusions are given in Chapter 7. 1.2 Mathematical Description In this section the equation and the different types of boundary conditions associated with the boundary value problems of heat con- duction for two-dimensional bodies are stated and briefly discussed. Consider the boundary value problem of heat conduction for a two dimensional region R, with the boundary S as shown in Fig. 1.4. The body is assumed to have temperature-independent thermal properties and to be homogeneous and isotropic with no heat generation. The Figure 1.4 Geometry Showing a two-dimensional heat conduction problem. .../- Tshht) S (R) y Y X (a) . (b) S T s [*5 f— 5 7’s [— ] Ts(s.t) 11 S (R) 4 Ir T =0 5 (R) 52 y L—w n. S 3 3T -1 n 1L 5n—=k Q(S.t) x s -k a =h (T '7...) s (C) "S S S (d) -k 3%.: ghS(TS-T°°) Figure 1.5 Geometries illustrating various types of boundary conditions in a two-dimensional region. differential equation is, 2 = l_3T x, ,t) V T(x,y.t) a at (1.2.1) with the initial condition of T(x,y.t) + Ti(x..v) (1.2.2) as t + 0 where T is the temperature, x and y are the space variables, t is the time, and a is the thermal diffusivity. Type of Boundary Conditions In this thesis the primary concern is with linear boundary condi- tions. Linear boundary conditions, in general, can be separated into four different types. I. Prescribed Surface Temperature (First Kind) The surface temperature of the boundaries can be specified to be a constant, a function of time or position, or a function of both position and time, T = Ts(s,t) on S (1.2.3a) where 5 denotes location along the surface boundary S. See Fig. 1.5a. If the temperature at the boundary surface vanishes, the homogeneous boundary condition of the first kind is obtained, T = O on S (1.2.3b) II. Prescribed Heat Flux Across the Surface (Second Kind) In the second type of boundary conditions the normal derivative of temperature is Specified to be a constant, or function of time t, 10 or position 5, or both 5 and t at the boundary, 3T S’t = k-1q(S,t) on S (1.2.4a) ans where §%-denotes differentiation with respect to outward pointing nor- 5 mal to the surface boundary S. See Fig. 1.5b. This boundary condition is equivalent to prescribing the magnitude of the heat flux along the boundary. If the normal derivative of the temperature at the boundary surface vanishes, the homogeneous boundary condition of the second kind is obtained (insulation), 33-3—5151: 0 on S (1.2.4b) 305 III. Heat Transfer to the Ambient by Convection (Third Kind) If the flux across the surface boundary is proportional to the temperature difference between the surface boundary and the ambient, the boundary condition is, -k filgfi—tl = hS[TS(s,t) - Too(s,t)] (1.2.5a) S where T00 is the ambient temperature and hs’ the proportionality con- stant, is the heat transfer coefficient. See Fig. 1.5c. This equation can also be written as, 3T(S,t) k an + hsTS(s,t) (1.2.5b) = hSTm(s,t) s f(s,t) on S As hS + 0 this case tends to the boundary condition of case 11, and as 11 hS + w it tends to the boundary condition of case I. If the ambient temperature is zero, the homogeneous boundary condition of the third kind is obtained, k 11%;), hsTs(s,t) = 0 on S (1.2.5c) 5 IV. Mixed Boundary Condition In this case all or some of the previous types of boundary condi- tions may be specified over different portions of the boundary S. T = Ts(x,t) on S], k 2%IELEI = q(s,t) on S , IIS 2 1.2.6) aT(s,t) = ( k ans + hSTS(S,t) hSTm(S,t) on 83, T = O on $4, and so on. See Fig. 1.5d. The boundary conditions described above cover most cases of physi- cal problems. There are also the radiative boundary conditions with heat flux obeying the fourth-power temperature law, the natural convec— tion boundary condition, and the interface conditions associated with phase change. Such boundary conditions are nonlinear and are not considered. However, it Should be noted that the surface element method is capable of treating nonlinearities at the boundaries and approximately inside domains for certain cases. 12 1.3 Duhamel's Theorem Duhamel's theorem provides a means for solving boundary value problems of heat conduction with time dependent boundary conditions. It utilizes the solution to a corresponding fundamental problem with time-independent boundary conditions. This method is a useful tool for obtaining the solution of problems with transient boundary condi- tions whenever the solution to the corresponding fundamental problem is available. Duhamel's theorem has been expressed in different forms. Briefly it states that if w(F,t) is the solution to a linear system initially at zero temperature, due to a unit stepwise input, then the solution to the same system initially, at zero temperature due to a time- varying input F(t) (instead of a unit step) is, t Of F(x)ip(F,t-X)dx (1.3.1) — -a_ T(r’t) — at where F’is the position vector, t is the time, and X is a dummy variable. The input function, F(t), can be any type of time-dependent boundary condition (for instance, prescribed surface temperature, ambient tempera- ture or heat flux), or heat generation. Using Leibnitz's rule for dif- ferentiation of an integral, an alternative form of (1.3.1) can be obtained, 1: _ T(T~',i) = of m) Mfg—{fill d). - (1.3.2) It also has been conventional to treat problems with arbitrary space-variable conditions by using Duhamel's method with integrating over space [20]. In the following section it is shown that Duhamel's method can be extended to simultaneous variation of both time and space 13 conditions. 1.4 Derivation of Duhamel's Integral Equation for Arbitrary Time and Space Variable Heat Flux Boundary Conditions in a Two- Dimensional Region Consider the two-dimensional heat conduction problem given in Section 1.2 with the heat flux boundary conditions given by equation (1.2.4a) for an arbitrary geometry as shown in Fig. 1.5b. This inves- tigation is restricted to systems initially at zero temperatures (Ti = 0), because most practical problems can be reduced to this case. For simplicity it is assumed that q(s,t) is non-zero only over the portion of the boundary S from s = O to s = L, and the other portion of the boundary is insulated (q(s,t)=0, for s > L). See Fig. 1.6. The objective is to find the expression for the solution of the above problem using Duhamel's integral method. In the first step the solution to the fundamental problem is found. The fundamental problem is identical to the above problem with the excep- tion that the variable flux boundary condition q(s,t) is replaced by a special unit step function. It is described by the following equations: av 2 =_1._El V Wq a at (1.4.1) wq(x.y.0) = 0 (1.4.2) 3P __fl.= k ans 0 for t < O or s < n (1.4.3) = 1 for t > O and n < s < L where n is a point along S between s = O to s = L, and wq(x,y,n,t) is the temperature rise at position (x,y) and time t caused by a unit 14 q(s,t) S S ’1 y / X "S q=0 (a) F(SJ) s,t) dn q( _-> S M /\< (R) (b) Figure 1.6 Geometry Showing a two-dimensional region heated by an arbitrary heat flux only over a portion of its boundary. 15 step change of heat flux at time t = 0, from n to L (n < s < L), which is called the Flux Based Fundamental Solution (FBFS). It is indicated in Fig. 1.6 by the cross-hatched portion. Notice that, for fixed (x,y) and t, wq(x,y,n,t) decreases as n increases, WQIXSYSUIt) > wq(xay,n+dnat) (1-4°4) Temporarily let, n be fixed and consider the variation of the heat flux with time only, From the fundamental solution the temperature rise at position (x,y) and time t due to a unit step change of heat flux at time X is, vq(x.y.n.t-x) (1.4.5) where (t-X) is the time that has elapsed since the step at X. Also the temperature rise at time t due to a unit step change of heat flux at time X+dX is, wq[x.y.n.t-(A+dl)] ’ (1.4.6) Then from (1.4.5) and (1.4.6), the temperature rise at position (x,y) and time t due only to a unit step change in q for X < t < X+dX is 'dxw (XSYSUSt'A)=wq(XSYSnat'A)'wq[X,y,n,t-(A+dA)] (104-7) 9 where dX is a differentiation operator for X. Notice that wq(x,y,n,t-X) is greater than uq[x,y,n,t-(X+dX)]. Using (1.4.7), the temperature rise at position (x,y) and time t due to the value q(n,t) for A < t < X+dX and n being fixed is given by 3¢q(Xsyin.t-1) TQ(n9A)dqu(XSYSnat-A)='q(naA) 3A dA (1.4.8) 16 for small dX. Since the problem is linear, superposition can be em- ployed and the total effect of all step changes of heat flux over small dX's from time zero to time t is Simply found by integrating (1.4.8) from O to t. Hence one can write, t 3Wg(X.Y9n,t-X) wg(x,y.n.t) = - i q(n.X) 31 d1 (1.4.9a) 0 Since, 3Wq(X,y,nst-A) 3W 31 ‘ ' a t-X and a a at- awq(x y n A) = 84g at aTt-x) it is evident that, 3¢q(fsysnst-41 = - 8wq(xay:nat'x) 8X 8t Therefore (1.4.9a) may be written as, t 3¢q(X:y,n,t-A) w'(X.y.n.t) = I q(n.A) at dX (1.4.9b) q o where Ta is the temperature rise for the case that q is zero for s < n, and is uniformly distributed over space for n < s < L. Further, one can Show that the temperature rise for the case that the heat flux, q, is zero for s < n+dn and is uniformly distributed for n+dn < S < L, is, 17 t 3Wq(xay,n+dnat'x) w;(x.y.n+dn.t) = of Q(n.X) at dX (1.4.10) Using (1.4.9b) and (1.4.10), the temperature rise due to a uniform heat flux q, between s = n and S = n+dn and for t > O is, 'dan(XSYanat) = wq(X.y,n,t) ' wq(x,yan+dnat) (1.4.11) _ awé(xay’nat) - 3n dn where dn is a differentiation operator for n. (Notice that wé(x,y,n,t) is greater than wé(x,y,n+dn,t).) Introducing (1.4.9b) and (1.4.10) into (1.4.11), one can Show, t 32Wq(XSYSnat'A) -dan(X.y,n.t) = -0f H(n.l) atan dldn (1.4.12) Again superposition can be employed and the total effect of the varia- tion of heat flux from S = O to s = L can be found by integrating (1.4.12) over space from S = O to s = L. L t 32¢q(XSY.n,t-X) T(x.y.t) = -0; OI q(n,A) Eta” dan (1.4.13) In this problem it was assumed that only a portion of the surface boundary is exposed to heat flux with the remainder being insulated. However, if none of the boundary S is insulated, the first integral in (1.4.13) extends over entire boundary S. Furthermore, if the ini- tial temperature of the system is Ti instead of being zero, the solu- tion becomes 18 t 32w (XSYSnat'A) - = _ q T(x.y,t) Ti g of q(n.A) atan dan (1.4.14) In the above problem the input function is the heat flux along the boundary which varies with both space and time, and the solution is in terms of the FBFS, wq. If, however, the surface temperature were known along the boundary as the input function, (instead of heat flux), then in a similar manner to that described above, the solution in terms of the Temperature Based Fundamental Solution (TBFS), OT, can be ob- tained as, t 132wT(x,y,n,t-A) T(x.y.t)-Ti- - é of [Ts(n.x)-TiJ anat dan (1.4.15) Equations (1.4.14) and (1.4.15) are rather general expressions for the case that the input function varies with both space and time in a two-dimensional region. To the author's knowledge, they do not appear anywhere in the open literature. Colladay [21] has shown the Duhamel's integral equation with flux-based kernel for flow between heated paral- lel plates. His equation, however, is different from what is given here. In the next section the application of the above equations to the solution of two-dimensional transient heat conduction problems is discussed. 1.5 Integral Equation Formulation The Duhamel's integral equations given in Section 1.4 are now employed to obtain the solution to equation (1.2.1) subject to the various types of boundary conditions discussed in Section 1.2. The solution can be determined in two different ways, either a) using wT, 19 the TBFS, or b) using wq’ the FBFS. To compare the two approaches and discuss their utility for each particular type of boundary condition, both form of solutions given by (1.4.14) and (1.4.15) are considered. For the problems with the boundary conditions of the first type, where the temperature is solely specified everywhere along the boundary, (1.2.3a), the right hand side of (1.4.15) is known, and one can solve for the temperature history of any interior point of R, by direct integration. If, however, (1.4.14) is employed instead of (1.4.15), the direct evaluation of T(x,y,t) is not possible, because of the un- known heat flux, q, in the right hand side of this equation. In this case an inverse problem must be first solved for the unknown surface heat flux which is the required information needed by (1.4.14) to find T(x,y,t) at any interior point. (The inverse problem is discussed further in the solution of the mixed boundary-value problems.) There- fore, in the problems with the first type of boundary condition, (1.4.15), is more appropriate than (1.4.14). On the other hand if the boundary condition is of the second type where q is solely specified along the boundary S, (1.2.4a), then the right hand side of (1.4.14) is known which leads to evaluation of a direct integral. In this case (1.4.14) is more appropriate than (1.4.15). For the boundary conditions of the third and fourth type where neither surface temperature nor its normal derivative are completely known over the entire boundary S, none of the above equations can be used directly to obtain temperature history for any interior point. An example is given to illustrate this case better. Consider the homogeneous convective boundary condition given by (1.2.5c) where the temperature and its normal derivative are related through a linear expression along the boundary as, 20 k §I%%SEI.+ hSTS(s,t) = 0 on S (1.2.5c) S substituting for q in (1.4.14) from (1.2.5c) one can write, t 32W (XSYSUSt'A) T(x,y,t)-Ti: - f I h T (n,1) q dan (1.4.16) 5 0 S S atan Equation (1.4.16) cannot directly be integrated for T(x,y,t), Since Ts(s,t) inside the integral is unknown. In other words the number of unknown functions in (1.4.16) is more than one, Ts(s,t) and T(x,y,t). However, for a point along the boundary S, (1.4.16) reduces to, 32 W (nat'x ) T (s, t)- -j of h ST w(n, q dan (1.4.17) S s atan which is a Volterra integral equation of the second kind with the only unknown function, Ts(s,t), both inside and outside the integral. As an inverse problem (1.4.17) can be solved numerically for Ts(s,t). Once the surface temperature, TS(s,t), has been determined, the solu- tion to the interior temperature history, T(x,y,t), can be obtained by substituting Ts(s,t) into (1.4.16). Hence, for the problems with mixed or convective boundary condi- tions the temperature history at any interior point (x,y) can be determined in two steps: a) Find the boundary information by solving an inverse integral equation, and b) using the boundary data obtained in a), find the interior temperature history by using a direct integration. In this thesis the solution of two-dimensional heat conduction 21 problems only based upon wq(FBFS) is considered. Equation (1.4.14) is the basic starting point in the development of the SEM formula in the following sections. Furthermore, the problems with the first and second type of boundary conditions are not discussed herein, since their solutions can be found explicitly by direct integration. CHAPTER 2 SURFACE ELEMENT METHOD FOR TWO BODIES IN CONTACT 2.1 Introductory Remarks The heat transfer between two bodies with perfect or imperfect contact at the interface is of fundamental importance and it has accord- ingly received considerable attention over the last two or three decades. It is important in the problems involving electric contact, electronic cooling, welding, fins, contact conductance, and many other applications for which two similar to dissimilar bodies are attached one to the other over small parts of their surface boundaries. Some examples are given in Chapter 1. In general it is difficult to ob- tain analytical solutions to such problems. In this Chapter a transient multinode surface element method for two arbitrary geometries contacting over part of their surface boun- daries is. developed and formulated. Both cases of perfect contact and imperfect contact are considered. The method starts with the Duhamel's integral equation given by (1.4.14) which is then approximated numerically in a piecewise manner over time, and the boundaries and the interfaces of interest. The method is superior to the other numerical procedures such as FDM, FEM, or BIEM, for the particular problems mentioned above. In this method only certain parts of boundaries (interface) need to be discretized as opposed to the whole body discretization required in the FDM and the FEM or discretization of the whole boundaries needed in the BIEM. 22 23 2.2 Problems to be Considered To illustrate the capability and limitations of the method, the three different geometries given below are investigated. The first geometry involves two semi-infinite bodies initially at two different temperatures suddenly brought together over a small circular region and insulated elsewhere (contact conductance problem). See Fig. 1.1a. The second geometry has a semi-infinite insulated cylinder attached to a semi-infinite body (the intrinsic thermocouple problem). See Fig. 1.2a. The third geometry is a semi-infinite body with the mixed boundary condition of a step Change of the surface temperature over an infinite strip and insulated elsewhere. See Fig. 1.3a. In each case the solution to the interfacial heat flux is empha- sized. For the third geometry, however, the solution to the interior temperature history is also formulated. 2.3 Discretization over Space In order to solve numerically the integral equation given by (1.4.14), the surface boundary is divided into N finite surface elements Asj, as shown in Fig. 2.1. (Notice that only the parts of the boundary with non-zero values of heat flux need to be discretized). Equation (1.4.14) can be written as, t N 3 W (XSYSnst‘A) T(x,y,t)-T =- I [.2 I q(n.x) q~ at,“ danx (2.3.1) 24 “(s,t) Asj ' Qj(t) ——’ s / i M T I 2/1/ 1 SO=0 s1 52 Sj-lsj SN (R) Figure 2.1 Geometry showing discretization over the heated portion of the surface boundary. (ch-(t) ' ___. q ]q1 q2 3 I N l 2 3 N -—-i-S Figure 2.2 Uniform heat flux assumption over each element. 25 In the simplest form of approximation, the heat flux is assumed uniform over each element (but different, in general, at each element. See Fig. 2.2), so that, t N a Si T(x,y,t)-T1.=- f { 2“. qJ-(I) gi- [wq(x.y.n.t-A) I 11d) 0 i=1 5. 3-1 (2.3.2) It N ()3 ( ) =- { Z q. A -—— [Aw x,y,t-X ]}dX 0 jzl J at qj where Avqj (x.y.t) = wq(x.y.sj,t)-wq(X.y.sj_1,t) (2.3 3) Further, if the temperature rise at position (x,y) due to a unit step increase in heat flux q at the element j is denoted ¢j(x,y,t), it can be shown that -Av (x,y,t) = ¢-(X.y.t) (2.3.4) Qj J (See Sec. 1.4). Using (2.3.4) in (2.3.2) results in t 3¢'(Xay’t'k) I qj(A) J at dX (2.3.5) IIMZ T(x,y,t)-Ti= j 1 O which gives the temperature rise at location (x,y) and time t due to the effect of N heat flux histories q1(t), q2(t), ..., qN(t). The func- tion 4j(x,y,t) is the basic building block solution needed in the above expression and is termed an influence function. 26 2.4 Surface Element Formulation for Two-Bodies in Contact In this section the surface element method formulations for two arbitrary geometries in contact over part of their boundaries are de- veloped and described. a) Perfect Contact Two different geometries of a semi-infinite body and a semi- infinite slab initially at uniform but different temperatures are brought together into perfect contact over the interface of width w. One side of the slab, x = 0, is held at zero temperature and the other parts of the boundaries are assumed to be insulated. See Fig. 2.3. The bodies may have different thermal conductivities, k, and density- specific heats, pc The upper body is referred to aS region 1(y>0) p' and the other body as region 2(yO, as [x|+ w (2.4.3a) T2=O for t>0, x=o, y<0 8T2 , (2.4.3b) T: o for t>O, x=w, y<0 T1=T2 3T1 8T2 k1 5y—-- k2 3y for t>0, nggw, y=O (2.4.4a) 8T1 F: 0 for t>0, x<0 and x>w, y=0 T1=Ti1 for t>0, as y+ w (2.4.4b) T2=T12 for t>O, as y+~m where T1 and T2 denote the temperature distibutions, k1 and k2 repre- sent the thermal conductivities, a1 and 02 refer to the thermal diffu- sivities, x and y are the Space coordinates, and t is the time. Analysis First, the interface is divided into N finite surface elements (each being an infinite strip), ij, as shown in Fig. 2.4. (If the system is symmetric only half of the interface needs to be dis- cretized.) It is assumed that there is no spatial variation of tempera- ture or heat flux over each element (uniform approximation). The heat fluxes q1(t), q2(t), ..., qN(t) are arbitrary functions to be found. The heat flux qj(t) which leaves body 2 in Fig. 2.4 is the same heat flux that enters body one over the region x=xj_1 to x=xj, that is, 28 Semi-infinite body 1,757/Ufl/' ' 7'77’?’ x Semi-infinite slab )— insulated ///)'I w :0 . Figure 2.3 Geometry of semi-infinite slab attached to a semi-infinite body. 1y insulated Figure 2.4 Possible distribution of surface elements for connected semi-infinite body and slab. -k ———-- -k ———-= q (t) for t>O, xj_1§x5xj, y=O (2.4.5) Using (2.3.5), the temperature at element k in body one and at time t can be given by, Tk1(t)=Til+jgl oftqj(t) 33113;;l32,dx (2.4.6a) where Tkl E T(xL,O,t) . x) = xk-Axk/Z (2.4.6b) and 4);)(t-1) E ¢§1)(x;,0,t-X) (2.4.6c) is the temperature rise at element k and time t due to a unit step heat flux at element j of surface 1. ¢£;)(t) is the basic building block needed for body 1. It Can be found from the known solution given in Chapter 6, for a constant heat flux over an infinite strip of a semi- infinite body. See Fig. 2.5a. Similar to (2.4.6a), an integral equation can be given for the kth surface element of body 2, (2) N t 36k. (t-X) Tk2(t)=T12-_§ f qj(X) -——l—5$—-—-dx (2.4.7a) 3-1 0 where Tk2 E T2(x&,0,t) (2.4.7b) 30 and ¢£§)(t—X) E ¢§2)(x£,0, t-X) (2.4.7c) is the temperature rise at element k and time t due to a unit step heat flux at element j of surface 2. The minus sign before the summation in (2.4.7a) is used because the heat flux is pointing outward from body 2. The basic building block for body 2, 6:?)(t), can be found from the solution of a semi-infinite slab heated over an infinite strip with zero temperature on one side and insulated on the other side. See Fig. 2.5b. This solution can be obtained from the known solution for a con- stant heat flux over an infinite strip of a semi-infinite body by using the method of images (see Chapter 6, and [22].) For the case where the bodies are in perfect contact, one can write, Tk1(t) = Tk2(t) for k=1. 2, ..., N (2.4.8) By introducing (2.4.6a) and (2.4.7a) into (2.4.8), a set of integral equations for k=1, 2, ..., N can be found, N -T. = 2 oft q. (t) )-——l—————-dx (2.4.9) where = (1) (2 ) iijt) ikj (t)+ ¢k (t ) (2.4.10) Equation (2.4.9) represents a set of Volterra equations of the first kind which can be solved simultaneously for N unknown heat flux his- tories q1(t), q2(t), ..., qN(t). 31 heat flux q=l insulated insulated 1 heat flux q=l I (a) (b) Figure 2.5 Basic building blocks for the geometries of semi-infinite body and semi-infinite slab. quh) qk.M-l 72 I I \ I I /—'qu I I (he I l/ I I I t4 ' J I I;1 t 't t t M-l M 0 1M 2 Figure 2.6 Geometry illustrating uniform heat flux assumption over each time interval. 32 b) Imperfect Contact Even though the perfect contact is a common interface assumption, it will only be valid for very intimate contact, such as a soldered joint. For imperfect contact (2.4.8) is replaced by qk(t)=hk(t)[Tk2(t)-Tk1(t)] for k=1, 2, .... N (2.4.11) where hk(t) is the time-variable contact conductance for surface ele- ment k. The above relation tends to the case of perfect contact (given by (2.4.8)), as hk+ m. It also includes the cases of convection, prescribed heat flux and prescribed temperature. By introducing (2.4.6a) and (2.4.7a) into (2.4.11), a set of integral equations for heat fluxes qk(t), k=1, 2, ..., N, can be obtained, qk(t) N t 36 .(t-X) - = .__EL_____ T1.2 T1.1 EETETA-jil of qj(t) 3t dX (2.4.12) for k=1, 2, ..., N Equation (2.4.12) represents a set of Volterra equations of the second kind with the unknown heat fluxes, qk(t)'s, appearing both inside and outside the integrals. These integral equations can be solved simul- taneously in an inverse fashion for the unknown heat fluxes. The method of the solution is described for the case of imperfect contact which includes the other cases as well. 2.5 Solution of the Simultaneous Integral Equations The Volterra equation (2.4.12) can be approximated by a system of linear algebraic equations by replacing the integrals with suitable quadrature formulas. As the first step, the time region 0 to t is 33 divided into M equal small time intervals, At, so that tM represents the value of t at the end point of the Mth interval, d 1H tM = MAt (2.5.1) Equation (2.4.12) can be written as qk(tM) N M ti 36k.(tM-x) T. -T. = -+ z 2 f q.()) .J dX 12 11 hkitM) j=1 t=1 ti-l J at (2.5.2a) for k=1, 2, ..., N where t0 5 O (2.5.2b) In the simplest form of approximation the heat flux histories qj(t) are assumed to have constant values in each time interval (see Fig. 2.6) so that qu N M Ti==fiRMA-jil iEI qji A¢kj,M-i for k=1, 2, ..., N (2.5.3a) where Ti = TiZ-Ti1 (2.5.3b) qji E qj[(l-1/2)At] (2.5.3c) A¢kj,M-i = ¢kj,M+1-i ' ¢kj,M-i (2-5°3d) In the form given by (2.5.3a), the heat fluxes qu's (for j=1, 2, .., N) can be determined at different time intervals one after another, 34 by marching forward in time for, M = 1, 2, 3, ... While calculating each new time component, the fluxes at previous times, qjl’ qj2’ qj3’ qj,M-2’ qj,M-1 are known for j=1, 2, ..., N. Thus for each time step, equation (2.5.3a) represents a system of N equations with N unknowns qlM’ q2M’ q3M, ..., qNM' The objective is to solve this system for the unknowns qu, for j=1, 2, ..., N. For convenience define, ¢kji E ¢j (x xk,0, t ) (2.5.4) Introducing this definition into (2.5.3a) and noting that, ¢kj0=0’ gives M+gq.¢.=1.+gM21q¢ hKM j=1 3M k31 l j_1 1__1qji kj M- 1 N M-l - jzl 1:1qj1¢kj M+1' _1 for k=1, 2, ..,, N (2.5.5) Equation (2.5.5) is written in standard form with unknowns on the left and knowns on the right. Expressing (2.5.5) in matrix form gives, M-1_ _ M-l = =T'i + Z 3 q - (11+M 2 (2.5.6) i=1 ‘ i=1 M =1)EM with ‘6'" I” P¢111 ¢211 c¢N1i -qNU ¢121 ¢22i ¢N2i _ . 1 HM : diag [fi——- 1M .4 35 ¢ini ¢2Ni (2.5.7a) (2.5.7b) (2.5.7c,d) where 31 and fiM are called the influence matrix and the conductance matrix, respectively. matrices, where on II C” u "N *N Then (2.5.6) III -fl M-1 i=1 M '9'" M-1 - i=1 2 M-i If further E, and 5 are defined to be the nu °w14 5” gnu-‘0 q. 1 can be written as, (2.5.8a) (2.5.8b) (2.5.8c) (2.5.8d) 36 E =‘fi (2.5.9) The solution of (2.5.9) is -1 qM fill fi- ' (2.5.10) The E matrix, multiplier of EM, has to be calculated at each time step if the diagonal matrix fiM’ is a function of time. However, for the case that the contact conductances do not change with time, the E matrix needs to be calculated only once during the entire solution and an al- ternative form of solution can be given as, (see Appendix A) -1._ Q1 = E T, (2.5.11) _ _ = M'l _ =-1_ QM = Mq1 + B [ z qi] - c F for M=2, 3, ... (2.5.11b) i=1 where fi- fi-1 = - cpl (2.5.12) (Since i is not a function of time in (2.5.12), the subscript M is dropped.) For the case of perfect contact where hkMV w, the diagonal conductance matrix, HM, becomes zero which implies that, E = 5 (2.5.13) 1 Introducing (2.5.8a,c,d) and (2.5.13) into (2.5.10), results in a simpler form of solution as, q1= C T'i (2.5.143) 37 qM = M6, - E‘l'F for M=2, 3, ... (2.5.145) The elements of the [NxN] influence matrix 51 are: 1 2 ' ¢kji = ¢£j2 T ¢kjg (2°5'15) If the two bodies in contact have the same geometry and thermal pro- perties then, z 2(13(1) = 2¢(2) (2.5.16) ¢kji kJi kJi It is helpful to display the expression for GM more explicitly. To il- lustrate, the case of perfect contact at the interface with only two elements is considered (N=2). In other words there are two heat flux histories, q1(t), q2(t) to be determined. For simplicity only three time steps are considered (M=3). At the first time step, (2.5.14a) becomes: c c '1 T q11 11 12 i = (2.5.17) q21 C21 C22 Ti where = = (1) (2) ij ¢kj1 - ¢kj1 T ¢kj1 (2°5°18) Solving the above system for q11 and q21 yields T.(C -C ) - ‘ 22 12 (2.5.19a) q11 ’ A 38 = Ti(C11'C21) q21 A where A = C11C22‘C12C21 For the second time step, M=2, (2.5.14b) becomes: "' r- ” ”1" '1 q121 qlll C11 C121 F1 1922- (Q21- -C21 C22- th. Solving (2.5.2 ) for q12 and q22 yields z (2Ti-F1)C22-(2Ti-F2)C12 q12 A z 2 _ F1C22’F2C12 q11 A = (2Ti-F2)C11-(2Ti-F1)C21 q22 A _ F2C11'F1C21 ' 2q21 ' A where 1 = $111q11 + ¢121q21 2 ‘ 2211q11 + ¢221q21 (2.5.19b) (2.5.19c) (2.5.20) (2.5.21a) (2.5.21b) (2.5.22a) (2.5.22b) In a similar manner, for the third time step, M=3, one can write 39 F c -F.C - _ 1 22 2 12 q13 ‘ 3q11 A F c -F c - 2 11 1 21 q23 ‘ 3q21 ‘ A where 2 F1 = ,El[¢1l’3-i qli+¢12.3—i q21'] 2 F2 = iEIE¢21.3-i q11+¢22,3-1 Q21] (2.5.23a) (2.5.23b) (2.5.24a) (2.5.24b) Notice that F1 and F2 are the only terms that should be evaluated at each time step. Because of convolution behavior of the summations given in (2.5.8c) and (2.5.8d), the influence matrices, Ei's, need to be calculated at each time step. Consequently, most of the computation effort is in the evaluation of column matrix 5, particularly as the value of M becomes larger. CHAPTER 3 TWO SEMI-INFINITE BODIES IN CONTACT OVER A CIRCULAR AREA 3.1 Introduction In this chapter the transient thermal response of two semi- infinite bodies, brought suddenly into the perfect thermal contact over a small circular area is considered. The rest of the areas of the contacting planes are insulated. See Fig. 1.2a. The temperatures of the bodies are initially at uniform but different values. The SEM is employed to obtain the transient solutions for the interface heat fluxes and temperatures, and the thermal constriction resistance of the contact area. The solutions cover the entire range of dimension- less times. The results are compared with those obtained by other investigators. Two different cases are investigated; a) identical materials on both sides of the contact plane, b) different materials on the two sides of the contact plane. The former case is similar to the problem of a uniform step temperature change over a disk on the surface of a semi-infinite body and insulated elsewhere (see Fig. 3.1). Previous Work The steady-state problem for this geometry has been previously examined and is well known [22,23]. The transient problem for Case a has been analyzed by several authors by considering a single semi- infinite body with isothermal disk on its surface [24-29]. Normington and Blackwell [24] and later Blackwell [25] were the first to seek the solutions in oblate spheroidal coordinates, in order to eliminate the mixed boundary conditions which occur on the surface if the problem is set up in cylindrical coordinates. They developed an approximate 40 4] solution by using Laplace and Legendre transform techniques. However, the solutions were valid only for long times (t+>4) in the first paper and for short times (t+<.1) in the second one. Keltner [26] used the same coordinates to obtain a one-dimensional approximate solution using the heat balance integral method. He has also examined a finite dif- ference solution in cylindrical coordinates for a one-dimensional averaged model [27]. His solutions are more appropriate for early to middle times. Schneider, et al [28] have developed finite difference solutions in oblate spheroidal coordinates for the two-dimensional axisymmetric case. In a recent paper [29], Marder and Keltner examined the problem by using the method of separation of variables. The composite Case b, where the materials on each side of the contact area are not identical. was first studied by Heasley [30] in an approximate manner. He replaced the region of contact by a per- fectly conducting sphere and solved a one-dimensional problem in spher- ical coordinates. His solution is approximately valid for long times. Other work has been done by Schneider, Strong and Yovanovich [3]], Sadhal [32], and Beck and Keltner [l9]. Schneider et al [3T], have employed finite difference techniques posed in oblate spheroidal coor- dinates to obtain the numerical solutions for the two-dimensional axi- symmetric case. They have used about 200 nodes within each of the bodies. Their results as reported by Beck and Keltner [19], are too high at the early times. Sadhal [32] has solved the problem analytically by using Laplace and Legendre transform techniques. His solution is valid for large values of dimensionless time (t+>10). Beck and Keltner [l9], were the first to employ the surface element method to solve the problem. They have considered only one element across the interface and 42 [— L] U m step change of temperature ;--r Figure 3.1 Geometry of a semi-infinite body with a step change of the surface temperature over a circular area. ' 1 plane of aN symmetry insulated Figure 3.2 Distribution of surface elements for two semi-infinite bodies in contact over a circular area. 43 solved the problem analytically by utilizing Laplace transform techniques. Two expressions have been recommended. One is based on a heat flux form of Duhamel's integral (equation (22) of [19]), which is for early times and the other is based on a temperature-kernel form of the Duhamel's in- tegral which is for late times (equation (56) of [19]). Most of the above-mentioned solutions are restricted to either the time or the spatial validity [24727,30,32,19]. The finite difference approaches by Schneider et. al. [28,3l] also have some difficulties such as: the effort in setting up large grids, and the restricted dimen- sionless time step that can be used. 'They require a large number of nodes and consequently required considerable computer time if the accurate results are desired at early times and large times. 3.2 Statement of the Problem Two semi—infinite bodies initially at uniform but different tempera- tures are brought together into perfect contact over a Circular region of radius a. The rest of the surfaces of the contacting planes are assumed to be insulated. See Fig. 1.2a. The bodies can have different thermal conductivities k, and density-specific heats, pcp. The upper body is referred to as region 1 (z>0) and the other body as region 2 (z<0). The initial temperatures are denoted by T1.1 and Ti2 for region 1 and 2, res- pectively. The bodies are assumed to have temperature~independent thermal properties. The describing differential equations are 2 2 a T1 +.l 8T1 + 3 T1 = l..3:l. 3r2 r 3r 322 a1 a (3.2.1) 321 3T 321 BT 2+1 2+ 2=L—.2_ 2 r 3r 2 a 3 ar 32 .2 44 subject to, T1 = T1.1 for t=O, r39, zip (3.2.2a) T2 = T1.2 for t=0, r39, 259 (3.2.2b) and T1 - T1.1 for t>0, as r+w and z+w (3.2.3a) T2 = T].2 for t>0, as r+w and 2+-00 (3.2.3b) 8T1 3T2 T = 5'".— = 0 at r=0 f0)" t>0 (3.2.4) T1 = T2 for t>0, Ogrga, z=0 (3.2.5) 3T1 3T2 k -—-= k 1 32 2232 8T1 37 = 0 (3.2.63) for t>0, r>a, z=0 .52. _ 0 (3.2.6b) where T1 and T2 denote the temperature distributions, k1 and k2 represent the thermal conductivities, a1 and a2 refer to the thermal diffusivities, r and z are the space coordinates, and t is the time. 3.3 Surface Element Solution The surface element method which is presented in Chapter 2 is now employed to obtain the solution to the above problem. 45 3.3.1 Discretization of the Contact Area Due to the axisymmetry nature of the problem, the contact area can be divided into N annular surface elements (not necessarily equal) with each of these elements having inner and outer radii denoted by aJ._1 and aj, respectively (j=1,N and a0=0). See Fig. 3.2. In general the heat flux and the temperature vary across the interface (as well as with time), but they are approximated to be constant over each surface element (and time interval) and are specified at the points; r1 = 0 (3.3.1) a.-a._1 rJ = J__§J___; J = 2, N Since the two bodies are in perfect contact (hrw), the simplified form of solution given by (2.5.14) can be used. Using (2.5.8d) and (2.5.13) in (2.5.14) yields: T. (3.3.26) ll '6'" H l cl1 __ __ =_ M-l qM Mq1 - o1 '[iEI ¢ M+1-i qi] for M=2,3,... (3.3.2b) where g and T, are the influence matrix and the initial temperature vec- tor defined by (2.5.7a) and (2.5.7d), respectively. At each time step the system (3.3.2) can be solved for unknowns qlM’ qZM’ ..., qNM provided the influence functions, ¢kji's’ and'Ta are known. The ¢kji functions are given by; 2 m ¢ .. = z ¢ .. (3.3.3) k31 m=1 kJT 46 wherecpEji is the temperature rise at element k due to a unit step heat flux at element j of surface m, at time ti = i-At. It is the building block solution and can be evaluated from the solution for the problem of a semi-infinite body heated with a constant_heat flux over an annular area (see Fig. 3.3). The latter solution can be found from the known available solution for the case of semi-infinite body heated by a constant disk heat source [33]. See Fig. 3.4. 3.3.2 Evaluation of the Influence Functioncbrji Consider the geometry of body m shown in Fig. 3.5. Let the tempera- ture rise at element-k(r=rk) due to unit heat flux at the disk with radius a. and at time ti be denoted by; J m z m (3.3.4) ekji _ 6 (rk, aj, ti) By applying simple superposition one can show that: ¢E1i = 9:1, (3.3.5a) ' k=1,2,...N m _ m m .= ¢kji - ekji - ek,j-1,i J 2,3,...N (3.3.5b) which can be conveniently calculated utilizing, j-l j=1,2,...N m m m 6 .. = B .. - ¢ . (3.3.6) "31 "3‘ n=0 km k=1,2,...N 3.4' Temperatures at the Surface of‘a Semi-Infinite Body with 3 Disk Heat Source In reference [33] a solution for the local surface temperature his- tory for a semi-infinite body exposed to a circular heat source is provided 47 uniform heat qo Figure 3.3 Semi-infinite body heated over an annular-shaped region aj_]1 (3-4-3b) = ;,2- for r+=1 (3.4.3c) The functions K(-) and E(-) are the complete elliptic integrals of the first and second kinds, E. K(n) = [2 [I-n2 sin 2y] 1Nady (3.4.4a) o 50 1T 2' E(n) = I [1-n2s1n2y11/2dy (3.4.4b) 0 These functions are tabulated in [34] (where the definition is slightly different) and are available in computer libraries. Equation (3.4.1) provides a very accurate solution for large dimen- sionless times depending on the radial locations. It can extend down to dimensionless time t+:.04 for small r+'s, less than unity. At the early times the problem of semi-infinite body heated over a disk area can be approximated by the case of a semi-infinite cylinder insulated on the sides and heated over a disk area centered at the end. See Fig. 3.6. In reference [35], the temperature solution for the sur- face of the semi-infinite insulated cylinder exposed to a circular heat source is provided in dimensionless form which is valid for all radii and times. The solution is, + + + 2 t+ 1/2 ” T (r ,0,t ) = -——-(——) -2 2 b+2 " ‘1=1 (3.4.5a) erfc(th+1/2/b+)a (A]r+/b+)J (AT/6+) 1 0 1 1 1 + + 4. + 2 + 1(1" 9b ) [1,4009] where + + + + + co J (A.r /b )J (x./b ) I(r+,b+)52 2 0 l 1 ‘ (3.4.5b) i=1 [1130(11112 b = b/a (3.4.5C) 51 . . . + and the d1men51onless e1genvalues Xi are found from; 31(1:)=o (3.4.6) The functions erfc(o), J0(-), and J1(-) are the complementary error func- tion, and the Bessel functions of zero and first orders, respectively [34]. Except for extremely small dimensionless times (t+<.01) or large b+(b+>10) the evaluation of (3.4.5a) does not involve a large number of terms in the summation unless one tries to evaluate the steady state part, I(r+,b+), directly. Some simplified indirect procedures of evaluat- ing I(r+,b+) are provided in [35]. A new direct method is also developed by the author which evaluates this part more effectively. In this method the tail of the summation is replaced by an integral in an approxi- mate manner. This integral can simply be integrated over the correspond- ing limits. (See Appendix B). The evaluation of the explicit summation (3.4.5a) requires only less than 12 terms for t+/b+2 being greater than .33 and about 121 terms for extremely small values of t+/b+2(t+/b+22.0001). This is because erfc (n) decreases very rapidly with n and also because A: monotonically increases with i. The solutions given by (3.4.1) and (3.4.5a) can be used to evaluate etji (defined by (3.3.4)) for large and small dimensionless times, respectively. Then (3.3.6) and (3.3.3) can be employed to determine the influence functions, okji's. Once the influence functions and consequently the influence matrices, 31's, have been determined, the solutions to the heat fluxes qlM’ q2M’ ... qNM’ can be obtained by substituting the results into (3.3.2a) and (3.3.2b) for M=1 and M>1, respectively. 52 3.5 Thermal Constriction Resistance of the Two Solids In reference [9 ] the transient thermal constriction resistance is defined as "the difference between the average temperature of the contact area and the temperature far from the contact area divided by the total instantaneous heat flow through the contact area." Based on the above definition one can write TC(t)'T11 712']E(t) RC1(t)= w, Rc2(t)= W (3.5.1a,b) where TE(t) is the average temperature of the contact area, Qc(t) is the total heat flow through the contact area, and Rc1(t) and Rc2(t) are the thermal constriction resistances for body 1 and body 2, respec- tively. The total thermal constriction resistance for the two semi- infinite bodies can be determined by T. -T. Rc(t) = Rc1(t) + RC2(t) = {i163— V (3.5.2) The average contact area temperature, Tg, can be obtained by summing the products of the elemental temperature and the fraction of the total contact area occupied by the element. N T (tM) = Z T (Aj/Ac) (3.5.3) i=1 3” where TM is the temperature at the center of element j at time t = M-At and M A = n32, A. = n(a. - aj_1) (3.5.4a,b) 53 The total heat flow through the contact area, Qc’ can be determined by summing up the heat flows over all elements, N QC(tM) = .21 quAj , (3.5.5) J: Utilizing (3.5.3-3.5.5) into (3.5.1a,b) and (3.5.2) yields, 1 N AE'jEI TjMAj'Til RC1(tM) = N (3.5.6a) 1:1 “111% 1 N T. - ——- T. . 12 Ac jfl JMAJ RC2(tM) = N (3.5.6b) 3'51 quj T. -T. RC(tM) = %——‘-1— ‘ (3.5.6c) jfl quA1 3.6 Results and Discussions The two cases a and b mentioned earlier in subsection 3.1 have been solved. In each case the contact area has been divided into ten variable-spaced elements with smaller elements being closer to the edge of the contact area. The exact solutions given by (3.4.5a) and (3.4.1) were used to calculate the influence functions at "early" and "late" times, respectively. The values of the influence functions for dimensionless times .001, .1, and 1000 are listed in Tables 3.1 through 3.3 These Tables also represent the influence matrices. Notice that 54 TABLE 3.1 Values of the influence function, ¢k' , for semi-infinite body at . . . + 3” d1mens1onless t1me tM = .001. k/j 1 2 3 4 5 6 7 8 9 10 1 .0357 .0 .0 2 .001 .0356 .0001 .0 3 .0 .0012 .0332 .0014 .0 4 .0 .0012 .0332 .0013 .0 5 .0 .0012 .0332 .0013 .0 6 .0 .0012 .0332 .0013 .0001 7 .0 .0045 .0265 .0044 .0003 8 .0003 .0042 .0265 .0044 .0003 9 .0003 .0042 .0265 .0044 10 .0003 .0042 .0265 TABLE 3.2 Values of the influence function, ¢k3 M’ for semi- -infinite body at dimensionless time tM= k/j 1 2 3 4 5 6 7 8 9 10 1 .1887 .1665 .0750 .0697 20656 .0596 .0280 .0268 .0257 .0245 2 .0600 .2312 .0909 .0791 .0710 .0643 .0299 .0285 .0272 .0259 3 .0346 .1254 .1433 .1009 .0823 .0717 .0328 .0310 .0294 .0279 4 .0260 .0880 .0830 .1447 .1003 .0809 .0361 .0338 .0318 .0299 5 .0203 .0663 .0572 .0853 .1453 .0995 .0415 .0380 .0352 .0329 6 .0162 .0519 .0431 .0595 .0866 .1454 .0530 .0456 .0408 .0373 7 .0138 .0439 .0359 .0483 .0658 .0988 .0837 .0591 .0484 .0426 8 .0124 .0394 .0320 .0426 .0567 .0789 .0558 .0836 .0589 .0482 9 .0112 .0355 .0287 .0379 .0500 .0666 .0432 .0558 .0835 .0589 10 .0102 .0321 .0258 .0339 .0439 .0576 .0360 .0433 .0558 .0834 55 TABLE 3.3 Values of the influence function, oij, for semi-infinite body at + ‘l dimensionless time tr = 1000. x L; H ouooouowm-bwmi—n .1996 .0708 .0453 .0366 .0308 .0266 .0241 .0227 .0214 .0203 2 .1989 .2635 .1584 .1197 .0977 .0830 .0746 .0699 .0658 .0622 3. .0992 .1148 .1671 .1066 .0805 .0662 .0588 .0547 .0512 .0482 4 .0990 .1082 .1297 .1733 .1136 .0875 .0760 .0701 .0652 .0610 5 .0988 .1050 .1160 .1338 .1785 .1194 .0983 .0890 .0817 .0758 6 .0987 .1031 .1102 .1191 .1374 .1828 .1359 .1158 .1032 .0940 7 .0493 .0511 .0538 .0569 .0621 .0735 .1039 .0759 .0632 .0559 8 .0492 .0508 .0531 .0557 .0597 .0671 .0804 .1048 .0767 .0642 9 .0492 .0506 .0526 .0548 .0580 .0634 .0709 .0812 .1057 .0778 10 .0491 .0504 .0522 .0540 .0567 .0609 .0660 .0715 .0819 .1064 56 the diagonal terms of these matrices are greater than the off diagonal terms. The influence matrices were used in (3.3.2a,b) to find the un- known elemental heat fluxes. Equations (3.3.2a,b) were evaluated with three computer programs developed in the course of the project, RAN3, RAN2, and BBY9. Each used two routines called LINIV_and VQNLL_from the ISML Computer Library to perform the inversion and multiplication of the influence matrices. The results obtained from the surface ele- ment solution in this chapter, were compared to those given by the other investigators [24,26,29,3l,32,19] on the bases of thermal con- striction resistance and heat flux across the contact area. 9.2.1. The first problem examined is that of an isothermal disk on the surface of a semi-infinite body. With values of k=1, a=1, a=1, and Ti=2’ the problem was solved for elemental heat fluxes for various values of dimensionless time, from t+=.001 to t+=104. These results are tabulated in Table 3.4a, and 3.4b. Table 3.5 provides the normalized spatial variation of the surface heat flux at different times. Normalization obtained by dividing elemental values by the value of the centerline element which covers the area 0 5_r+ 5_.2. It can be seen that all the elemental heat fluxes approach the steady state condition and after dimensionless time about 20, the normalized heat flux distribution remains constant. This indicates that for t+>__20, the heat flux across the disk can be approximated by a pro- duct of a function of t+ and a function of r+. From the late time asymptotic solution given by Normington and Blackwell [24], it can be shown that the heat flux across the disk has a steady state distribu- +2)-1/2. tion of (1-r This implies that heat flux goes to infinity as TABLE 3.46 Results for elemental heat flux, qj, for an isothermal disk region on the surface of a semi-infinite body (t+ = at/az). + r = 0 .3 .45 .55 .65 t+ El. #1 4 .001 18.30 8.30 8.30 8.30 8.30 .005 8.184 8.184 8.184 8.185 8.185 .01 5.787 5.788 5.789 5.789 5.791 .05 2.588 2.596 2.613 2.645 2.718 .1 1.846 1.871 1.914 1.970 2.072 .5 1.031 1.075 1.141 1.210 1.320 1. .8979 .9395 1.003 1.067 1.170 10. .7184 .7534 .8062 .8598 .9459 100. .6669 .6995 .7486 .7986 .8787 1000. .6510 .6823 .7308 .7796 .8578 10000. .6460 .6776 .7252 .7736 .8513 TABLE 3.4b Results for elemental heat flux, qj, for an isothermal disk region on the surface of a semi-infinite body (t+ = at/az). r+ = .75 .824 .875 .925 .975 t+ El. #6 7 8 9 10 .001 18 32 18.32 18.32 18.33 25.75 .005 8.195 8.237 8.424 8.851 16.19 .01 5.830 5.966 6.270 6.864 13.51 .05 2.903 3.198 3.570 4.204 9.164 .1 2.284 2.586 2.939 3.526 7.875 .5 1.519 1.781 2.067 2.532 5.824 1. 1.354 1.595 1.856 2.281 5.269 10. 1.098 1.297 1.513 1.864 4.318 100. 1.021 1.206 1.407 1.733 4.015 1000. .9963 1.177 1.373 1.692 3.920 10000. .9886 1.168 . 363 1.679 3.890 58 TABLE 3.5 Normalized spatial variation of the surface heat flux at different times for the geometry of Case a, see Fig. 3.1. (q+ = q/chvs r+ = r/al- + r t+ 0 .3 .45 .55 .65 .75 .825 .875 .925 .975 .001 1. 1. 1. 1 1. 1. 1. 1. 1. 1.41 .01 1. 1. 1. 1. 1. 1.01 1.03 1.08 1.19 2.33 .1 1. 1.01 1.04 1.07 1.12 1.24 1.40 1.59 1.41 4.22 1. 1. 1.05 1.12 1.19 1.30 1.51 1.78 2.07 2.54 5.87 10. 1. 1.05 1.12 1.20 1.32 1.53 1.81 2.11 2.60 6.01 20. 1. 1.05 1.12 1.20 1.32 1.53 1.81 2.11 2.60 6.02 w 1. 1.05 1 12 1.20 1.32 1.53 1 81 2.11 2 60 6.02 Averaged steady state values from 1/J1-r+2 based on the elemental area m 1.05 1.12 1.20 1.32 1.52 1.77 2.08 2.67 6.35 TABLE 3.6 Normalized area averaged interface heat flux for the geometry of Case a, see Fig. 3.1. 'q'+= E/ 300. Beck & Keltner [19] Marder & Keltner [29] t TSEM Eq. (56) Eq. (22) s = 2 s = 5 .001 15.057 12.357 14.522 13.088 8.646 .002 10.872 9.030 10.434 10.408 7.950 .005 7.169 6.079 6.792 7.063 6.485 .01 5.310 4.591 4.969 5.238 5.145 .02 3.998 3.539 3.688 3.948 3.943 .05 2.843 2.606 2.573 2.813 2.813 .1 2.268 2.136 2.038 2.248 2.248 .2 1.871 1.803 1.702 1.858 1.858 .5 1.533 1.508 1.530 1.515 1.526 1. 1.371 1.359 1.667 1.297 1.366 2. 1.259 1.254 2.334 1.104 1.257 5. 1.163 1.161 1.005 1.155 10. 1.114 1.113 1. 1.082 20. 1.080 1.080 1. 1.023 100. 1.036 1.036 1. 1000. 1.011 1.011 10000. 1.004 1.004 m 1.0 1. 59 r+ goes to 1. The steady state values of the elemental heat flux given in Table 3.5 are compared with the corresponding values obtained from the asymptotic solution, (1-r+2)’1/2 , based on the elemental area. There is good agreement between the two solutions for all ele- ments but the outermost element, which indicates an error of 6%. Fig. 3.7 shows the normalized flux distribution across the disk at several times. It can be seen that the region of uniform heat flux shrinks as t+ goes to infinity, which implies that the assumption of uniform heat flux across the disk (made by Beck and Keltner [19] in their q- based solution) is only good at early times. Fig. 3.8 and 3.9 illustrate comparisons of area averaged inter- face heat flux]r calculations performed using SEM with the results ob- tained from the two solutions given by Equations (22) and (56) in reference [19] and the two other solutions given in reference [29] for values of s=2 and s=5. In Fig. 3.8 normalized values of averaged heat flux (6(t)/6(m)) is plotted versus dimensionless time t+. While Fig. 3.9 shows the percent of deviation in the values of the heat flux obtained by the above-mentioned alternative solution from that found by SEM solution. Values of fiTt)/fi(m) as a function of time are also tabulated in Table 3.6 for the same five solutions. The Eq. (22) solution [19] is based on a heat flux form of Duhamel's integral. This solution uses approximate expressions for influence functions at early dimensionless times. It closely matches the SEM solution up to time t+ =.005 (less than 5% deviation, see Fig. 3.9). + The area averaged interface heat flux, ETt), was obtained by summing the products of the elemental heat flux and the fraction oftfiwztotal interface area occupied by the element. 60 q/qE q: + O 1 A J 1 l 1 4 l 1 0 .5 l r+-r/a Figure 3.7 Normalized heat flux distribution across the interface at various times for the geometry of Case a, see Fig. 3.1. 61 SEM _ . . Eq.(22),[l9] ---Eq.(56),[19] O O ch,[29]' ° ° s=5,[29] Figure 3.8 Normalized area averaged interface heat flux versus dimension- less time for the geometry of Case a, see Fig. 3.l. 62 20 T T I ' I l =5,[29] /5 15 - ‘- Eq.(56).[191 s=2.[291 Eq.322) [19 10 -- .- Percent of Error 0 1 J 1 J 10'3 10'2 10" 10° 101 102 10 t+= at :7 L. w 104 Figure 3.9 Percent of relative error in the area averaged interface heat flux with respect to SEM solution as a function of time for Case a, see Fig. 3.l. 63 Between t+ = .005 to t+ = .5 the deviation increases up to 10%, and for t+ > .5 the solution ceases being useful. This solution is only appro- priate for the small dimensionless times. The Eq. (56) solution is in good agreement with SEM solution for mid-to-late times (t+ > .1). It matches exactly the SEM solution for times greater than 20 (see Fig. 3.8). The percent of deviation in the surface heat flux starts with the maximum value of 17.9 at t+ = .001 and decreases as t+ becomes larger, it reaches a value of 5.8 at t+ = .1 and approaches to zero at t+ z 20. See Fig. 3.9. The validity of the s = 2 solution is in question for very short times t+ :_.001 (13% deviation for t+ z .001 shown in Fig. 3.9). However, it is in good agreement with the SEM solution for short and intermediate times between t+ 2 .002 to t+ z 1. The s = 5 solution remains good until about t+ 2 10. Even though reference [29] states that the s = 5 solution is very good down to nearly t+ = .001, the results from Fig. 3.8, 3.9 and Table 3.6 show that it is in poor agreement with SEM solution for t+ < .01 (Table 3.6 indicates a difference of 42.6% at t+ = .001). Because of the effect of the finite artificial boundary location assumed in [29] (the assumption made in order to use the sepa- ration of variables technique [29]) the s = 2 and s = 5 solutions both reach the steady state value too early; at t+ 2 10 and at t+ 2 100, respectively. 9355.9 The composite Case b is considered with the glass being the material of one body and copper of the other. The glass-copper combination has also been investigated in references [31,32,19]. The thermal conducti- vities, k's, are 1.03 and 381 W/m-k and the thermal diffusivities, 64 0'5 are .06 x 10'5 mZ/s and 13.2 x 10'5 mz/s for glass and copper, respectively. The initial temperature of glass is T1.1 = 0, and that of copper is T1.2 = 2 K. With value of a = 1 m, the problem was solved for the elemental surface heat fluxes and temperatures. Results for the sur- face heat fluxes are displayed in Table 3.7. Column one in this table contains dimensionless times, t+, based on the lower thermal diffu- sivity, 01 < 02, ranging between .001 and 10000. The results for spatial variation of the surface heat flux, given in Table 3.5, for Case a, is almost the same for Case b. That means after dimensionless time about 20 the normalized heat flux distribu- tion across the interface remains constant. In Fig. 3.10, the normalized interface temperature distribution is plotted versus dimensionless radius, r+, at several times. Normalization is obtained with respect to the initial temperature of body 2 (copper). It can be seen for t+ = .1 that the interface temperature distribution is almost uniform. This supports the validity of the T-based solution given by Beck and Keltner for the late times. Notice that the region of uniform temperature shrinks at t+ goes to zero which implies that the T-based solution is not appropriate for short times. Figures 3.11 and 3.12 compare the results obtained for Case b with those found from the T-based and the q-based solutions of reference [19], based on the area averaged interface heat flux and the area .1. averaged interface temperature , respectively. Again it can be seen that the q-based solution given in [19] is in good agreement at early T The area averaged interface temperature, T(t), is obtained by summing the products of the elemental temperature and the fraction of the total interface area occupied by the element. TABLE 3.7a Results for elemental heat flux, qj, for a disk-shaped interface of two semi-infinite bodies, one of copper and the other of glass (t+=alt/a2)+. 10. 100. 1000. .45 r = 0 .3 .55 .65 El. #1 2 3 4 5 36.66 36.71 36.78 36.83 36.89 11.84 11.85 11.85 11.85 11.85 3.791 3.842 3.931 4.044 4.254 1.844 1.929 2.059 2.191 2.403 1.475 1.547 1.656 1.766 1.943 1.370 1.437 1g538 1.641 1.805 1.337 1.403 1.501 1.602 1.762 .1. TABLE 3.7b 01 = Thermal diffusivity of glass (body one). Results for elemental heat flux, qj, for a disk-shaped interface of two semi-infinite bodies, one of copper and the other of glass (t+=alt/a2)+. 10. 100. 1000. r+=.75 .825 .876 .925 .975 El. #6 7 8 9 10 36.98 37.05 37.07 36.49 52.16 11.93 12.21 12.83 14.03 27.58 4.687 5.307 6.030 7.231 16.15 2.780 3.274 3.811 4.684 10.82 2.256 2.665 3.108 3.829 8.870 2.096 2.477. 2.889 3.560 8.248 2.047 2.418 2.821 3.475 8.053 .1. 01 = Thermal diffusivity of glass (body one). 66 100 ' ' '11: 1 '1 99- - 01 / 961— .- Percent of Interface Temperature Change + r =r/a Figure 3.l0 Normalized interface temperature distribution for the geometry of Case b, see Fig. l.2a.(Normalized with respect to initial temperature of body 2, copper) 67 15 1 1 1 I l 1 10 Eq.(56).[191 \ E '3 '3 SEM |*: S!— — EQ-(221.[]9] 0 l J 1 1 j 3 - - - l 10 3 10 2 10 ‘ 100 101 102 103 10L 4, 0111: t=-—— 2 Figure 3.11 Dimensionless area averaged interface heat flux versus dimen- sionless time for the geometry of Case b, see Fig. l.Za, (Normalized with respect to the steady state value of q) 68 1.1 . T 1 . 1 T l,— EQ-(28).[l9] SEM 1.0- . /_Eq.(57),[l9] /_ 102 10 10 Figure 3.l2 Dimensionless area averaged interface temperature versus dimensionless time for the geometry of Case b, see Fig. 1.2a. (Normalized with respect to the steady state value of T3 69 times while the T-based solution is suitable for mid-to-late times, t+ > .1. Comparison Based on Thermal Constriction Resistance The SEM solution is also compared with the other available solu- tions on the bases of the dimensionless thermal constriction resis- tance across the contact area for both Cases a and b. Tables 3.8 and 3.9 provide these results. Table 3.8 is for Case a, two identical semi-infinite bodies, and Table 3.9 is for Case b, with copper being the material of one body and glass that of the other. The first column in each table is the dimensionless time, t+, which extends over many decades. The results from the finite difference solution of Schneider et al [31] are provided in the second columns which are least accurate at the early times and most accurate at the late times. The third columns come from the exact solution given by Sadhal [32] which is claimed to be very accurate at late times, t+ > 10. The fourth columns are for the T-based solution given by Beck and Keltner[19], which is most appropriate for large times but is remarkably accurate down to t+ 2 .1. Column 5 in each table is for the q-based solution [19] and is accurate at early times. The results obtained by the SEM solution are displayed in the 6th columns. The last (7th) column of Table 3.9 (for case a) was calculated from the one-dimensional approximate solu- tion given by Keltner [26]. which closely matches SEM solution at the short times and retains good accuracy for the mid-to-late times. Values of R: as a function of time are also plotted in Fig. 3.13 for the same solutions of Case b. Fig. 3.14 shows the percent of error between the above-mentioned 70 TABLE 3.8 Results for dimensionless constriction resistance, for an isothermal disk region on the surface of a semi-infinite body. R: = Rc-a-k. SSY Sadhal Beck & Keltner [19] Keltner t+ [31] [32] Eq. (56) Eq. (22) SEM [26] .001 .0386 .0202 .0172 .0166 .0162 .002 .0409 .0277 .0240 .0230 .0223 .005 .0463 .0411 .0368 .0349 .0340 .01 .0532 .0544 .0503 .0471 .0473 .02 .0637 .0706 .0678 .0625 .0641 .05 .0851 4.8988 .0959 .0972 .0879 .0914 .1 .1074 .2029 .1171 .1226 .1102 .1142 .2 .1336 .1685 .1386 .1468 .1336 .1382 .5 .1695 .1752 .1658 .1634 .1631 .1685 1. .1933 .1879 .1839 .1500 .1824 .1895 2. .2120 .2010 .1994 .1071 .1984 .2074 10. .2368 .2247 .2245 .2242 .2347 100. .2475 .2413 .2413 .2414 .2477 1000. .2495 .2472 .2472 .2472 10000. .2499 .2491 .2491 .2491 w .2500 .2500 .2500 .2500 74 TABLE 3.9 Results for dimensionless constriction resistance for a disk-shaped interface of two semi-infinite bodies, one of copper and the other of + glass. Rc - RC/Rc(w). SSY Sadhal Beck & Keltner [19] t+ [31] [32] Eq. (56) Eq. (56) SEM .001 .1549 .0820 .0701 .0676 .002 .1645 .1118 .0968 .0932 .005 .1870 .1655 .1462 .1404 .01 .2158 .2187 .1961 .1890 .02 .2595 .2833 .2570 .2507 .05 .3474 .3844 .3480 .3523 .1 .4382 .4688 .4101 .4413 .2 .5442 .5552 .4418 .5349 .5 .6869 .6637 .3857 .6610 1. .7805 .7564 .7362 .2539 .7300 10. .9495 .8992 .8982 .8973 100. .9903 .9654 .9654 .9654 1000. .9982 .9888 .9888 .9889 10000. .9997 .9964 .9964 .9964 0 72 SEM -- - EQ-(551.U9] .0 <0 quzzhiun o o [32] . o o [3]] l J 4 J J l l 19’3 10'2 10'1 10° 101 102 103 10 omit t: :2" Figure 3.l3 Normalized thermal constriction resistance across the inter- face versus time for the geometry of Case b. 73 10— .— L- 2 - 1.: ‘ . , . F 5 .‘ , l l . 0 l - [32] - 2- [31] . 3- EQ-(561.[19] 1 _ 4- Eq.(22),[19] . 5- [25] -5 1 J J 1 1 4 10'3 10’2 10‘] 10° 101 102 103 10‘ t _ at :2 Figure 3.14 P rcent of relative error in thermal constriction resistance, Rc’ with respect to SEM solution as function of time,Case a. 74 solutions [26,29,3l,32,l9] and the present solution for Case a. It can be seen that the solution given by Keltner [26] is in the best agree- ment, for the total time within approximately 5%. The two approximate solutions given by Beck and Keltner are the next closest. They provide remarkably accurate results for entire region except between t+ = .01 to t+ = .1 where the error is slightly high (6 to 10%). In particular the T-based solution is very accurate at the late times. It closely matches the present solution and the exact solution given by Sadhal [32] for times t+ > 2 (less than .5% error). The F0 solution given by Schneider et al [3T], is good for mid-to-late times but is in poor agreement at early times, less than .01. It shows 13% error at t+ = .01 which increases to 132.5% at t+ = .001. A study of the above results indicates that the multinode surface element solution presented here is the most accurate over the entire time region. CHAPTER 4 INTRINSIC THERMOCOUPLE PROBLEM 4.1 Introduction Intrinsic thermocouples are widely used for measuring rapid tran- sient surface temperatures of conducting solids. Such measurements are important in studies of nuclear reactors, laser heating, re-entry vehicle heating, and other applications where the response of the sur- face temperature is of interest. An intrinsic thermocouple as defined in [27] is "one in which the material whose temperature is to be measured (called substrate) forms part of the thermocouple circuit." The most common geometry for the intrinsic thermocouple problem is a semi-infinite cylinder (called wire) attached perpendicularly to a semi-infinite solid (substrate) as shown in Fig. 1.1a. Heat is conducted along the thermocouple wire away from the junc- tion if the substrate's temperature increases. This heat transfer between the substrate and the wire can introduce a significant error in measurement of the interface temperature. For single step change at t = 0, the error is maximum at the initial moment and decreases with time. In this chapter the transient thermal behavior of the interface between the thermocouple wire and the substrate due to a step change in the substrate temperature is analyzed. The surface element method is employed to obtain the transient solution for the interface heat fluxes and temperatures. The results are compared with those obtained by other investigators for the same problem. The wire and the sub- strate are considered to be in perfect contact. 75 76 Previous Work Intrinsic thermocouple problem was first studied by Burnett [36] in 1961. He developed an approximate analytical solution for a single wire attached normally to the rear (insulated) surface of a thin slab which was exposed to a constant heat flux on its front face. The wire was infinitely long and insulated on the sides. See Fig. 4.1. The effect of the wire was simulated by a uniform disk-shaped sink. The sink strength was assumed to be equal to the flux into a semi-infinite solid which had a surface temperature equal to that of the insulated face of the slab in absence of the wire. However, this temperature was higher than that at the actual slab-wire junction and consequently, the solution provided a conservative estimate for the error in the tem- perature measurements. Larson [37], and Larson and Nelson [38] im- proved Burnett's solution by applying a correction factor which pro- vides a more accurate value for the strength of the sink. Their solu— tion, like the Burnett solution, is only valid for the dimensionless times (with respect to the substrate) greater than one. In 1967, Henning and Parker [39] studied the transient thermal response of an intrinsic thermocouple due to a step change in the sub- strate temperature. They modeled the system as a semi-infinite cylin- der attached normal to the surface of a semi-infinite body (called the idealized model). See Fig. 4.2. Perfect thermal contact was assumed between the wire and the substrate and all thermal properties were considered temperature independent and constant (linear). Also it was assumed that the heat loss from the wire to the surrounding was negligible. To simplify the analysis, they introduced a hemispherical region shown in Fig. 4.2, between r = R and x = 0, with infinity thermal 77 F-—::E;”,, substrate(thin slab) r 2 6077//////' “311111111 A $7/7///( I ” wire ;; x w uniform q Figure 4.l The geometry for Burnett's problem. ' substrate Figure 4.2 Henning and Parker's idealized geometry for an intrinsic thermocouple problem. 78 ¢:onductivity and no heat capacity. This assumption resulted in the temperature distributions in both the wire and the substrate being functions of only one space dimension (axial in the wire and radial in the substrate) and time. Utilizing the above assumptions they derived analytical solutions for transient temperature distributions in the semi-infinite substrate and the semi-infinite wire by using the method of separation of variables. To investigate the validity of the analysis, they further per- formed a series of experiments by means of the capacitor pulse-heating method [39,40]. The experimental data were compared with the analy- tical results based on the response-time+ of the thermocouples. In order to improve agreement between analytical and experimental results, Henning and Parker introduced an empirical correction factor, "G", which is related to the shape of the isotherms in the substrate. This modification improved the analytical solution for the late times. How- ever, (as stated by other authors too [27,4l]) errors up to 20% were still observed for early-to—mid times (t+ < 4). Giedt and Nunn [42], further improved Henning and Parker's solution by applying a modifica- tion which provided a more accurate value for the early time response of a thermocouple. They stated that the early time solution of an in- trinsic thermocouple can be approximated by using the solution for two semi-infinite bodies (initially at different temperatures) that are suddenly brought together. The response-time of an intrinsic thermocouple is the time that the junction temperature reaches 95% of the initial step change in the substrate temperature [27,39]. 79 In 1971, Bickle [43] and Bickle and Keltner [44] examined the problem by using a combined experimental-numerical method. They em- ployed the "deconvolution" technique developed in [43,44] to predict the thermocouple response due to an arbitrary substrate temperature variation. 1 In 1973, Keltner [27] studied the response of a single wire in- trinsic thermocouple for the case where the substrate undergoes a step temperature change utilizing finite difference and quasi-coupling methods. He was first to employ the finite-difference procedure to solve the problem for both one-dimensional and two-dimensional models. He examined four finite difference models. In his first model (Two-D), Keltner used spherical coordinates in the substrate and cylindrical coordinates in the wire. Because of an unexpected inflection point which occurred in the interface temperature and heat flux distributions, he dropped this model in favor of the next two models which used cylindrical coordinates in both the wire and the substrate [29]. Models 11 and III were both axisymmetric and two-dimensional as shown in Fig- ures 4.3 and 4.4. Model III is essentially the same as model II ex- cept in the node size and the total number of nodes. Because of the finer node structure, model III provided better results than model 11, specially for the region near the outer radius of the interface between the substrate and the wire. However, it was found that the model III finite difference grid is still not fine enough to precisely model the +corner region near r+ = 1 [27]. The fourth FD model was one dimensional in oblate spheroidal coordinates for the substrate and in cartesian (:oordinates for the wire. Because of the deficiencies involved in the finite difference SSolution such as; the effort in setting up large grids, large computer mnmu mucmgmmmg cw cmgmuwmcoo HH vaoz 202 cmms mucmgmw$wu wumcwu m.¢ wgzmwu 8] Hum; mucwemmme cw uwewuwmcou HHH quoz Low :mme mucmcmmwwv waver; ¢.¢ mgsuvu 222,... mmmmmm m 82 expense and restricted dimensionless time that could be covered, Keltner developed an alternative method that he called "quasi-coupling". In this procedure-he assumed one-dimensional temperature distributions in both the substrate and the wire. To provide a uniform heat flux into the end of the cylinder (wire), a disk-shaped infinitesimally thin region having zero heat capacity and infinite thermal conductivity was inserted between the substrate and the wire. (In other words, the interface temperatures were assumed to be only a function of time; the interface temperature could be considered to be the average across the interface.) Utilizing this method Keltner obtained remarkably good agreement with his finite difference solutions. The quasi-coupled solutions unlike his finite difference solutions were simple and inexpensive to obtain. In 1976, Shewen [41] examined the problem by utilizing a two- dimensional finite difference model in oblate spheroidal coordi- nates for the substrate and in cylindrical coordinates for the wire. The two dimensionality of his solution allowed a more detailed study of the temperature variations along the interface, and, consequently. provided more accurate results particularly at the early times. Keltner and Beck [18] were the first to employ the surface element method to solve the problem. They have considered only one element across the interface and solved the problem analytically by utilizing Laplace transform techniques. Some approximate solutions based on the heat flux and temperature kernel forms of Duhamel's integral were presented for early and late times. All of the above-mentioned solution methods (except the two- dimensional finite difference solutions given by Keltner and Shewen) ignored the two-dimensionality of the problem which is especially im- portant at early times. Even though the two-dimensional finite 83 difference solutions of Keltner and Shewen provided more accurate re- sults, they are not entirely satisfactory. They have some difficulties due to the necessity of setting up extremely fine grids near the inter— face and many large grids further from the interface. This in turn increases the number of numerical calculations and computer time. 4.2 Statement of the Problem The idealized geometry for the intrinsic thermocouple problem is shown in Fig. 1.1a. The semi-infinite cylinder referred to as region 1 and the semi-infinite substrate as region 2. At zero time the entire substrate undergoes a step change of temperature to T1.2 while the wire is at its initial temperature of T11. The describing differential equations, initial and boundary conditions are the same as those for two semi-infinite bodies attached over a circular area given in Section 3.2 of Chapter 3, except for the initial and some boundary conditions of region 1, given by (3.2.2a), and (3.2.3a) and (3.2.6a), respectively. These conditions for the wire (insulated cylinder) are changed to, T1 = T1.1 for t = 0, Ogrga, 230 (4.2.1) and T1 = T1.1 for t > 0, Ogrga as 2 + 2 (4.2.2) 3T1 ‘T‘ = ° [0" t > 0+ r = a, 2:0 (4.2.3) 4.3 Solution The solution to the intrinsic thermocouple problem can be found using the surface element method in the same manner as presented in 84 Chapter 3. The contact area between the wire and the substrate is divided into ten annular variable-spaced elements with smaller elements being closer to the corner region. See Fig. 4.5. The temperature and the heat flux are assumed uniform over each element. Equations given in Sections 3.3 and 3.4 are valid and can be used for the intrin- sic thermocouple problem. The only exception is that in Chapter 3 both bodies were semi-infinite and thus the required influence functions were both evaluated from the solution of a semi-infinite geometry heated by a constant heat flux over an annular region on its surface. In the intrinsic thermocouple problem only one of the bodies has semi-infinite geometry (substrate); the other one is a semi-infinite cylinder (wire), 1 2 (bkjl ‘7‘ ¢ij (4-3-1) To evaluate the influence function for the wire the exact solution for the case of a semi-infinite insulated cylinder heated by a constant heat flux over a disk area centered at the end, given by (3.4.5a), can be utilized in conjunction with (3.3.5). The problem is analyzed for a chromel substrate and an alumel wire. The chromel and alumel combination has also been investigated in references [18,27,4l]. The thermal conductivities, k's, are 19.21 and 29.76 w/m-K, and the thermal diffusivities, a'S, are .492 x 10'5 and 5 .663 x 10' mZ/s for chromel and alumel, respectively. The initial temperature between the substrate and the wire is, Ti2 - T 1K. i1 = The problem is solved for the elemental surface temperatures and heat fluxes for various values of dimensionless time, ranging between + 4 t = .001 to t+ = 10 . The dimensionless time is based on the thermal diffusivity of the substrate (body 2), 85 insulated Figure 4.5 Discretization of the interface between the semi-infinite thermocouple wire and the semi-infinite substrate. 86 t = ;:_ (4.3.2) The results are compared with those obtained from the Henning and Parker's solution [39], Keltner's model 111 finite difference solution [27], Shewen's finite difference solution [41], and Keltner and Beck's single node surface element solution [l8] on the bases of the area averaged interface temperature. 4.4 Results and Discussion The results of the surface element solution are presented in terms of the temperature and heat flux distributions across the interface, and the area averaged interface temperature and heat flux as functions of time. Fig. 4.6 shows the normalized spatial variation of the inter- face temperature at different dimensionless times. These results are also provided in Table 4.1. Normalization is obtained by subtracting 11’ It can be seen from Fig. 4.6 that the temperature gradient in the region T1.1 from elemental values, Tj’ and dividing the results by T1.2 - T near the corner of the interface is very large at early time and ap- proaches to zero as t+ goes to infinity. This justifies the two dimen- sionality of the problem and also indicates that more accurate studies are required in this region, especially at the early times. Keltner in his model III finite difference solution used more than 400 nodes in the semi-infinite body to represent the corner region effectively. He found that even his fine grid structure introduces some errors. Employing oblate spheroidal coordinates Shewen [4l] ob- tained approximately the same results with a smaller number of nodes (120) and consequently less computational effort than that used by 87 Hwhnmwp u ah pruwh + maaa. maaa. maaa. maaa. maaa. maaa. maaa. maaa. maaa. maaa. 50H meaa. Haaa. oaaa. amaa. mmaa. mmaa. mmaa. emaa. mmaa. Nmaa. .ooooH omma. efima. ona. moma. Noma. mama. vama. HaNa. mmma. mmma. .oooH aqua. mmva. aoca. mama. emma. amma. mmma. avma. oema. mmma. .ooH mfiem. ommm. aHmm. mmmm. ccmm. Homm. “mam. Heam. mHHm. Haom. .cH mean. ummn. Hams. mmun. HHNN. Nmmm. momx. Hnmu. Nmmn. momm. .m mmmm. mamm. ammm. mmvm. ommm. ammm. NmHm. cmom. HHom. mmam. .H mmmm. momm. mmom. mmam. mmmm. Namm. Nmmm. maem. Noam. Hmmm. m. mHom. nuam. comm. mmqm. mamm. vfifim. mnae. nmmv. ease. amme. N. mmmm. mmmm. ammm. mmfim. amav. mmme. Name. mama. aeve. mmme. fl. ammm. nmmm. maom. mHac. mHNe. mama. mmev. mmmv. Name. vamv. mo. «Nam. mmom. mane. omme. mmev. name. meme. flame. mmme. mmme. No. ommm. Hmmv. mama. mavc. mmmw. mama. mmmc. mmmw. mmme. mmmv. Ho. umfim. cmmv. omqv. mmme. mama. mmmw. mmme. mmma. mmwe. mmmv. moo. mmav. mmva. NHme. amwv. mmma. mmmv. mmme. mmmv. mmma. mmmv. moo. mane. mmmv. mmmv. mmme. mmmw. mmmw. mmme. mmmv. mmme. mmmv. Hoo. oH a m m m m w m N a» .Pm +p mma. mma. mum. mmm. mu. mm. mm. me. m. o. +2 .Nm\u~a u +p we?» mmmrcowmcmewc mo cowpocsw 6 mm .mp .wsspmgmaEmu mummemucw Faucmsmpw mo mm3~m> umNPFmELoz Hé WEE. 88 tflnoa 20. be 1 + = (Tc’Ti11/(T12’Ti1) \J I L U "’ .6 , " .5 ° .5.. - l J . l .4.. 0 t - t+= 2%. SEM Solution a 0 FDM Solution[27] 0 FDM Solution[4l] .3i I 1 J 1 l 1 1 1 1 0 .5 l. r+=r/a Figure 4.6 Normalized interface temperature distribution for intrinsic thermocouple with chromel half-space and alumel wire. 89 Keltner. Their results are also shown in Fig. 4.6 for time t+ = .1. It can be seen that the agreement between the above finite difference solutions and the present surface element solution is very good. Notice that the number of nodes used in the SEM solution is only ten, 12 times less than that used by Shewen, and about 40 times less than that used by Keltner. When the substrate initially undergoes a step change of tempera- ture, there is an instantaneous change to a common interface tempera- ture which depends upon the thermal properties of the substrate and .1.. the wire . The normalized value of this initial interface temperature++ is given by [18], T. -T. + 1c 11 -1 T. =-——17——-= (1+3) (4.4.1) 'C T12 T11 where B = [kIDICI/k202C2]1/2 ‘ (4.4.2) For a chromel substrate and an alumel wire the value of TTc is equal to .4285. After the initial moment the effect of the edge starts to penetrate toward the center of the contact area. This is evident from the results given in Table 4.1. It can be seen that at t+ = .001, the interface temperature is almost at its initial value (.4285) for + Notice that at the initial moment there is no spatial variation in the interface temperature. ++ This temperature corresponds to the surface temperature of two semi- infinite bodies initially at different temperatures suddenly brought into the perfect contact. 90 r+ < .825. The only part of the contact area which is disturbed by the edge effect is between r+ - .825 and r+ = 1. The edge effect pene- trates further toward the center as t+ becomes larger, reaching r+ = .45 at t+ = .01 and covering the entire area for t+ > .02. The difference between the centerline and the outermost element temperatures is about 10.7% at time t+ = .001 and decreases as t+ becomes larger. For t+ > 20 the interface temperature distribution becomes almost uniform which indicates that the one-dimensional approximate solution given by Henning and Parker and the T-based solution of Keltner and Beck [18] are appro- priate for the late times. Fig. 4.7 shows the normalized flux distribution across the inter- face, q+(t), at several times. It can be seen that the region of uni- form heat flux shrinks as t+ increases. After dimensionless time about 20 the q+ distribution remains constant and a universal curve is obtained. Fig. 4.8 shows the normalized area averaged interface temperature versus the dimensionless time. Results for the model 111 finite dif- ference solution of reference [27], and the finite difference solution of reference [41] as well as the SEM solution are presented. Table 4.2 provides the results of the above-mentioned three solutions along with the results from the approximate analytical solution of reference [39], and the T—based and the q-based solutions of [TB]. The first column in this table is the dimensionless time which extends from t+ = .001 to t+ = 500. The results of the finite difference solutions of ref- erences [27] and [41] are given in the second and the third columns, respectively. The fourth column is evaluated from equation (22) given by Henning and Parker, which is only good for the late time, t+ > 20. The early and the late times results of the T-based solution and the 91 3 I 1 1 j 1 1 l 1 1 2 P- -I 01 -E c- +" t+= . 0] c t+=l t*=.2 t+_>_20 / J 1 d :/ Z .- t+=.001 4 l 1 l l 1 4 4 1 0 . 5 l r =r/a Figure 4.7 Normalized heat flux distribution across the interface at various values of time for intrinsic thermocouple. 92 fi I 1 l I 1 .9- - 5 .8- - T N PP S 15" 7- :0 II |+’_u ,6. _ .5. _ 0 SEM Solution ‘ o FDM Solution[27] 4” ' FDM Solution[41] 3 1 1 J 1 I 1 10'3 10‘2 10" 10° 10' 102 103 104 + azt t "—2 a Figure 4.8 Normalized area averaged interface temperature histories for intrinsic thermocoupie.problem. 93 q-based solution of Keltner and Beck, are displayed in the next four columns. The last column represents the present SEM solution. As can be seen from Fig. 4.8 and Table 4.2, there is very good agreement be- tween the finite difference solutions and the present solution for the time range covered. For a given time between t+ = .01 and 10, the SEM values are between the two finite difference solutions. The F0 solu- tion of reference [41] is more accurate at early times than that of reference [27]. This is because of that the oblate spheroidal coordi- nates used in the former model describes the geometry more closely than the cylindrical coordinates used in the latter model. As mentioned earlier, however, both solutions have difficulties regarding the com- putational effort and cost, particularly for the early times, t+ < .01. The T-based and the q-based single node surface element solutions given by Keltner and Beck [18] are convenient in that the mathematics is not difficult and the expressions are simple to evaluate. Each solution provides two expressions; one for early times and the other for late times. The Q-based solution is more appropriate for the early times. It approaches the exact solution (.4285) as t+ + 0, and closely matches the SEM solution up to dimensionless time t+ = .1. It also provides relatively good results for the late times t+ > 10. The T-based solution does not approach the exact solution as t+ + 0, and consequently it is less accurate than the q-based solution for the early times. Because of the constant interface temperature assumption, however, it yields very good results for the late times of t+ > .1. It should be noted that neither the T-based solution nor the q-based solution is solely suitable for the complete time domain. However, a combination of the early-time q-based solution (equation (41a) of [18]) and the late-time T-based solution (equation (33a) of [18]) provides TABLE 4.2 94 Normalized area averaged interface temperature histories for chromel substrate and alumel wire. + .001 .002 .005 .01 .05 10. 50. 100. 200. 500. T-based Solution [18] Solution q-based FD Solutions Henning and early late early late [27] [41] Parker [39] time time time time SEM .6084 .4489 .4342 .4335 .4421 .6118 .4500 .4366 .4364 .4480 .6185 .4521 .4413 .4422 .4402 .4510 .6257 .4545 .4467 .4488 .4505 .4599 .6356 .4546 .4581 .4700 .4782 .6540 .4709 .4765 .4916 .4991 .6731 .4907 .4904 .4973 .5215 .5283 .6972 .5280 .5263 .5770 .5826 .7373 .5910 .5805 .6302 .6338 .7729 .6452 .6328 .6896 .6921 .8109 .7042 .6915 .7688 .7714 .8602 .7810 .7700 .8202 .8246 .8933 .8327 .8091 .8236 .8694 .9207 .8757 .8614 .8687 .9139 .9482 .9186 .9108 .9137 .9382 .9629 .9417 .9365 .9381 .9735 .9585 .9550 .9559 .9832 .9737 .9715 .9719 95 very good results over the entire time range. These two solutions match very closely at dimensionless time t+ = .1. Fig. 4.9 shows the percent of error between the above-mentioned solutions [27,41,39,18], and the present solution. It can be seen that the finite difference solution given by Shewen [41] is in the best agreement for t+ > .002 within approximately 1%. The model III finite differerence solution of Keltner has a deviation about 2% at t+ - .01 which decreases to less than 1% at late times. The solution given by Henning and Parker is good for late times but is in poor agreement for early-to-mid times less than 20. It shows 6% error at t+ = 20 which increases to 40.3% at t+ = .001. All approximate solutions presented by Keltner and Beck lie within 2% of the present solution over their range of validity except the early-time T-based solution which shows an error of about 4.2% at time zero. This solution can be modified by letting the factor of 1.90484 (given in equation (34) of [18]) be replaced by nl/Z. The modi- fied solution provides better results for the early times (t+ < .005) and approaches the exact solution as t+ + 0. See Fig. 4.9. In Fig. 4.10 normalized values of averaged interface heat flux, 6:, is plotted versus dimensionless time. Normalization is obtained by dividing the heat flux values by k2°a'(Ti2'Ti1)' A study of the above results shows that the multinode surface element solution provides an accurate representation of the idealized intrinsic thermocouple. It is superior to other availble solutions in terms of accuracy and ability to treat the complete time domain. The method is most suitable for calculating the interface temperature and heat flux, particularly at early times when the two—dimensionality of the problem is significant. Figure Relative Error, Percent -l. 4.9 96 ”T l l T l ’— Early time [39] .- T-based[13] $333254 1:12.811. FD[4l] ‘\ ”" \ ,’ Early\\ .. time / .. Late time q-based \ I _ [18] I q based[18] I _ FD[ 27] _ I J 1 l I J 10‘3 10'2 10" 10° 10' 102 103 104 a t 4+ a Percent of relative error (with respect to SEM solution) in interface temperature as function of time for intrinsic thermocouple problem. 97 1.2 .a/k2(TiZ-Ti]) lo- 11 q Fisgure 4.10 Normalized area averaged interface heat flux histories for intrinsic thermocouple problem. CHAPTER 5 SEMI-INFINITE BODY WITH MIXED BOUNDARY CONDITIONS 5.1 Introduction In this chapter the transient thermal response of a semi-infinite body with the mixed boundary conditions of a step change of the surface temperature over an infinite strip and insulated elsewhere is considered. See Fig. 1.3a. This problem is similar to the problem of two semi-infinite bodies with identical pr0perties having different initial temperatures brought suddenly into the perfect thermal contact over an infinite strip and insulated over the rest of the contacting planes. The surface element method is employed to obtain the transient sollrtions for the interface heat fluxes and the thermal constriction resistance of the contact area. The solution has applications in the problems involving electronic cxxfling, strip welding, fins, and thermal contact conductance. To the autfiuar's knowledge there is no solution available for the above problem in the open literatures. Sadhal [45] has examined the related problem of two semi-infinite bodies having perfect contact over a series of equally-spaced infinite strips. The regions between these strips were insulated. By considering the planes of symmetry between the strips (see F“ig. 5.1) he solved the problem for large times by utilizing the Laplace: tranform technique. However, his solution is not valid for the sithations in which there is a small fraction of the interface in contact: (C 5 2), and consequently cannot be applied to the problem of two senri-infinite bodies with a single strip contact. 98 99 lines of symmetry 1 I insulated region L 412 V 11 1 i NJ Figure 5.] Geometry for Sadhal's problem. 100 5.2 Statement of the Problem The geometry being considered is shown in Fig. 1.38. A semi-infinite body is subjected to a uniform step temperature change, Tc’ over an in- finite strip with width 2a, on its surface. _The rest of the surface is insulated. The body is initially at the uniform temperature of T0, and the thermal properties are assumed to be independent of temperature. The describing equations are: T T g;3,+ géi,. %,%%. (5.2.1) T = T0 for t = 0; |x| 3_0; Z 3_0 (5'2°2) T = TC for t > 0; -a §_x 5_a; Z = 0 (5-2-33) %;-= 0 for t > 0; [X] > a; Z = 0 (5°2°3°) T = To for t > 0 as x +1:2 and z + 2 (5-2-4) inhere 'T denotes the temperature distribution, 0 represents the thermal diffusivity, x and z are space coordinates, and t is time. 5.3 Solution 3)! considering the problem symmetry about the x axis, the solu- tion is; required only for the region x 3_0, and consequently only half 0f the (:ontact area need to be discretized. To apply the surface ele- ment method, the surface region between x = 0 to x = a is divided into ten eleunents (each being an infinite strip) over each of which the heat flux iss uniform and at the center of which the prescribed temperature 15 TC. See Fig. 5.2. Equations (3.3.2a) and (3.3.2b) given in Chapter lOl / 311:3 I I Plane of symmetry [2 T0 Semi-infinite body Figure 5.2 Possible distribution of surface elements for the problem of semi-infinite body with constant surface temperature, Tc’ over an infinite strip. 102 3 can be used to determine the elemental heat fluxes at various values of time. It should be noted that for this case where only a single body exists, the elements of the initial temperature vector, T}, are the same and equal to TC - To. The influence functions for the geometry shown in Fig. 5.2 can be evaluated from the exact solution for the problem of a semi-infinite body heated with a constant heat flux over an infinite strip by apply- ing simple superposition. The procedure is the same as that presented in Section 3.3.2. It should be noted, however, that for the geometry of Fig. 5.2, the function ekji represents the temperature rise at element k due to a unit heat flux over the elements 1 through j (x = 0 tt) x = aj) at time ti’ while the influence function °kji denotes the temperature rise at element k due to a unit heat flux at only element j at the same time. The solutions for the local and the spatial average surface tem- perathre histories for a semi-infinite body exposed to‘a constant heat flux (aver an infinite strip are provided in Chapter 6. The solutions are griven by (6.4.3) and (6.4.4a), respectively. In the SEM solution Presefrted here, both local and average values of the influence function are used and the results are compared. A‘Leraged Interface Heat Flux 'Thea spatial average of the interface heat flux can be obtained by 5”"“109 the products of the elemental heat flux and the fraction of the total ccnqtact area occupied by the element. N E (tM) = 32-2-1 qu (Aj/AC) (5.3.1) 103 where qu is the heat flux at element j at time tM = M-At, Aj is the area of element j, Ac is the total contact area, and N is the number of elements. For the case where all elements are semi-infinite strips, (5.3.1) can be written as __ N q (tM) ‘ jgl qu ' [(aj ' aj-1)/a] (503-2) where aj, aj_1, and a are shown in Fig. 5.2. If further the elements are equally-spaced, the average heat flux can be given by, _ 1 N q (tM) = N321 qu (5.3.3) Thermal Constriction Resistance Based on the definition given in Section 3.5, the thermal constric- ticn1 resistance (per unit length of the strip) for the above problem can be written as, T - T RC (tM) = -9—:—°— (5.3.4) 2a q (tM) which 'is related to dimensionless heat flux by, _ 1 R (t ) - (5.3.5) C M 'f: 2k q (tM) where '7: _ . a q ’ k(%€:fi67 (5.3.6) 104 Equation (5.3.5) can be written in dimensionless form as, = __________ R 5 R .k (5.3.7a,b) Analytical Solution for the Average Interface Heat Flux By considering only one element along the interface an approxi- mate analytical solution for the average interface heat flux can be obtained by utilizing the Laplace transform technique. This method was first used by Keltner and Beck [l8] for the intrinsic thermocouple problem. Their solution was discussed in Chapter 4. Starting with Duhamel's integral for the average interface temperature one can write, T- = 03?? ft" (0, t- x) d). (5.3.8) 0 where the influence function $(O,t) is the spatial average temperature rise over the element (z=0) for a unit heat flux. Taking the Laplace transform of (5.3.8) gives [T -T 0] = $061M (o, t x) ax] (5.3.9) or [T -To] = 5 3(5) .2315) . (5.3.10) where the functions 3(5) and $(s) are the transforms of q(t) and $(t) respectively. Solving for'3(s) provides 105 31s) = ___£;::___ (5.3.11) 52 °'$(S) This equation can be written in dimensionless form as ...: = 1 + 2 ... (5.3.12) q (5) (5+)‘ - ¢+(s+) where +=0tt 7.51; ”1.25: t 51-, o a , s a (5.3.13a,b,c) and q+ is given by (5.3.6). Early Time Solution For small dimensionless times, the average influence function ¢+ is given by (6.4.4c) as + '1/t+ 1 ¢+|= 2 (t+/n)1/2 _ %__+ %—(t+)3e (1-§-t+) (5.3.14) For an error less than .033% the exponential term can be dropped for t+ < .3, and thus + <1" 2 2 (t‘L/w)“2 - ,E— (5.3.15) Taking the Laplace transform of (5.3.15) yields -:: + -3/2 1 = (S ) _ (5.3.16) substituting (5.3.16) into (5.3.12) and taking the inverse transform 106 gives -- t /TT2 1/2 q (1*) = (n 1+)‘1/2 + %-e erfc (-t+ /n) (5.3.17) This equation can be further simplified to q (t+) = (5 t+)'1/2 + 1- (5.3.18) 11 which provides less than 3% error for 1* < .01. Late Time Solution For late dimensionless times, the function ¢+ is given by (6.4.4b). For an error less .2% the second term on the right hand side of this expression ( ) can be dropped for t+ > 100, and one can write + 3nt 5+ [1n 1* + 3 - y], (y = 5772...) (5.3.19) =l|o—' Taking the Laplace transform of (5.3.18) gives ¢+ = _l:.(3 - 2y - in 5*) (5.3.20) 115 substituting (5.3.19) into (5.3.12) and taking the inverse transform by utilizing the approximate method (which is accurate for f(s) func- tions that vary slowly with in 5) given in [46] yields '77 + q (t ) = n(£n 2t+ + 3 - 2))”1 (5.3.21) which shows that the average surface heat flux increases as the inverse of the logarithm of time for large times. 107 5.4 Solution to the Interior Region Once the elemental heat flux histories have been determined, the solution to the interior temperature history, T(x,z,t), can be obtained by superimposing the total effect of all these heat fluxes. That is N M (x,z) T(x,z,tM) = TO + .2 .2 q.. A ¢j,M-i (5.4.1) where (x,z) (x,z) (x.Z) A ¢J.M-i = ¢j,M+1-i ' ¢j,M-i (5.4.2) (X.Z) and the influence function, ¢jM , is the temperature rise at point (x,z) due to a unit heat flux at element j at time tM = M-At. Equation (5.4.2) can be written in a more convenient form as, M-l TM(x,z) = M T0 - kgl Tk(x,z) (5.4.3) N M (x,z) + Z Z qji ¢j,M+1-i i=1 i=1 (X.2) The solution for the influence function, ¢jM , can be obtained from (6.3.18) given in Chapter 6. 5.5 Results and Discussion Two cases of equally-spaced elements and variable-spaced elements are examined. The first case is utilized to learn how the accuracy varies with the size of element, while the second case is used to obtain accurate results in the corner region near x = 1a. For both cases, elemental surface heat fluxes are determined for various values 108 of dimensionless time, ranging between t+ = .01 to 1000. At each time, the elemental heat fluxes are evaluated in twenty time steps. This means that for larger times, larger time steps are considered. For instance, to determine, qj(t+)'s, at dimensionless times of t+ = .01, 1, and 1000, the time steps of At+ = .0005, .05, and so are used, respectively. No. of Time Steps (NTS) = ——4——-= ———-= -——- = 20 (5.5.1) This substantially reduces the computational effort and provides more uniform accuracies for results (with respect to NTS) than the case where a small constant time step is used for the entire time range. To show how accuracies of results change with NTS, the case of equally-spaced elements is also examined with values of NTS being equal to 10, 5, and 2. In Fig. 5.3, the normalized area averaged interface heat flux is plotted versus 1/NTS, for different values of dimension- less times. Study of this figure leads to the following observations: 1. The accuracy of the heat flux histories vary linearly with fi%§3 and become more accurate as NTS increases. (The most accurate values can be obtained at fi%§-= 0 by employing linear interpolation). For NTS = 2, however, slight deviation (from a straight line) can be seen at early times, t+ < .1. The deviation is less than 2% at t+ = .01, and approaches zero as t+ becomes larger. 2. The slope of the straight lines shown in the figure are large at early times and become smaller for late times. This implies that, in order to have good accuracy at early times, NTS should be large, while the same accuracy can be obtained with smaller values of NTS, at later times. For instance, the required NTS for less than 3% errors in heat 109 .02 °' .05 l0 ' ' ' ' loo 0 1 o .25 5 1-55: NTS t+ Figure 5.3 Variation of the normalized area averaged interface heat flux with number of time steps used in the calculations. 110 flux histories are; 20, 10, and 5 for dimensionless times of t+ = .01, 1, and 10, respectively. To show how accuracy changes with the size of elements, the case of equally-spaced elements is considered with different number of ele- ments along the interface (N = 1, 2, 5, 10). The number of time steps (NTS) used was fixed and equal to 20. The results are shown in Figures 5.4 and 5.5. Fig. 5.4 shows the normalized interface heat flux, q+, versus Ax+ = Ax/a for various values of t+, while in Fig. 5.5, q: is plotted versus t+ for different cases of N = 1, 2, 5 and 10. A similar behavior to that of Fig. 5.3 is observed in Fig. 5.4. It can be seen that q: varies linearly with Ax+, and becomes more accurate as Ax+ approaches zero, especially for early times. The slopes of the lines decrease as t+ goes to infinity. This indicates that even few elements along the interface can produce good accuracy for large times. (This can also be observed in Fig. 5.5.) In order to show precisely the heat flux distribution across the strip, especially in the corner region near x = ta, ten variable-spaced elements were used with NTS = 20. The elements near the corner were smaller (about 1/4) than those close to the center line. Fig. 5.6 shows the normalized spatial variation of the heat flux across the strip at several times. Normalization is obtained by dividing elemental values by the value of the center line element which covers the region 0 5_x+ 5_.2. It can be seen that the region of uniform heat flux shrinks as t+ increases. At a dimensionless time about 30, the norma- lized heat flux distribution remains constant, which indicates for large times that the heat flux across the strip can be approximated by a product of a function of t+ and a function of x+. This behavior was 111 I I I I 1 I I I I 6. .02 4:— o .5- ! hf) .2 n *0— .05 l+ + 0' ‘1 w o] 2. r *4 "—H-o— . .2 A —-—o—E A TP—‘—f 4' JToo. J 1 1 1 4 . 1~ 1 1 1 0 o .5 l Ax+=Ax/a Figure 5.4 Variation of the normalized area averaged interface heat flux with the size of elements. 112 O 000' 00. 2222 u u u u -‘NU'I—" 0 1 1 1 24 10‘2 To”1 100 To1 102 103 t+- 9i"- a2 Figure 5.5 Results of average interface heat flux using different number of elements along the interface(local o's used). 113 6" g T T f“ 4.1- - OJ 0- \ o- 3.)- 2w- 1. 0 J 0 .5 1 + x =x/a Figure 5.6 Normalized heat flux distribution across the strip at various values of time. 114 also observed in the problems considered in Chapters 3 and 4. As mentioned earlier,both local and average values of the influence functions are examined in this problem. Table 5.1 shows the results obtained from the local-¢ solution and those found from the average-o solution. Both solutions use ten equally-spaced element with NTS = 20. The first column is the dimensionless time ranging between .01 to 100. The next two columns show the normalized values of the area averaged interface heat flux resulting from the local-¢ and the average-¢ solu- tions, respectively. The last column provides the corrected values of the average heat flux for the case where N+w (corresponding to the values of q: at Ax+ = 0, in Fig. 5.4). Comparison of the results from the second and third columns with the corrected values in the last column indicates that the average-¢ solution gives more accurate results than the local-o solution. The errors associated with the average-o solution are about 2/3 of those related to the local-o solution, which implies that the former solution is more appropriate than the latter one, par- ticularly when the number of elements is small. In Fig. 5.7 the results obtained from the multinode surface element solution (with 10 equally-spaced elements and NTS = 20) are compared with those evaluated from the early time analytical solution and the late time analytical solution given by (5.3.17) and (5.3.21), respec- tively. As it can be seen the multinode SEM solution is in very good agreement with the early time solution up to dimensionless time about 1., and matches closely the late time solution for times greater than 5. Notice that the early and late times analytical solutions match sur- prisingly well near a dimensionless time of 1. 115 TABLE 5.1 Comparison between the results obtained from the local-¢ solution and the average-o solution. q+ = ETijiigT Corrected Solution + local-¢ solution Av.-¢ solution .N = (w) t .01 6.1886 6.2228 6.2710 .05 3.0424 3.0578 3.0856 .1 2.2750 2.2858 2.3058 .5 1.3028 1.3076 1.3172 1. 1.0582 1.0616 1.0684 10 .6340 .6354 .6380 100: .4415 .4428 .4436 116 7 l I I I 6. _. 5. " O 4 ‘SEM Solution _ rot- I. 'U I... .¥ +1! 1 0' 3." " 2." " O 0 . E 1 _ Late-time ar y Solution 1 __ time . ° Solution 0 l 1 l 10'2 10" To0 To1 102' 103 + at t = a2 Figure 5.7 Comparison of the multinode surface element solution and the approximate analytical single node solutions for the early and the late times. CHAPTER 6 TEMPERATURES IN SEMI-INFINITE BODY HEATED BY CONSTANT HEAT FLUX qO OVER HALF SURFACE 6.1 Introduction This chapter presents the analytical solution for the transient temperature distribution in a semi-infinite body that is heated by a constant flux over half the surface and is insulated over the other half. The solution is a basic one in heat conduction and has a number of direct applications. It is important for the problems involving electronic cooling, contact conductance, fins and oil and gas wells. However, the prime motivating force behind the development of this solution is its application as a basic kernel for the surface element technique. The solution for the surface temperature is given in Carslaw and Jaeger [22]. See also [47]. Grimado, [48], considered the temperature distribution in a composite semi-infinite body subject to constant heat fluxes at the surface. His solution is expressed in terms of integrals of the form given by (6.3.11) which were solved numerically. Interior temperatures are more difficult to obtain than for the sur- face and are not given in a convenient closed form in the literature. The general solution presented here is valid for all times and any location in the entire region. It is given in terms of an integral which is shown to be conveniently and efficiently evaluated by two series expressions. 117 118 6.2 Mathematical Description The describing differential equation, initial and boundary condi- tions for an isotropic, homogeneous, semi-infinite body that is heated by a constant heat flux over half its surface are (see Fig. 6.1) 2 2 8 T a T laT ___+_—= —— (6.2.1) 8x2 822 a at T(x,z,O) = Ti (6.2.2) T -k g—Z-l z-o = q0 x < o (6.2.3a) = o x > o (6.2.3b) T(x,z,t) = T. t > 0 (6.2.4) 1 as 2 and x + m There is no heat generation inside the body and the thermal properties are independent of temperature. 6.3 Derivation of the Solution The exact transient temperature solution for the problem described by equations (6.2.1) through (6.2.4), can be obtained starting with the equation (1) given on page 258 of Carslaw and Jaeger [22], which gives the temperature rise for an instantaneous line source parallel to the y axis and passing through the point (x', z') for an infinite region. By considering the semi-infinite geometry of (z > 0) and integrating the effect of line sources at z=2' for the range of negative x (-00 < x < 0), the solution for the above problem is 119 uniform heat flux qo over l/2 space, x<0 v v U1VTWVIV AMAUDQAMAU/ 2=x Figure 6.l Geometry of a semi-infinite region heated by a uniform heat flux, qo, over half space m < x < O and z=0. 120 i: _9_ :0 e-U2 + (MW/480.13.! (5.3.1) 2nkt -w Equation (6.3.1) gives the temperature rise for an instantaneous plane source of strength 0 (J/mz) at time t=0; the source is parallel to the plane z=0, passes through the point z' and is over negative x's, that is, -w < x < 0. (The Caret on T denotes the instantaneous pulse case). Integrating (6.3.1) yields _ Q, -(z-z')2/4at f x ) 6.3.2 2pcp(01t)1/2 8 er C(2(at)1; ( ) —I > I If the source is at the plane 2' = 0, then (6.3.2) reduces to (6.3.3) Continuous heat flux For a time variable heat flux q(t) at the surface, the temperature at any time t can be obtained using Duhamel's integral, t T(x,z,t)-Ti = 5 q(1)$ (x,z,t-1) 61 (5.3.4) where 8 is T for a unit 0 at the surface. Then from (6.3.4) for the constant heat flux, qo, over the negative x surface, the solution for the temperature rise is 121 40(a/r)1/2 t _ dA T(X,Z,t)"Ti - 2k 6 W (6.3.5) X )dx 2[a(t—A)]1/2 2 e.Z /4a(t-A)erfc( This expression is valid for all x values less than, equal to, and greater than zero, and all 2 values equal to and greater than zero. Special cases a) For the special case of z=0, the solution is q 1/2 T(X,0,‘12)--T.i = T? (91%) [afflm 2 (6.3.6) - ———17—2—" E (l—H which is also given for zero initial temperature in Carslaw and Jaeger [22], equation (1), page 246. b) For the special case of x=0, temperature at the center line, the solution from (6.3.5) is ”2 ierfc (L) (6.3.7) Zvat i which is exactly one half of the solution for the body heated over the entire z=0 surface. General case The general case of z > 0 is not available in the literature. Introducing u as 122 z u = (6.3.8) 2[e(t-1)]1/2 into (6.3.5), the general solution becomes q0z w du -u2 xu T(x,z,t)-T. = f -—— e erfc (——) (6.3.9) ‘ 2kn112 z/2(6t)1/2 u2 2 To better describe the results, define the dimensionless groups as follows x = x , z = -—Ze———- (6.3.10a,b) 2(et)172 2(61)“2 T-T. 2 Z 1 p = — = "' a T = (6.3.10C,d) x X (co/0mm)”2 Notice that p is independent of time. Using these definitions (6.3.9) can be written in dimensionless form as 2 -u T(p.Z) = z fm-ly e erfc (g) du (6.3.11) 2 u The number of independent variables in (6.3.9) has been reduced from three, (x,z,t), to two dimensionless variables, (Z,p), in (6.3.11). It can be shown that the integration of (6.3.11) by parts yields 2 erf(X) 2 2 +2pZ [me-p y erf(y)dy X T(X,Z) = (11)“2 ieric(Z)-e'Z (6.3.12) x "““T72’ 21 (N) E1(X2+Z (For details see Appendix C.) 123 The functions erf('), erfc(-), ierfc(°), and E1(-) which appear in the equations (6.3.2) through (6.3.12) are the error function, the comple- mentary error function, the integrated error function, and the exponen- tial function, respectively defined by, f( ) 2 In _y2 d (6 3 13) er = e y . . Til/2 0 2 w - 2 erfc(n) = 1- erf(n) = f e y dy (6.3.14) 1/2 TT 11 . 2 ” Terfc(n) = f erfc(Y) dy (6.3.15) 517?.n -yd ”I -- (6.3.16) n y These functions are tabulated in [34], and are also available in com- puter libraries. The integral in the last term in (6.3.12) can be represented by a function H defined as H(X,p) = 1—1175- {” e"p y erf(y) dy = H(X,Z/X) (6.3.17) 11 (See [49] for related integrals). Then the general temperature solu— tion for a constant heat flux and 230 can be written as 2 T(X,Z) = (5)”2 ierie(z)-e‘Z erf(X) - (6.3.18) X 2 2 1/2 E (X +2 )+ (n ) ZH(X,Z/X) (10m 1 124 This expression is skew-symmetric with respect to x-axis and it can be shown that, TL(X,Z) = 2 n1/2 ierfc(Z) - 1+(x,2) (6.3.19) where the plus and minus subscripts indicate that the solution is for positive and negative values of X, respectively. The first term on the right hand side of (6.3.19) is the solution to the same problem if the entire surface was heated by a constant heat flux. The function, H(X,p), can be represented in series form for the three different regions indicated in Fig. 6.2. a) For the Region]p[_> 1 _2 .. (-1)"en(pzx2) _. ’) n-O an+1(2n+1)ep X or m n 2 2 = %. (5:)1r(n+1.p X ) (6.3.20b) n=0 p (2n+1) n where the functions en(~) and f(n,-) are the truncated exponential function and the incomplete gamma function, respectively, and are de- fined in [34] by 0 m' = n! e (6.3.21) (6.3.22) 125 Fig. 6.2 Geometry showing various regions of |p|l. 126 b) For the Region [pl < 1 Equation (6.3.20a) cannot be used for the case for which |p| < 1 since p appearing in the denominator of (6.3.20a) causes the summation to diverge. For this case of |p| < 1 the following expression is provided, n 2n+1 2 _ 2 °° ('1) p en(X) H(X,p)-1-erf(X)erf(pX)- 5' Z 2 (6.3.23a) "‘0 (2n+1) ex or w 2n+1 2 = - _ Z. ('1)':p F(n+lsx ) H(X,p) 1 erf(X)erf(pX) fl n20 (2n+1)n! (6.3.23b) c) For the lines [pl = 1 It can be shown for p=1 that H(X,p) is given by H(X,1) = -H(X,-1) = 1/2[1-erf2(X)] . (6.3.24a) which results in erfc z -22 1/2 1+(x=z,2) = [e +5 ierfc(Z)] Z 2 - ;T72'E1(ZZ ) (6.3.24b) The general solution given by equation (6.3.18), is valid for all times and any location in the entire region. Even though this expres- sion is valid for both positive and negative values of X, it is recommended only for X>0. For X<0 the complementary expression given by (6.3.19) can be employed. In the next section the behavior and 127 evaluation of the function H(X,p), is discussed. 6.3.1 Evaluation of the Function H(X,p) The evaluation of the function H(X,p) can be accomplished using (6.3.20a) and (6.3.23a) for the regions |p| >11 and |p| < 1, respectively. Efficient subroutines can be provided by utilizing several recur- sive relations. The summations in these equations are well-behaved for all values of X and p (within the appropriate p range). Table 6.1 provides values for H(X,p) for various X and p values. In Fig. 6.3, the function H(X,p) is plotted versus X for various values of p. There are several interesting points regarding the behavior of H(X,p) which can be seen from (6.3.17), Table 6.1 and Fig. 6.3; some of these are: a) Regardless of values of p, H(X,p) + 0 as X +.m. b) For large values of p compared to unity, H(X,p) goes to zero even for small values of X. c) For finite X, H(X,p) + 1, as p + 0. d) From (6.3.17) one can write, 1) H(X.-p) -H(X.p) (6.3.25a) 2) H(-X,p) H(X,p) (6.3.25b) e) The series expressions given by (6.3.20), and (6.3.23) for H(X,p) are valid only for positive values of p(X>0). For negative values of p(X<0), (6.3.25a) can be used. f) For X equal to zero, one can show that MNNNN-‘d-‘d wNN—l—l—l—l 128 TABLE 6.1 Values of function H(X,p) X p=0.l p=0.4 p=0.8 p=1 p=l.5 p=2 p=4 .00l .936549 .757762 .570446 .499999 .374333 .295l66 .155956 .0l0 .936543 .757737 .570396 .499936 .374239 .295040 .l55704 .050 .936390 .757l26 .569l75 .4984ll .37l954 .292001 .l49720 .100 .9359l3 .755222 .565378 .493676 .364907 .282707 .l32464 .200 .934020 .747676 .550466 .475202 .338043 .248406 .08ll70 .400 .926634 .7l857l .494960 .408240 .249144 .l48l82 .0ll453 .600 .9l495l .6736l6 .4l55l0 .3l7679 .l50202 .062803 .000439 .000 .88l962 .554237 .241324 .l44928 .030329 .004109 .000000 .200 .862366 .488908 .l682l8 .085664 .0l025l .000640 .000000 .400 .84l662 .424669 .lll029 .046577 .002883 .000072 .000000 .600 .820358 .363883 .069597 .023372 .000677 .000006 .000000 .000 .777l92 .257694 .023607 .004667 .000022 .000000 .000000 .200 .755665 .2l3246 .0l2800 .00l86l .000003 .000000 .000000 .400 .734287 .l74557 .006620 .000688 .000000 .000000 .000000 .600 .7l3096 .l4l345 .003265 .000236 .000000 .000000 .000000 .000 .67l373 .089686 .000689 .000022 .000000 .000000 .000000 TABLE 6.2 Values of the functions H(X,p) and ERFC(pX) for various X values p=0.l p=0.4 p= .8 X H(X,p) ERFC(pX) H(X,p) ERFC(pX) H(X,p) ERFC(pX) .00 .88l962 .887537 .554237 .57l608 .24l324 .257899 .20 .862366 .865242 .488908 .497250 .l682l8 .174576 .40 .84l662 .843053 .424669 .428384 .lll029 .ll32l2 .60 .820358 .820988 .363883 .3654l4 .069597 .070266 .00 .777l92 .777297 .257694 .257899 .023607 .023652 .60 .7l3096 .7l3l00 .l4l345 .l4l350 .003265 .003266 .00 .67l373 .67l373 .089686 .089686 .000689 .000689 TABLE 6.3 Number of terms in series of function H(X,p) given by equations (6.3.20a) and (6.3.23a) to obtain 8 decimal places accuracy (UN-”CO >< . O O p=0.6 p=0.8 p=0.98 p=l.02 p=l.2 p=l.5 .01 l l l l l l .l 2 3 3 3 3 3 7 8 10 10 l0 9 10 l4 l8 l8 l7 l5 13 20 27 28 25 l9 129 H(XJD Figure 6.3 Function H(X,p) versus X for different values of p. 130 g-arctan (l) for |p| > 1 (6.3.26a) H(0.p) = 1 - g-arctan (p) for |p| < 1 (6.3.26b) 9) For large values of X, say X > 3, function H(X,p) can be approximated by the complementary error function of pX, 2 2 H(X,p) 2 1—1172 I: e"p y dy = erfc (pX) (6.3.27) 71 In the next section this result is utilized to obtain a simplified ex— pression for the solution of T(X,Z). In Table 6.2, the values of H(X,p) and erfc(pX) (with 6 significant digit accuracy) are compared for var- ious X and p values. The number of terms required to obtain 8 decimal places accuracy, in the series expression for H(X,p) is shown in Table 6.3. It can be seen that the number of terms increases as X increases, and also as lp-ll decreases (p close to one). The number of the terms increases rapidly as X becomes greater than 10. But fortunately, the evaluation of the series is needed only for X < 3. (For values of X greater than 3, the H function can be approximated by (6.3.27)). Hence, no more than 30 terms are required to obtain 8 decimal places accuracy for evaluation of the series in the H function for the large range of any X and Ip-ll > .02. h) For small values of pzxz, H can be approximated by H(X,p) = g-[arctan(%) - sz] for p > 1 (6.3.28a) 131 and H(X,p) 2 1-erf(pX)erf(X) (6.3.28b) - g-[arctan(p)-pX2] for p < 1 6.3.2 Simplified Relations for |X[> 3 For |X| > 3, the general solution for the dimensionless tempera- ture given by (6.3.18) can be reduced to a simpler form. By introduc- ing (6.3.27) into (6.3.18) and noticing that pX = 2 one can write 22 x T+(X,Z) = e' erfc(X)-E;;T7§ E 1(x2+22) for x > 3 (6.3.29a) Equation (6.3.29a) is very accurate (to 6 decimal places), for X > 3. It also gives good accuracy down to X > 2 (to 3 decimal places). This accuracy (for 2 < X < 3) can be further improved to 5 decimal places, by subtracting an additional term, E-E1(X2+Z2), from the right hand side of the (6.3.27), which results in (6.3.29a) becoming: 2 2) erfc(X) - Z X 1+ 2 2 (TI)1 ) 1+(x,2) = e' E x +z (6.3.29b) 1( Furthermore, due to the nature of the terms exp(-Z2) erfc(X) and X(fl)-1/2 E1(X2+22) in equation (6.3.29a), there are negligible contri- butions of these terms for X > 4. It can be shown that for X > 4 the 8 and the second term is less than first term is less than 1.6 x 10' 1.5 x 10-8. Hence, for say 8 decimal places accuracy, these terms can be dropped and the equation becomes independent of X. Then one can write 132 T+(X,Z) 0.0 for X > 4 (6.3.30a) T_(X,Z) 2ierfc(Z) for X < -4 (6.3.30b) These equations also give excellent accuracy down to |X| > 3. 6.3.3 Evaluation of T(X,Z)and Discussion of the Results The evaluation of T(X,Z) given by equation (6.3.18) is not diff- cult if a computer is employed. The first three terms in the equation can be directly calculated using a computer library for the error and the exponential functions. The evaluation of the last term (H(X,p)) also can be easily obtained as discussed in the previous section. For the special cases of Z=0, X=0, p=1, X > 3, and |X| > 4 the simplified equations (6.3.6), (6.3.7), (6.3.24b), (6.3.29a), and (6.3.30) are provided, respectively. Table 6.4 provides values of the dimensionless temperature T(X,Z), for various values of the dimensionless X and Z. In Fig. 6.4, T(X,Z) is plotted versus X for various Z values. It can be seen that the solution is skew-symmetric with respect to X-axis. The solution is well-behaved for all times and any location over the entire region. Fig. 6.5 shows the lines of the dimensionless isotherms in the X and Z coordinates. It can be seen that for |X| > 1.5 the solu- tion is almost independent of X. 6.4 Other Cases of Boundary Conditions The exact solution given herein is important because it is a basic geometry in heat conduction. There are many other possible cases that can be obtained using this solution. It can be utilized as a building block for a number of other boundary conditions for both semi-finite awmada 133 TABLE 6.4 Values of dimensionless temperature for solution given by equation (6.3.l8), for different values of X and Z 4.000 3.000 2.000 1.000 -.600 -.200 .100 .050 .010 .000 .010 .050 .100 .200 .600 .000 .000 .000 .000 Olll DOOM-H dddddddddm Z=0.0 .000000 .999999 .999587 .966475 .866022 .525251 .340279 .209176 .059991 .000000 .940009 .790824 .659721 .474749 .133978 .033525 .000413 .000001 .000000 ddddd Z=0.2 .370489 .370488 .370095 .338908 .246380 .959476 .832317 .760328 .700368 .685245 .670122 .610162 .538172 .411013 .124109 .031582 .000394 .000001 .000000 Z=0.6 .552776 .552776 .552506 .532887 .481679 .360987 .319722 .298190 .280757 .276388 .272019 .254587 .233054 .191790 .071098 .019889 .000270 .000001 .000000 TABLE 6. 5 Z=1.0_ .178148 .178147 .178019 .169798 .151056 .113288 .101382 .095254 .090312 .089074 .087836 .082894 .076765 .064860 .027092 .008350 .000129 .000000 .000000 Z=1.2 .092341 .092341 .092264 .087610 . .077569 .058273 .052313 .049253 .046788 .046171 .045553 .043088 .040029 .034068 .014773 .004731 .000078 .000000 .000000 Z=2.0 .003467 .003467 .003463 .003249 .002849 .002153 .001946 .001840 .001755 .001734 .001712 .001627 .001521 .001314 .000618 .000218 .000004 .000000 .000000 Values of dimensionless temperature, T(x ,q+ ,t+ ), for various values .000 .010 .100 .200 .600 .800 .000 .200 .600 .000 .000 .000 of x+ and 2 due to a strip source, at time t ddddd—J z+=0.0 .932950 .932908 .928771 .915898 .743545 .524128 .999587 .474607 .133964 .033524 .000413 .000001 z+=0.2 1.307326 1.307287 1.303467 1.291610 1.137423 .958406 .684851 .410877 .124096 .031581 .000394 .000001 z+=0.6 .512998 .512977 .510862 .504375 .429491 .360259 .276118 .191696 .071088 .019888 .000270 .000001 z+=l.0 .161448 .161440 .160683 .158384 .133704 .112948 .088945 .064814 .027087 .008350 .000129 .000000 .25 z+=l.2 .082878 .082874 .082475 .081265 .068534 .058070 .046093 .034040 .014770 .004731 .000078 .000000 z+=2.0 .003031 .003031 .003016 .002970 .002508 .002142 .001729 .001312 .000617 .000218 .000004 .000000 134 j l j T l 2 ‘0 .. h .1 d .2 .. |— .3 d :I .2 1__ __ E? .4 -.5 - I- .6 —I '.8 _ l. #2. _ l.5 0 s; I \\__ 1 -3 -2 -l O l 2 3 Fig. 6.4 Dimensionless temperature T(X,Z) versus X for different values of Z in semi-infinite body heated by a uniform heat flux over half space X<0, Z=0. 135 r —+ N —-—1 J 2 Figure 6.5 Lines of dimensionless isotherms in semi-infinite body heated by a uniform heat flux over half space x < 0. 136 and finite geometries. Fig. 6.6 is given to illustrate some of these, for both semi-infinite and finite geometries. For most semi-infinite cases the solution can be obtained by simple superposition. For example, the solution for Fig. 6.6b, is the solution for a semi- infinite body that is uniformly heated minus the solution for the semi- infinite body heated across an infinite strip with the same heat flux. For the finite cases the solutions can be found by using the method of images. In this method, a number of sources and sinks are distributed inside the body such that the corresponding boundary conditions are satisfied (see [22] and [50] for more details). 6.4.1 Application to Strip Heat Source As suggested by the foregoing there are other boundary conditions that can be obtained using the above results. One case of particular interest is for the semi-infinite geometry heated by a constant heat flux over an infinite strip with the width 2a and insulated elsewhere. See Fig. 6.7. The solution to this case is needed in Chapter 5 in order to find the required influence functions. This solution can be found from the solution given by equation (6.3.18) by applying simple superposition to get: - 1 41+ +-1 1 1 T(X+,Z+, t+)=e Z / [-erf (21:111721 + erf (2(t11T721] + + x-1 (x+-1)2+ + x+1 E (x+1)atz 2(n t )1/2 E1[ 41+2 2] 2(fl t )1/2 1 41+ 2] 1T21/2+ + ( p)- H( " +1111 (6.4.1) 2(1)11H77[ 2(1”)171’pz1’(*)171 137 a) b) c) | [Ill/lllll‘ ‘1 ll \\\ ‘ \\ \ ’7////7///////// T=0 g) h) 1) Figure 6.6 Various possible cases that can be treated using the solution givenfor figure 6.1 as a building block. 138 uniform heat flux go ///////1111///" I V l I '//////j////// ; 2a); Figure 6.7 Geometry of a semi infinite body heated by a uniform heat flux, qo, over an infinite strip with width 2a. 139 where - x/a, z+ E z/a, (6.4.2a,b) X | I ‘ at/az, p E z/x C” II z+/x+ (6.4.2c,d) and T is defined by equation (6.3.10d). For the special case of 2+ = 0, the solution becomes: + + + + x +1 x -1 T(x ,0,t ) = erf - erf 2(1")172 2(1:“)172 (6.4.3) + x++1 E [(x++1)2]_ x+-1 E [{x+-l)2] 1 4 2(111”)1/2 41 2011+)”2 1 1 which is also given in [22], equation (3), page 264. Equation (6.4.1) is valid for all times and any location in the entire region. Table 6.5 provides values of the dimensionless tempera- ture T(x+,z+,t+), for various values of x+ and 2+ at time t+ = .25. Results are also shown in Fig. 6.8. It can be seen that for large values of 2+ (deep in the body), the temperature distributions be- come flatter as expected. In Fig. 6.9, some values of dimensionless surface temperature T(x+,0,t+) are plotted versus x+ for various values of t+. Another way of illustrating the results is provided in Fig. 6.10, which shows the isotherms in the x+ and 2+ coordinates for time t = 1. In the application of the surface element method the solution for the spatial average temperature over the heated region is sometimes utilized. At the heated surface this average temperature is 140 2 l I I I l z+=0 _ + t=a25 1 .l ...) .2 Q L _ +O N. — 4-" 1,— 3: .4 .6 ":8 1—1.0 \ - \, o 1.5 . j f§\\ I o 1 x+ 2 3 + + + + . Figure 6.8 Dimensionless temperatures T(x ,z ,t ) versus x , for various values of 2+ at t+=.25. l4] 1 —] T I I T t*=.0525 .25 ‘ .5 I— 1. - P _ +A *1 c .— +;' 1"' 10. E: 100. o 1 \j\\:::::::EEEIEEEEE::::33;Z:;;:-—-—-+_.______ o 1 2 3 o O O + + O D Fig.6,9 DimenSionless surface temperature T(x ,0,t ) versus dimenSIOnless + O O I I x , for various values of d1menionless times. 142 ._u+u Low .ou+~ um nwgum muwcwmcw cm Lm>o xzym “mm; Ecomwcs w x3 vmuow; Avon muwcwecw1wsmm cw macmcuomw mmwpcowmcoswo op.o mesmwu 143 El0,t+) = 2+(t+)1/2{21erfc(t+)'1/2 -fi'1/2[1+EZ(1/1+)]} (6.4.4a) For large times the average temperature is found from (6.4.4a) to be 4 1 , n1 - —+3-v], qu/k TI“ 3.1:"- Tlo,t+)-Ti 1 [1 + (6.4.4b) (y=.5772...) which shows that the temperatures near the heated surface increases as the logarithm of time for large times. For an error less than .2% the second term in (6.4.4b) can be dropped for t+ > 100. For small times (6.4.4a) can be approximated by TLO,t+)-T1 qoa/k + + = 2(t+/n)1/2-%—+%(t+)3 e"1/t (1%?) (6.4.4c) For an error less than .033%, the exponential term can be dropped for 1* < 0.3. CHAPTER 7 SUMMARY AND CONCLUSIONS A transient multinode surface element method for solution of two and three dimensional heat conduction problems with linear boundary conditions has been presented. The method uses Duhamel's integral and involves the inversion of a set of Volterra integral equations, one for each surface element. The use of Duhamel's integral requires that the describing differential equations be linear. In Chapter 1, Duhamel's theorem was introduced and integral equa- tions with temperature-based kernel and heat flux-based kernel were presented and discussed. Though both types of kernels can be employed in the SEM, only heat flux-based kernels were used in the examples given in this dissertation. The method is applicable to homogeneous and composite geometries with perfect or imperfect contact. In Chapter 2, the multinode sur- face element method formulations for two arbitrary geometries in per- fect (or imperfect) contact over part of their boundaries were developed and discussed. Through the use of piecewise uniform approximations of time and space variables, it was demonstrated how the integral equations presented in Chapter 1 can be transformed to a set of algebraic equations. To show the flexibility and applicability of the method to two- dimensional homogeneous and composite bodies, the multinode surface element formulations developed in Chapter 2, were applied to three dif- ferent problems given in Chapters 3, 4 and 5. In Chapter 3, the SEM was employed to solve the problem of two semi-infinite bodies (initially at two different temperatures) that are suddenly brought together over a small circular region and insulated 144 145 elsewhere. The intrinsic thermocouple problem was investigated in Chapter 4. Values for the thermal constriction resistance of the con- tact area (in Chapter 3), and the interface heat fluxes and tempera- tures (in Chapters 3 and 4) were presented for complete range of dimen- sionless time. The results obtained from the SEM solutions had excellent agreement with existing analytical and numerical solutions. In Chapter 5, the method was further applied to the problem of a semi-infinite body with mixed boundary conditions of a step change of surface temperature over an infinite strip and insulated elsewhere, for which analytical solutions do not currently exist. It was found that the accuracy of results varies linearly with the size of element and the size of the time step used in the numerical computations. For each of these problems the multinode surface element method performed very well. It was found that the heat flux and temperature at any point along the interface can be readily obtained utilizing this method. The method is most suitable for calculating interface tempera- tures and heat fluxes for the geometries connected over relatively small portion of their surface boundaries. The results obtained in Chapters 3, 4 and 5 showed that very high accuracy is attainable with a relatively small number of surface elements. This feature makes the surface element method superior to the alternative numerical methods such as FDM, FEM, or BIEM for the type of problems considered. In the SEM only the interface between the two geometries requires discretiza- tion as opposed to the discretization of the whole domain needed in the finite difference and finite element methods or discretization of the whole boundaries in the BIEM. This in turn reduces the size of the numerical calculations and computer time. The advantage of considering only interface nodes requires, 146 however, that the "building blocks" (or kernels, or the "influence functions"), ¢'s be known for the geometries under consideration. For many geometries the kernels are known or can be obtained by well-known analytical or numerical procedures. The influence functions required for the geometries of Chapters 3 and 4 were obtained from the known solutions given in [33,35], while those needed in Chapter 5 were found from the exact solution for the problem of a semi-infinite body heated by a constant heat flux over half the surface which was presented in Chapter 6. Three computer programs were written in the course of this research, RAN3, RAN4, and BBY9. (Documentation is available from the author.) Both uniform and nonuniform nodes can be employed in the programs. The programs have the capability of treating up to twenty interface nodes and can be applied to both symmetrical and nonsymmetrical geometries. They also can be adapted to be used for a number of different cases with different geometries, providing the proper influence functions, ¢'s, are used. Because of the convolution behavior of the integral equa- tions, the influence functions need to be calculated at each time step and stored for later use. This can result in a large computational effort, particularly if the number of the time steps (NTS) is large. However, it was found that (see Chapter 5) in most cases (except for very early times, t+ < .001) a value of NTS between twenty to forty can produce very accurate results. It was further shown that for large times (t+ > 10), even a few time steps (2 5) can provide good accuracy. Further improvement in terms of capability, accuracy and computer time is needed in the multinode surface element method. The following recommendations are given for future work: 1. Though the method has been applied to the problems with linear 147 boundary conditions it can be used for nonlinear boundary conditions. In particular, the radiative boundary condition can be investigated. The assumption of constant temperature and heat flux over each element and time interval as was done in this dissertation is not necessary and higher order interpolation functions (linear, quadratic, etc.) can be introduced to improve the accuracy. Work is needed to explore new methods for reducing the computer time. One possibility is to use the Laplace transform techni— que to avoid the convolution behavior of the numerical compu- tations. New influence functions for more basic geometries (such as semi-infinite body heated by a constant heat flux over a rectangular area or infinite cylinder heated over a portion of its circumferential area) need to be derived. The basic theory that is given in this dissertation is also applicable to three-dimensional problems and the problems with multiple interfaces. Some of these cases should be investigated. Finally, the method can be extended to include some specific problems of convection heat transfer and fluid flow problems. This in turn will make the method more competitive with already existing numerical procedures. APPENDICES APPENDIX A DERIVATION OF EQUATIONS (2.5.113) AND (2.5.11b) Equation (2.5.11a) can readily be obtained by considering fbr the first time step (M = 1), the vectors E and F (given by (2.5.8c) and (2.5.8d), respectively) are zero. E ='E = O for M = 1 (A1) Substituting (A1) into (2.5.8b) and then its result in (2.5.10) yields q1 = C T. (A2) which is the same as (2.5.11a). To show how (2.5.11b) is derived, equation (2.5.9) is expanded for different values of M. By introducing (2.5.8b,c,d) into (2.5.9), for M = 1,2,3,..., one can write atM= 1, Eq1=Ti (A31) atM=2, EEZ=T1+E|atM=2-F|atM=2 (A32) atM=3, EE3-Ti+-E—latM=3'FlatM=3 (A33) atM -1, '54}pr +E| at '44-?) at 11-1 (A3,M-1) atM , EEM=T1+E|atM-f|atM (A3M) By adding all these equations together, (A31) through (A3M), and notic- ing that 148 149 E I at M = F I at M-l + QI qM-I (A4) it can be shown that = M _ _ : M’l _ _ C [.E “1] ' MTi + 11 .5 Q1 F l'at M (A5) 1-1 1-1 or _ M'1_ ='1_ :‘1 : M-1_ =‘1 __ qM + .2 qi = M C Ti + C 61 Z q1 - C F I at M (A6) 1:1 1:1 -1 substituting for E,‘F, and E T} from (2.5.8a), (2.5.8d), and (A2), respectively yields _ _ = M-l =-1_ qM=Mql+B[.Zlq,l-c F (A7) 'I: Where the matrix B is defined as = :‘1 = B = H 4 (A8) Equation (A7) is the same as (2.5.1b) and is valid for M.: 2. Notice that the vector E'is a function of time and should be calculated at each time step. APPENDIX B EVALUATION OF I(r+,b+) GIVEN BY (3.4.5b) A new method is developed to evaluate the values of I(r+,b+) given + + + JO(Air /b ) J1 (AI/b ) NV 3b ) = 1 [A1- Jo (A312 (Bl) .i + + m E To evaluate this term effectively the summation is divided into two different parts: I(r+,b+) = I1(r+,b+) + 12(r+,b+) (82a) where I -1 + + + 1111+1b+1 ‘ mix J°(AI: /: )(:1)§:1/b ) (BZb) i=1 i o i and m J (A.r+/b+) J (A./b+) I(r+,b+) = 0 1 1 1 (82¢) 2 iglmax 111 Jo (11732 The first term in (82a), 11, is evaluated by direct series expan- sion. The second term (the correction term or the tail of the summa- tion), 12, can be simplified in the following manner. For large values of A, one can write [34] 150 151 2 1/2 n 2 1/2 1 J0(Ai) 2 (EA?) cos (Ki - a) - (EX?) ('1) (34) + 1/2 + + .3 2b + + 1 JO(Air /b ) — (NA r+) cos (Air /b - 4) (BS) 1 and + 1/2 + ~ .22. +-_3_11 J1(Ai/b ) ~ (“Ai) cos (xi/b 4 ) (36) + 1/2 _ Z_b__ 1 _ (“Ai) Sln (Xi/b - 4) Substituting (B3) through (B6) into (BZc) yields 00 + + + ~ b 12(1” ,b ) "' 2121 (Y‘+)1/2 max + + n . + n cos (Air /b - a) Sln (Xi/b - 1) A12 or + . l-r + m 51" ( + Ai) I (Y'+ b'I') : J— 2 b 2 1 (r+)1/2 1:1 112 max (37) 1+r+ cos ( b+ 1i) A12 From the composite midpoint rule in numerical analysis one can show that of 1(1):.117) mm +_rflix_ (BB) 152 Introducing (B8) into (B7) yields + + b sin AA 12(r ’b ) - “(r+)1 [IA )2 d A max (89) + + + _fm cos BA dx] + J0(Amaxr lb ) J1 ()‘max/b ) Amax A2 [Amax JoUmaxH2 where + + A = 1': , B = 131;— (BlO,a,b) b b The integrals on the right hand side of (B9) can be evaluated as w . sin (AA ) 1 £2754 dA = 1 "1‘1" - A Ci (Axmax) (Blla) Amax max w cos (BA ) cos BA _ max . Bn I> T d) - (max +8 51 (meax) - -2-— (3111)) max where Si(-) and Ci(-) are the sine and cosine integrals given by equa- tions (5.2.8) and (5.2.9) of [34], respectively. Finally, substituting (Blla) and (Bllb) into (89) gives, r+’b+) 1+r+ + b 51” (AAmax 2“(31/2 “(r+)1/2 [ x I ) - cos (BAmax):I l2 ( 2 max ;Z;;%T7§-[(1-r+)Ci(AAmax)-1(1+r) 51 (meax)] + + + Jo(xmaxr /b ) J1 ()‘max/b ) [Amax Jo (Amax)]2 (812) + 153 As it can be seen from (812), the correction term, 12, can readily be evaluated providing Amax’ r+, and b+ are known. To compare the new method presented here with the direct series expansion method, I(r+,b+) was evaluated for values of r+ = .1 and b+ = 10 with four deci- mal places accuracy. It was found that the number of terms needed in the direct series expansion approach is about 6-7 times larger than those used in the new alternative method. This in turn reduces the size of the numerical calculations and computer time. APPENDIX C EVALUATION OF THE INTEGRAL IN EQUATION (6.3.11) The integral in equation (6.3.11) can be written as I :12 UE-e' erfc (g) du * (c1) Integrating equation (C1) by parts, by letting dw = é3-du (C2a) -U2 v = e erfc (—) (CZb) yields I =-—-e ”2 erfc (—) (m — 2 fm.l.e‘U2(1 + %39du z P(n)1;z z ” w 2 u _ 1 -Z2 - 2 f e erfc (-) du — 7—e erfc (X) (C3) Z - - ———-Z—7—-J'm1—1--e'u2(1 + 530 du - 2 [m e"u erfc (E) du 90111 2 Z Z p where (6.3.10c) was used (p = Z/X). 1 By using the substitution of u2(1 + 63) = y and interchanging the order of integration, the first integral in (C3) can be wirtten as 2(1 + -—) = l—fm .EEX d = 1 E (X2+Zz) (C4) 2 y 2' 1 II = [Z e'u p2 22+x2 y It now remains to find an expression for the second integral on the 154 155 right hand side of (C3). By using the relation erfc(o) = 1 - erfc(-), the integral can be written as (D co _ 2 _ 2 f e u du - f e U erfc (% Z . III ) du (1:5) (111/2 ‘77- m -u2 u erfc(Z) - f e erf (E) du Z By using the substitution of %-= y and interchanging the order of inte- gration of the integral on the right hand side of (C5) one can write on 1/2 III = 1%1- erfc(Z) - p IX e_p2y2 erf(y) dy (C6) substitute (C5) and (C4) into (C3) yields )1/2 -Z2 I = £%—- ierfc(Z) - %- erf(X) --;Z%;T77 E1 (X2+72) ” 2 2 (C7) + 2p f e-p y erf(y) dy X Finally, using (C7) into (6.3.11) yields _ 2 T(X,Z) = (411/2 ierfc(Z) - e Z erf(X) m (C8) _ :X 2 2 'szz E (X +Z ) + 2pZ f e erf(y) dy (“)1/2 1 X which is the same as equation (6.3.12). LIST OF REFERENCES IO. 11. 12. 156 LIST OF REFERENCES Ozisik, M. N., Boundary Value Problems of Heat Conduction, Inter- nation Textbook Co., Scranton, Pennsylvania, 1968. Myers, E. G., Analytical Methods in Conduction Heat Transfer, McGraw-Hill, New York, 1971. Zienkiewicz, O. C., The Finite Element Method, Third Edition, McGraw-Hill Book Co., London, 1977. C. A. Brebbia, (Ed.), Recent Advances in Boundary Element Methods, Pentech Press, London, 1978. P. K. Bannerjee and R. Butterfield, (Eds.), Developments in Boundany Element Methods - 1, Applied Science Publishers, London, 1979. T. A. Cruse and F. J. Rizzo, (Eds.), Boundary-Integral Equation Method: Computational Applications in Applied Mechanics, American Society of Mechanical Engineers, New York, 1975. C. A. Brebbia, The Boundary Element Method for Engineers, Pentech Press, London, 1978. C. A. Brebbia and S. Walker, (Eds.), Boundary Element Techniques in Engineering, Butterworths & Co., London, 1980. Schneider, G. E., "Thermal Constriction Resistance Due to Arbitrary Contacts on a Half-Space-Numerical Solution of the Dirichlet Problem," Progress in Astronautics and Aeronautics: Thermophysics and Thermal Control, Vol. 65, edited by R. Viskanta, New York, 1979. PP. 103-119. Schneider, G. E. and LeDain, B. L., "The Boundary Integral Equation Method Applied to Steady Heat Conduction with Special Attention Given to the Corner Problem," AIAA paper No. 79-0176, presented at the 17th Aerospace Sciences Meeting, New Orleans, La., Jan. 15-17, 1979. Khader, M. S., "Heat Conduction with Temperature Dependent Thermal Conductivity," presented at the 19th National Heat Transfer Con- ference, ASME paper No. 80-HT-4, 1980. Khader, M. S. and Hanna, M. C., "An Iterative Boundary Integral Numerical Solution for General Steady Heat Conduction Problems," Journal of Heat Transfer, Vol. 103, pp. 26-31, Feb. 1981. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 157 Rizzo, F. J. and Shippy, D. J., "A Method of Solution for Certain Problems of Transient Heat Conduction," AIAA Journal, Vol. 8, No. 11, Nov. 1970, pp. 2004-2009. Shaw, R. P., "An Integral Equation Approach to Diffsuion," Inter- national Journal of Heat and Mass Transfer, Vol. 17, pp. 693-699, 1974. Chang, Y. P., Kang, C. S., and Chen, D. J., "The Use of Fundamen- tal Green's Functions for the Solution of Problems of Heat Conduc- tion in Anisotropic Media," Journal of Heat and Mass Transfer, Vol. 16, 1973, PP. 1905-1918. Wrobel, L. C. and Brebbia, C. A., "A Formulation of the Boundary Element Method for Axisymetric Transient Heat Conduction," Inter- national Journal-Of Heat and Mass Transfer, Vol. 24, No. 5, pp. 843-850, 1981. Yovanovich, M. M., and Martin, K. A., "Some Basic Three-Dimensional Influence Coefficients for the Surface Element Method," AIAA paper No. 80-1491, AIAA 15th Thermophysics Conf., July 14-16, 1980, Snowmass, Co. Keltner, N. R., Beck, J. V., "Unsteady Surface Element Method,":lL of Heat Transfer, Vol. 103, No. 4, pp. 759-764, Nov. 1981. Beck, J. V. and Keltner, N. R., "Transient Thermal Contact of Two Semi-Infinite Bodies over a Circular Area," Paper No. AIAA-81-1162, presented at AIAA 16th Thermophysics Conf., June 23-25, 1981, Palo Alto, CA, to be published in the Progress in Aeronautics and Astro- nautics Progress Series 1982. Kays, W. M. and Crawford, M. E., Convective Heat and Mass Transfer, Second Edition, McGraw-Hill Book Co., New York, 1980. Colladay, R. S., "Inverse Problems Related to Transient Convection in the Thermal Entrance Region between Parallel Plates," Ph.D. Dissertation, Department of Mechanical Engr., Mich. State Univer- sity, 1969. Carslaw, H. S. and Jaeger, J. C., Conduction of Heat in Solids, Cambridge University Press, Cambridge, 1959. Strong, A. B., Schneider, G. E. and Yovanovich, M. M., "Thermal Constriction Resistance of a Disk with Arbitrary Heat Flux - Finite Difference Solution in Oblate Spheroidal Coordinates," AIAA pro- gress in Astronautics and Aeronautics, "Heat Transfer with Thermal Control Application," MIT Press, Vol. 39, 1975, Editor, Dr. M. M. Yavanovich, University of Waterloo. Normington, E. G. and Blackwell, J. U., "Transient Heat Flow from Constant Temperature Spheroids and the Thin Circular Disk," Quarterly Journal of Mechanics and Applied Mathematics, Vol. 17, Part 1. pp. 65-72, 1964. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 158 Blackwell, J. H., "Transient Heat Flow from a Thin Circular Disk Small-Time Solution," Journal Australian Math. Soc., Vol. 14, pp. 433-442, 1972. Keltner, N. R., "Transient Heat Flow in Half-Space Due to an Iso- thermal Disk on the Surface," J. Heat Transfer, Vol. 95, Ser. 3, No. 3. PP. 412-414, 1973. Keltner, N. R., Heat Transfer in Intrinsic Thermocouples Applica- tion to Transient Measurement Errors, Report SC-RR-72-O719, Sandia National Laboratories, Albuquerque, NM, January 1973. Schneider, G. E., Strong, A. B. and Yovanovich, M. M., "Transient Heat Flow from a Thin Circular Disk,‘| AIAA 10th Thermophysics Conference, Denver, Colorado, May 27-29, 1975, AIAA paper No. 75-707. Marder, B. M. and Keltner, N. R., "Heat Flow from a Disk by Separa- tion of Variables," Numerical Heat Transfer, Vol. 4, pp. 485- 497, 1981. Heasly, J. H., "Transient Heat Flow between Contacting Solids," Int. J. Heat Mass Transfer, Vol. 8, 147-154, 1965. Schneider, G. E., Strong, A. B., and Yavanovich, M. M., "Transient Thermal Response of Two Bodies Communicating through a Small Cir- cular Contact Area," Int. J. Heat Mass Transfer, Vol. 20, pp. 301- 108, 1977. Sadhal, S. 5., "Transient Thermal Response of Two Solids in Contact over a Circular Disk," Int. J. Heat Mass Transfer, Vol. 23, pp. 731-733. 1980. Beck, J. V., "Large Time Solutions for Temperatures in a Semi- Infinite Body with a Disk Heat Source," Int. J. Heat Mass Transfer, Vol. 24, pp. 155-164, 1981. Abramowitz, M. and Stegun, I. A., Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, Natibnal Bureau off Standards, Applied Mathematics Series, Vol. 55, 1964. Beck, J. V., "Transient Temperatures in a Semi-Infinite Cylinder Heated by a Disk Heat Source," Int. J. Heat Mass Transfer, Vol. 24, No. 10. PP. 1631-1640, 1981. Burnett, D. R., "Transient Measurement Errors in Heated Slabs for Thermocouples Located at an Insulated Surface," Journal of Heat Transfer, Trans. ASME, Series C. Volume 83, p. 505, 1961. Larson, M. B., "Time-Temperature Characteristics of Thin-Skinned Models as Affected by Thermocouple Variables," University of North Dakota, 1967. Published as NASA-CR-91453. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 159 Larson, M. B. and E. Nelson, "Variables Affecting the Dynamic Response of Thermocouples Attached to Thin-Skinned Models," Journal of Heat Transfer, Trans. ASME, Series C., Volume 91, p. 166, 1969. Henning, C. D. and R. Parker, "Transient Response of an Intrinsic Thermocouple," Journal of Heat Transfer, Trans. ASME, Series C, Volume 39, p. 146, 1967. R. Parker, ”Rapid Phase Transformations in Titanium Induced by Pulse Heating," Transactions of the Metallurgical Society of AIMI, pp. 1545-1549, August, 1965. Shewen, E. C., "A Transient Numerical Analysis of Conduction be- tween Contacting Circular Cylinders and Halfspaces Applied to a Biosensor," MS Thesis, University of Waterloo, 1976. Geidt, W. H. and Nunn, R. H., comments given in reference [39]. Bickle, L. W., "A Time Domain Deconvolution Technique for the Cor- rection of Transient Measurements," Report SC-RR-710658, Sandia Laboratories. Bickle, L. W. and Keltner, N. R., "Techniques for Improving Effec- tive Response Times of Intrinsic Thermocouples," Report SC-RR-710146, Sandia Laboratories (Sept. 1971). Sadhal, S. S., "Unsteady Heat Flow Between Solids with Partially Contacting Interface," J. of Heat Transfer, Vol. 103, pp. 33-35, Feb. 1981. Schapery, R. A., "Approximate Methods of Transform Inversion for Viscoelastic Stress Analysis," Proceedings of the Fourth U.S. National Congress of Applied Mechanics, pp. 1075-1085, 1961. Jaeger, J. C., "Approximations in Transient Surface Heating," Australian Journal of Scientific Res., Ser. A, S, pp. 1-9, 1952. Grimado, P. B., "The Transient Temperature Field in a Composite Semi-Space Resulting from an Incident Heat Flux," Quarterly of Applied Mathematics, Vol. XXXI, No. 4, Jan. 1974. -x dx and Z Rosser, J. Barkley, Theory and Applications of I e 0 Z _ 2 2 y -x2 I e p y dy f e dx Mapleton House, Brooklyn, NY, 1948. O O Beck. J. V., "Average Transient Temperature within a Body Heated by a Disk Heat Source," Heat Transfer, Thermal Control and Heat Pipes, W. B. Olstand, Ed., Vol. 70, Progress in Astronautics and Aero- nautics, pp. 3-24, 1980.