H E A T F L O W I N T H E M E L T A N D C R U C I B L E F O R C R Y S T A L G R O W T H By R o w i n W i s a d i K h o B . A . Sc. , T h e Unive r s i ty of B r i t i s h C o l u m b i a , 1989 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F A P P L I E D S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S M E T A L S A N D M A T E R I A L S E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A August 1991 © R o w i n W i s a d i K h o , 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholariy purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Metals and Materials Engineering The University of British Columbia Vancouver, Canada Date August 22,1991 DE-6 (2/88) Abs t r ac t A magnetic field is imposed on the melt during Czochralski crystal growth to control fluid flow. The influence of an applied magnetic field on heat transfer in a germanium melt has been investigated in this work. In heat-transfer modelling of the Czochralski process, values of thermal conductivity of the crucible material, pyrolytic boron nitride, are required as a function of temperature. The thermal conductivity of a pyrolytic boron nitride crucible material has also been determined in this work. Both the investiga-tion of the influence of the applied magnetic field and the determination of the thermal conductivity have involved experimental temperature measurements and mathematical modelling of heat transfer. The finite element equations for the two-dimensional heat-conduction equation have been derived using the Galerkin method of weighted residuals. They have been validated by comparing the finite element solutions with the corresponding analytical solutions. The finite element results are in excellent agreement with the analytical solutions. A steady-state conduction-dominated mathematical model has been developed to an-alyze the temperature measurements obtained within the germanium melt in a Czochral-ski crystal growth configuration with and without an applied magnetic field of 0.099 tesla. The effect of the applied magnetic field on the heat transfer in the melt has been determined by fitting the model-calculated results to the measurements. The effective thermal conductivity of the melt has been found to decrease by a factor of seven due to the magnetic field. The thermal conductivity across a pyrolytic boron nitride crucible plate (through crucible thickness) has been determined from measurements of temperature responses i i in liquid gallium positioned on both sides of the plate, in conjunction with a transient mathematical model which simulates the thermal responses. By matching the simulated thermal responses with the measurements, the thermal conductivity of the pyrolytic boron nitride crucible plate has been obtained as a function of temperature. Table of Contents Abs t r ac t i i Table of Contents i v L i s t of Figures x L i s t of Tables x x Nomencla ture x x i i Acknowledgements x x i x Dedica t ion x x x 1 In t roduct ion 1 2 L i te ra ture Survey 5 2.1 Growth of Single Crystals by the Czochralski Process 5 2.2 Heat-Transfer Mathematical Modelling of the Czochralski Process . . . . 8 2.3 Effects of Magnetic Fields on the Czochralski Growth Process 14 2.4 Pyrolytic Boron Nitride as a Crucible Material for the Processing of Semi-conductors 16 2.4.1 Production, Crystal Structure, and Properties of Pyrolytic Boron Nitride 16 2.4.2 Applications of Pyrolytic Boron Nitride 17 iv 2.5 Experimental Determination of Thermal Conductivity 23 2.5.1 Steady-State Methods 25 2.5.2 Nonsteady-State Methods 26 2.5.3 Methods for Plate Specimens 27 2.6 Reported Thermal-Conductivity Values for Pyrolytic Boron Nitride . . . 29 3 F in i t e Element Solutions for the Two-Dimens iona l Hea t -Conduc t ion Equa t ion 33 3.1 Finite Element Equations 35 3.2 Comparisons with Analytical Solutions 44 3.2.1 Steady-State Conduction 44 3.2.1.1 One-Dimensional Solutions 45 3.2.1.2 Two-Dimensional Solutions 47 3.2.2 Steady-State Convection 58 3.2.3 Transient Conduction 59 3.2.4 Transient Convection 67 3.2.5 Discussion on the Finite Element Solutions 69 4 Effect of a Magne t i c F i e l d on Hea t Transfer i n a G e r m a n i u m M e l t 72 4.1 Temperature Measurements and Radial Temperature Gradients 72 4.2 Mathematical Model 83 4.2.1 Assumptions 83 4.2.2 Domain and Boundary Conditions 84 4.2.3 Finite Element Mesh of the Domain 86 4.2.4 Sensitivity Analysis of the Parameters 89 4.2.4.1 Side Surface Temperature, Twau 89 4.2.4.2 Ambient Temperature, Ta 91 v 4.2.4.3 Top Surface Biot Number, Bito 94 4.2.4.4 Bottom Surface Biot Number, Bi^ 94 4.2.4.5 Relative Sensitivities of the Parameters 94 4.3 Determination of the Magnetic Effect on the Effective Thermal Conduc-tivity of the Ge Melt 98 4.3.1 Estimation of the Biot Numbers 98 4.3.2 Insulated Bottom Surface 100 4.3.3 Non-Insulated Bottom Surface 104 4.3.4 Reduction Factor in the Effective Thermal Conductivity of the Ge Melt 112 4.4 Discussion on the Reduction Factor 113 4.5 Relative Thermal Resistance Analysis of the Mathematical Model . . . . 115 5 De te rmina t ion of the T h e r m a l C o n d u c t i v i t y of a P y r o l y t i c B o r o n N i -tr ide Cruc ib l e 119 5.1 Description of the Material 119 5.2 Methods of Measurement 120 5.2.1 Steady-State Method 120 5.2.1.1 Experimental Apparatus 122 5.2.1.2 Temperature Measurement Procedure 122 5.2.1.3 Preliminary Experimental Results and Discussion . . . . 125 5.2.1.4 Problems Associated with the Steady-State Method . . . 128 5.2.2 Transient Method 130 5.2.2.1 Experimental Apparatus 130 5.2.2.2 Temperature Measurement Procedure 132 5.2.2.3 Test Materials and Specifications 133 vi 5.2.2.4 Experimental Results and Discussion 133 5.2.2.4.1 Measured Results 136 5.2.2.4.2 Comparison of the Results for the Three Differ-ent Materials 138 5.3 Mathematical Model 142 5.3.1 Assumptions 143 5.3.2 Domain and Boundary Conditions 144 5.3.3 Finite Element Mesh of the Domain 146 5.3.4 Results of Sensitivity Analysis and Discussion 146 5.4 Model Analysis of the Experimental Results 153 5.4.1 Corrections and Utilization of the Temperature Measurements . . 155 5.4.1.1 Corrections for the Thermocouple Errors and Furnace Temperature Gradient 155 5.4.1.2 Utilization, of the Thermocouple Readings 157 5.4.2 Calibrations of the Overall Heat-Transfer Coefficient 157 5.4.2.1 Assumed Function of the Overall Heat-Transfer Coefficient 163 5.4.2.2 Discussion on the Discrepancy between the Experimental and Calculated Curves 163 5.4.2.3 Model Predictions Assuming Zero Overall Heat-Transfer Coefficient 165 5.4.2.4 Fitting the Temperature-Time Curves of Pyrex 166 5.4.2.5 Consistency of C\ and C2 Constants 166 5.4.2.6 - Fitting the Temperature-Time Curves of B N 174 5.4.2.7 Fitted Equations for d and C2 177 5.4.3 Thermal Diffusivity Fitting of the Temperature-Time Curves of the P B N Crucible 177 vii 5.4.4 Sensitivity of CXPBN, <*Ga, h and t 192 5.4.5 Sources of Errors in the Derived Thermal DifFusivity 192 5.5 Thermal Conductivity of the P B N Crucible 198 5.6 Comparison of the Thermal Conductivities 201 5.7 Relative Thermal Resistance Analysis of the Mathematical Model . . . . 204 6 S u m m a r y and Conclusions 209 7 Recommendat ions for Future W o r k 211 A Derivat ions of the F in i t e Element Equat ions 212 A . l Cylindrical Coordinates (r, z) 212 A.1.1 Type of Element and Shape Functions 215 A.1.2 Mapping Functions and Numerical Integration 217 A . 1.2.1 Conduction Matrix 221 A.1.2.2 Convective/Radiative Matrices and Vectors 223 A.1.2.3 Capacitance Matrix 225 A . 1.3 Time-Stepping Techniques 225 A.1.3.1 Two-Point Recurrence Scheme 226 A . 1.3.2 Three-Point Recurrence Scheme 229 A . 2 Cartesian Coordinates (r, z) 232 B A n a l y t i c a l Solutions 234 B . l Steady-State Conduction 234 B . l . l One-Dimensional Solution 234 B. l .2 Two-Dimensional Solutions 235 B . l.2.1 Boundary Temperature Profile 1 236 B.1.2.2 Boundary Temperature Profile 2 236 viii B.2 Steady-State Convection - Two-Dimensional Solution 236 B.3 Transient Conduction - One-Dimensional Solutions 237 B.3.1 Case 1 237 B.3.2 Case 2 238 B.4 Transient Convection. - One-Dimensional Solutions 239 B.4.1 Case 1: Cartesian System 239 B.4.2 Case 2: Cylindrical System 239 C Ca lcu la t ion of the Ra te of Temperature Increase i n the L i q u i d G a l l i u m inside the M i d d l e Compar tment A s s u m i n g N o Hea t Loss 241 D Ca lcu la t ion of the Rate of Temperature Increase i n the L i q u i d G a l l i u m inside the Side Compar tment A s s u m i n g U n i - D i r e c t i o n a l Hea t F low across the Specimen 242 B ib l iog raphy 243 ix Lis t of Figures 2.1 Schematic drawing of a CZ pulling apparatus [86] 6 2.2 Schematic heat flow during the CZ crystal growth [56] 9 2.3 Schematic diagram of the P B N deposition process [10] 18 2.4 Articles of P B N [55] 19 2.5 Hexagonal crystal structure of boron nitride [25] 20 2.6 Thermal conductivity of P B N reported in reference [25] for specimen B and C 30 2.7 Thermal conductivity of P B N reported in reference [25] for specimen C . 31 3.8 Element type: (a) 8-noded quadratic isoparametric element; (b)Quadratic isoparametric boundary element 38 3.9 Mapping of the quadratic isoparametric element into the (£, 77) coordi-nates: (a) domain element; (b) boundary element 40 3.10 Comparison of finite element solutions with 1-dimensional analytical solu-tions for steady-state conduction using r,- = 0.01 m and ra = 1.61 m . . . 46 3.11 Comparison of finite element solutions with 1-dimensional analytical solu-tions for steady-state conduction using r,- = 1000.01 m and ra = 1001.61 m 48 3.12 Comparison of finite element solutions with 1-dimensional analytical solu-tions for steady-state conduction using r, = 0.0001 m and r0 = 1.6001 m 49 3.13 Comparison of finite element solutions with 1-dimensional analytical solu-tions for steady-state conduction using rt- = 0.0001 m and rQ = 3.2001 ra 50 x 3.14 Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state conduction using temperature profile 1 along the z axis 53 3.15 Comparison of finite element solutions with 2-dimensional analytical so-lutions for steady-state conduction using temperature profile 1 along the midline between the z axis and boundary 54 3.16 Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state conduction using temperature profile 2 along the z axis 55 3.17 Comparison of finite element solutions with 2-dimensional analytical so-lutions for steady-state conduction using temperature profile 2 along the midline between the z axis and boundary 56 3.18 Comparison of finite element temperature distributions with 2-dimensional analytical solutions for steady-state conduction using temperature profile 2: (a) 2 x 4 elements; (b) 4 x 8 elements; (c) 8 x 16 elements; (d) analytical 57 3.19 Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state convection along the z axis 60 3.20 Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state convection along the midline between the z axis and the boundary 61 3.21 Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state convection along the boundary 62 3.22 Comparison of finite element temperature distributions with 2-dimensional analytical solutions for steady-state convection: (a) 2x 4 elements; (b) 4x8 elements; (c) 8 x 16 elements; (d) analytical 63 XI 3.23 Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient conduction at x — 0 and x = 1 m 65 3.24 Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient conduction at x = 0.01 m 66 3.25 Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient convection at x — 1.0 m 68 3.26 Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient convection at r = 0.001 m 70 4.27aMeasured temperatures at 5 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 74 4.27bMeasured temperatures at 10 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 75 4.27cMeasured temperatures at 30 ram below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 76 4.27dMeasured temperatures at 40 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 77 4.28aMeasured temperatures at 5 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] . . 78 4.28bMeasured temperatures at 10 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] . . 79 4.28cMeasured temperatures at 30 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] . . 80 4.28dMeasured temperatures at 40 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] . . 81 xu 4.29 Radial temperature gradients in the Ge melt at the four melt depths in a crucible rotating at 14.4 rpm with and without an applied magnetic field of 0.099 T [87] 82 4.30 Idealized Domain of the mathematical model for investigating the influence of the applied magnetic field on the effective thermal conductivity of the Ge melt 85 4.31 Effects of finite element mesh on temperature field: (a) 2 x 3 elements; (b) 4 x 6 elements; (c) 8 x 12 elements; (d) 12 x 18 elements [unit in °C] . . 87 4.32 Finite element mesh of the mathematical model for investigating the in-fluence of the applied magnetic field on the effective thermal conductivity of the Ge melt 88 4.33 Effects of constant Twau\ (a) reference; (b) TwaU = 1050 [unit in °C] . . . 90 4.34 Effects of linear TwaU: (a) reference; (b) linear TwaU: 1025-975 [unit in °C] 92 4.35 Effects of Ta: (a) reference; (b) Ta = 800 [unit in °C] 93 4.36 Effects of Bito: (a) reference; (b) Bito = 0.115 [unit in °C] 95 4.37 Effects of Biho- (a) reference; (b) = 0.0575 [unit in °C] 96 4.38 Relative effects of the model parameters (temperature profiles at the mid-dle height of the melt) 97 4.39 Comparison of the measured and calculated temperature profiles at 5, 10, 30 and 40 mm below the melt surface under no magnetic field assuming bottom insulation 101 4.40 Comparison of the measured and calculated radial temperature gradients as a function of melt depth under no magnetic field assuming bottom insulation 103 X l l l 4.41 Comparison of the measured and calculated temperature profiles at 5, 10, 30 and 40 mm below the melt surface under an applied magnetic field of 0.099 T assuming bottom insulation 105 4.42 Comparison of the measured and calculated radial temperature gradients as a function of melt depth under an applied magnetic field of 0.099 T assuming bottom insulation 106 4.43 Comparison of the measured and calculated radial temperature gradients as a function of melt depth with and without an applied magnetic field of 0.099 T assuming bottom insulation 107 4.44 Comparison of the measured and calculated temperature profiles at 5, 10, 30 and 40 mm below the melt surface under no magnetic field without bottom insulation 108 4.45 Comparison of the measured and calculated temperature profiles at 5, 10, 30 and 40 mm below the melt surface under an applied magnetic field of 0.099 T without bottom insulation 109 4.46 Comparison of the measured and calculated radial temperature gradients as a function of melt depth with and without an applied magnetic field of 0.099 T for both the insulating and non-insulating bottom boundary conditions I l l 4.47 Thermal resistance circuit in the model domain used for investigating the influence of the magnetic field on the effective thermal conductivity of the Ge melt 116 4.48 Thermal circuit with relative resistances of the mathematical model used for investigating the influence of the magnetic field on the effective thermal conductivity of the Ge melt 117 xiv 5.49 A broken edge of the P B N crucible revealing the layers that make up its thickness (magnification = 16 x) 121 5.50 Boron nitride boat used in the experimental apparatus 123 5.51 Schematic diagram of the experimental apparatus for the steady-state measurements 124 5.52 Temperatures as a function of time for T C I , TC2 and TC3 positioned in the three cells as shown in Figure 5.51 126 5.53 Rates of temperature increase as a function of time for T C I , TC2 and TC3 positioned in the three cells as shown in Figure 5.51 127 5.54 Schematic diagram of the modified experimental apparatus for the tran-sient measurements 131 5.55 Typical output of the four thermocouple readings (taken from experiment no. 9 of Pyrex) 137 5.56 Temperature responses for Pyrex (taken from experiment no. 9 of Pyrex) 139 5.57 Temperature responses for B N (taken from experiment no. 2 of BN) . . . 140 5.58 Temperature responses for the P B N crucible in the thickness direction (taken from experiment no. 10 of the P B N crucible) 141 5.59 Idealized Domain of the Mathematical Model for the determination of the thermal conductivity of the P B N crucible 145 5.60 Finite element mesh of the model domain for the determination of the thermal conductivity of the P B N crucible 147 5.61 Time-varying input temperature profile used for the sensitivity analysis (taken from experiment no. 14 of the P B N crucible) 148 5.62 Effects of the thermal diffusivity of specimen 150 5.63 Effects of the thermal diffusivity of liquid gallium 151 5.64 Effects of the overall heat-transfer coefficient 152 xv 5.65 Effects of the specimen thickness 154 5.66 Corrected thermocouple readings 156 5.67 Reference flowchart for the procedure used to obtain the thermal conduc-tivity of the P B N crucible 158 5.68 Thermal diffusivity of Pyrex [88] used for model calculations in the fitting process (equation 5.49 ) 159 5.69 Thermal diffusivity of B N [88] used for model calculations in the fitting process (equation 5.50 ) 160 5.70 'Stagnant' and effective thermal diffusivity of liquid gallium [88] used for model calculations in the fitting process (equation 5.53 ) . . . . . . . . . 162 5.71 Typical fit between the experimental and calculated temperature-time curves (taken from experiment no. 9 of Pyrex) 164 5.72 Comparison of the model predictions with experimental curves assuming zero overall heat-transfer coefficient (insulated boundary) 167 5.73aTime-varying input temperature profiles of Pyrex used for model calcula-tions (experiment no. 1 to experiment no. 5) 168 5.73bTime-varying input temperature profiles of Pyrex used for model calcula-tions (experiment no. 6 to experiment no. 9) 169 5.74aComparisons of the experimental and calculated curves at the locations of T C I and TC2 for Pyrex (experiment no. 1 to experiment no. 5) 170 5.74bComparisons of the experimental and calculated curves at the locations of T C I and TC2 for Pyrex (experiment no. 6 to experiment no. 9) 171 5.75 Comparisons of the experimental and calculated curves at the locations of T C I and TC2 for B N using C\ and C 2 values from Pyrex for normal fc^/v and reduced k s N (0.5 x normal ksN) • 172 xvi 5.76 Normal and reduced (0.5 x normal) thermal conductivity of B N [88] (equa-tion 5.51 ) 173 5.77aTime-varying input temperature profiles of B N used for model calculations (experiment no. 1 to experiment no. 3) 175 5.77bTime-varying input temperature profiles of B N used for model calculations (experiment no. 4 to experiment no. 5) 176 5.78aComparisons of the experimental and calculated curves at the locations of T C I and TC2 for B N (experiment no. 1 to experiment no. 3) 178 5.78bComparisons of the experimental and calculated curves at the locations of T C I and TC2 for B N (experiment no. 4 to experiment no. 5) 179 5.79 C\ as a function of ambient temperature (combined Pyrex and B N values) (equation 5.55 ) 180 5.80 C 2 as a function of ambient temperature (combined Pyrex and B N values) (equation 5.56 ) 181 5.81aTime-varying input temperature profiles of the P B N crucible used for model calculations (experiment no. 1 to experiment no. 4) 183 5.81bTime-varying input temperature profiles of the P B N crucible used for model calculations (experiment no. 5 to experiment no. 8) 184 5.81cTime-varying input temperature profiles of the P B N crucible used for model calculations (experiment no. 9 to experiment no. 12) 185 5.81dTime-varying input temperature profiles of the P B N crucible used for model calculations (experiment no. 13 to experiment no. 16) 186 5.82aComparisons of the experimental and calculated curves at the locations of T C I and TC2 for the P B N crucible (experiment no. 1 to experiment no. 4)187 5.82bComparisons of the experimental and calculated curves at the locations of T C I and TC2 for the P B N crucible (experiment no. 5 to experiment no. 8)188 xvn 5.82cComparisons of the experimental and calculated curves at the locations of T C I and TC2 for the P B N crucible (experiment no. 9 to experiment no. 12) 189 5.82dComparisons of the experimental and calculated curves at the locations of T C I and TC2 for the P B N crucible (experiment no. 13 to experiment no. 16) 190 5.83 Thermal diffusivity of the P B N crucible in the thickness direction obtained from the fitting process (equations 5.57 and 5.58 ) 191 5.84 Effects of a P B N 193 5.85 Effects o f o G a 194 5.86 Effects of h 195 5.87 Effects of t 196 5.88 Specific heat capacity of P B N [2] (equation 5.59 ) 199 5.89 Thermal conductivity of the P B N crucible in the thickness direction cal-culated using apBNi PPBN and C P P B N (equations 5.61 and 5.62 ) 200 5.90 Comparison of the thermal conductivities of Pyrex [88] , B N [88] , P B N in the "a" direction [2] , P B N in the "c" direction [2] , and the P B N crucible in the thickness direction (this work) (equations 5.63 , 5.51 , 5.64 , 5.65 , and 5.61 and 5.62 , respectively) 203 5.91 Thermal resistance circuit in the model domain used for the determination of the thermal conductivity of the P B N crucible 205 5.92 Thermal circuit with relative resistances of the mathematical model used for the determination of the thermal conductivity of the P B N crucible . 208 A.93 Element type: (a) 8-noded quadratic isoparametric element; (b)Quadratic isoparametric boundary element 216 xviii A.94 Mapping of the quadratic isoparametric element from the (r, z) coordi-nates into the (£, rj) coordinates: (a) domain element; (b) boundary element219 A.95 Linear shape functions (two-point recurrence formulae) for the time do-main [57] 227 A.96 Quadratic shape functions (three-point recurrence formulae) for the time domain [57] 230 xix List of Tables 2.1 Properties of P B N : (a) from reference [10] ; (b) from reference [2] : . . 21 2.2 Thermal conductivity of P B N reported in reference [2] 32 2.3 Thermal conductivity of P B N reported in reference [10] 32 3.4 Different pairs of values of r; and ra used for comparisons of steady-state conduction finite element results with 1-dimensional analytical solutions 45 3.5 Two different temperature profiles along the axial boundary for compar-isons of steady-state conduction finite element results with 2-dimensional analytical solutions 51 3.6 Size of domain, boundary conditions and heat-transfer properties used for comparison of steady-state convection finite element results with 2-dimensional analytical solutions 59 4.7 Reference parameter values used in the sensitivity analysis 89 4.8 Parameter values used in the calculations for the normal case 110 4.9 Parameter values used in the calculations for the magnetic case 110 4.10 Biot numbers used in fitting the measurements assuming bottom insulationll2 4.11 Biot numbers used in fitting the measurements accounting for bottom heat loss 112 5.12 Dimensions of the test samples 133 5.13 Adjusted furnace (ambient) temperatures for the Pyrex test sample . . . 134 5.14 Adjusted furnace (ambient) temperatures for the B N test sample 134 xx 5.15 Adjusted furnace (ambient) temperatures for the test sample of the P B N crucible in the thickness direction 135 5.16 Reference values used in the sensitivity analysis of the model parameters 146 5.17 C\ and C% values for Pyrex at different ambient temperatures 174 5.18 C\ and C2 values for B N at different ambient temperatures 174 5.19 Comparison of the thermal conductivity values of P B N in the "a" and "c" directions with the values obtained in this investigation for the P B N crucible in the thickness direction 202 A.20 Sampling locations and weighting coefficients of the Gaussian-quadrature formula 220 xxi Nomencla ture L a t i n S y m b o l A area b radial length B magnetic flux density Bi Biot number Bifnagnet Biot number used under an applied magnetic field Biot number used under no magnetic field C axial length C heat capacity per unit volume cp specific heat capacity C i constant used for calibrating the overall heat-transfer coefficient c2 constant used for calibrating the overall heat-transfer coefficient d thickness det determinant f function F radiation view factor H) function h overall heat-transfer coefficient hc convective heat-transfer coefficient K radiative heat-transfer coefficient H h k xxu i,j unit vectors / content of inverse of Jacobian matrix IQ() modified Bessel function of order 0 of the first kind J 0 () Bessel function of order 0 of the first kind Ji () Bessel function of order 1 of the first kind k thermal conductivity; thermal decay constant kmmagnet effective thermal conductivity of Ge melt under an applied magnetic field kmnormal effective thermal conductivity of Ge melt under no magnetic field Ka thermal conductivity in basal plane / mean free path L length m number of sampling points in one direction; mass; number of elements axially; highest order of derivative n number of nodes in element; unit vector normal to surface; number of elements radially N trial (shape) function N time derivative of shape function Nz derivative of shape function with respect to z q heat (thermal energy) qcond conductive heat qconv convective heat qrad radiative heat qin heat input Q rate of heat generation or consumption; magnitude of thermal energy flux per unit area xxiii 9 rate of heat flow convective rate of heat loss radiative rate of heat loss 9in rate of heat input <?out rate of heat output r,z coordinates in radial and axial directions, respectively R residual; thermal resistance; radius t time; thickness time at half maximum temperature value T temperature Ta ambient temperature Tar ambient temperature for radiation Ts surface temperature Twau temperature at the side surface of Ge melt TCI thermocouple no. 1 TC2 thermocouple no. 2 TC3 thermocouple no. 3 TC4 thermocouple no. 4 Too ambient temperature T O O l ambient temperature for convection T •*• 0 0 2 ambient temperature for radiation V velocity V volume W weighting function coordinates in horizontal and vertical directions, respectively xxiv [B] B matrix [C] global capacitance matrix [Cv] global convective capacitance matrix {/} global thermal load (force) vector I M global convective vector Ur} global linearized radiative vector [J] Jacobian matrix [K] global thermal stiffness matrix [Kc] global conduction matrix [KH] global convection matrix [Kr] global linearized radiation matrix m global temperature vector {f} global vector of time derivatives of nodal temperatures At time step size AT temperature difference Ax distance between 2 points of temperature measurements M a t h e m a t i c a l Symbo l d differential I integral V differential operator [ ] square or rectangular matrix [ f transpose of matrix [ r1 inverse of matrix {} column or row vector X X V { } transpose of vector Greek S y m b o l a thermal diffusivity /3 root 7] root 6 point spacing 7 parameter in time-stepping scheme T domain boundary e emissivity A boundary Jacobian 7T 3.1415926536 9 parameter in time-stepping scheme p density a Stefan-Boltzman constant S summation u angular frequency £, f] natural local coordinates in horizontal and vertical directions, respec-tively C dimensionless time variable Subscr ipt bo bottom surface i inner i,j,l,m,n indices m melt max maximum xxvi mid middle height o outer r, z in r and z directions, respectively s specimen to top surface x,y in x and y directions, respectively £,77 in £ and 77 directions, respectively Superscr ipt (e) elemental A c r o n y m A S T M American Society for Testing and Materials CZ Czochralski L E C liquid encapsulated Czochralski M B E molecular beam epitaxy M O C V D metal organic chemical vapor deposition SI semi-insulating V G F vertical gradient freeze U n i t kg kilogram °K degree Kelvin m meter mm milimeter rpm revolution per minute s second T tesla xxvii W watt °C degree Celcius Acknowledgements I would like to extend my thanks and gratitude to Dr. Fred. Weinberg and Dr. Carlos Schvezov for their suggestions, support and encouragement over the course of this work. The emotional support given by Carlos is worth a special note. The help of Chris Par-feniuk on many aspects of this project is greatly appreciated. Thanks are also extended to many fellow graduate students and research engineers in the department, particularly Bernardo Hernandez and Willie Fajber for their technical assistance. I would like to take this opportunity to express my indebtedness to my parents and grandmother for their sincere love and caring, support and encouragement that they always have for me. I would also like to thank two very good friends, Patsy and Seniyanti, for their caring, support and encouragement. The friendship from Irawan Soeharjono, who has provided me with various entertaining activities, is noteworthy. Finally, I would like to acknowl-edge the Natural Sciences and Engineering Research Council of Canada for the financial support. xxix Dedica t ion I dedicate this thesis to my parents, brothers and sister. xxx Chapter 1 Introduction The thermal system in a crystal growth environment in which a semiconductor crystal is pulled from its melt is complex. Heat transfer in solid, liquid and gas phases has to be considered, complicated by fluid flow in the melt and gas due to natural and forced convection. The quality of the crystal produced is related to the thermal system. For example, thermal stresses resulting from temperature gradients in the crystal can generate dislocations which are detrimental to the crystal quality. In this investigation, two aspects of the thermal system are considered, heat transfer in a germanium melt in which fluid flow is occurring, and heat flow through the wall of a pyrolytic boron nitride crucible. Characterization of the heat flow in the melt and crucible requires thermal properties of these materials. Of particular importance is thermal conductivity, which is defined as "the heat flow across a surface per unit area per unit time, divided by the negative of the rate of change of temperature with distance in a direction perpendicular to the surface". It is also known as coefficient of conductivity or heat conductivity [74]. As the demand for reliable microelectronic devices continues to increase and the pro-cess of miniaturization pushes the physical size of the devices ever smaller, the quality requirements of semiconductor materials have become ever more stringent. The disloca-tion density in the crystal must be low, the crystal must have uniform solute distribution and shape, be striation-free, and the impurity in the crystal must be at a controlled level, to name but a few. As the properties and quality of the crystal are closely related to the processing conditions, meeting such stringent requirements is a major driving force 1 2 in improving the crystal growth processes, particularly Czochralski (CZ) growth which has been used successfully to grow elemental and compound semiconductors such as Si, GaAs and InP. The improvement of the CZ process has been mostly based on the method of trial and error. The fundamental principles of the process are not well understood quantitatively. As a consequence, the extrapolation of past empirical techniques for the CZ growth of silicon to the growth of gallium arsenide by the liquid encapsulated Czochralski (LEC) process has yielded poor results. Accordingly, intensive theoretical and experimental studies have been carried out in trying to understand the process better. One of the most extensive areas of research in crystal growth is heat transfer since it is a basic cornerstone in understanding the overall growth process. As a consequence of the experimental problems encountered in the hostile thermal and chemical environments present in the growth system and as an alternative way to understand the heat transport phenomena in the system, a large number of mathematical models of heat transfer has been developed to gain insight into certain areas of the complex CZ system and to link the final crystal properties to the processing conditions. They have ranged from simple analytical relationships to system-oriented models which require sophisticated numerical techniques and large computer memory space. Numerical techniques such as the finite element method and the Newton-Raphson method have been extensively used to obtain approximate solutions to the governing partial differential equations, and supercomputers are often utilized to solve for large systems of matrix equations. Another major area of investigation is the effect of a magnetic field on the process. As a way of obtaining better quality crystals, a magnetic field is increasingly being applied in the CZ process. Mathematical modelling has also played an important role in this area of research. The type of crucible material used is crucial to the quality as well as the properties 3 of the crystal. The effects of crucible material have also been studied quite extensively. New materials resulting in better crystal quality are replacing existing crucible materials. A material which has been increasingly used for the processing of semiconductors is pyrolytic boron nitride. The research carried out on the CZ process has resulted in better system design and improved the quality of the crystals. However, a thorough quantitative description of the process still requires appreciable fundamental work. Aside from the complex science and transport phenomena involved in modelling, a major problem is the lack of heat-transfer properties to characterize the CZ system. The heat-transfer properties are the thermal properties of the materials involved in the system and the heat-transfer coefficients necessary to build the mathematical models. Moreover, the reliability of the values which are available for some of these properties is being questioned [73]. Hence, obtaining these properties, particularly the more important ones such as thermal conductivity, has a high priority in this field of research, especially when new materials like pyrolytic boron nitride are introduced into the growth system. The need to measure the thermal conductivities of the various constituents relevant to the CZ growth system is clear. It is also equally important to examine the change in thermal conductivity due to external effects such as that of a magnetic field. Accordingly, the objectives of this work are to investigate the influence of an applied magnetic field on the effective thermal conductivity of a germanium melt in a CZ growth configuration, and to determine the thermal conductivity of a pyrolytic boron nitride crucible. Unfortunately, no simple standard methods to measure the thermal conductivities are available. There is no known standard method to examine the effect of a magnetic field on the effective thermal conductivity of a melt. Standard methods to obtain the thermal conductivity of the material, especially at high temperatures, require expensive or sophisticated experimental apparatus which is not readily obtainable or built in a 4 limited amount of time. Consequently, a method for investigating the effect of a magnetic field on the effective thermal conductivity must be devised and an alternative route to obtain the thermal conductivity must be taken. A way to achieve these goals is to utilize a combination of mathematical modelling and experimental measurements. In either case, a mathematical model of heat transfer is employed in the analysis of experimental temperature measurements. By matching the model-calculated results with the measurements, the influence of an applied magnetic field on the effective thermal conductivity of the germanium melt will be investigated and the thermal conductivity of the pyrolytic boron nitride crucible in the thickness direction as a function of temperature will be determined. Chapter 2 Literature Survey This thesis touches several essential topics such as growth of single crystals by the Czochralski process, heat-transfer mathematical modelling of the process, effects of mag-netic fields on the process, pyrolytic boron nitride as a container material for the pro-cessing of semiconductors, and experimental determination of thermal conductivity. A literature review of each of these topics will be given in separate sections of this chapter. Reported values for pyrolytic boron nitride materials are included. 2.1 Growth of Single Crystals by the Czochralski Process There are several methods to grow single crystals. These include the Bridgman, the floating zone, and the Czochralski method. The most significant commercial method is the Czochralski (CZ) process. The CZ growth process is generally more complex than the other growth techniques. A typical CZ crystal pulling apparatus is shown in Figure 2.1 [86]. The CZ growth system consists of a crucible supported by a pedestal which can be rotated and moved vertically during crystal growth. The crucible is surrounded by a graphite resistance heater or a radio-frequency coil with a graphite susceptor acting as a heat source. A pull rod, which also travels vertically, is located above the center of the crucible. This pull rod is used to lower the seed to the surface of the melt at the start of the process and subsequently to pull the crystal from the melt. Heat shields are sometimes placed above the crucible to achieve the desired radiative heat-transfer characteristics. The entire 5 front opening door front opening chamber ;ure 2.1: Schematic drawing of a CZ pulling apparatus 7 assembly is housed in a large temperature-controlled enclosure. The procedure of growing a crystal starts by filling the crucible with high purity polycrystalline material and often a small quantity of impurity or dopant, which alters the electronic properties of the grown crystal. The enclosure is filled with an inert atmosphere such as argon or nitrogen, and the crucible is heated to melt the material. A precisely oriented single crystal seed attached to the pull rod is lowered to the surface of the melt. After the seed and melt have reached thermal equilibrium, the seed is slowly withdrawn from the melt and a single crystal is pulled upward. A meniscus, whose shape is determined by the action of surface tension against the force of gravity, forms about the edge of the crystal by surface tension as the melt is pulled upward. The crystal is usually grown through a "necking" operation during the beginning of growth where a very small diameter crystal of some length is solidified in an attempt to eliminate dislocations. After necking the crystal is grown to the desired diameter by reducing the power output of the heater, and subsequently, system parameters such as the growth rate and heater output are adjusted to achieve a constant-diameter crystal. A n active control mechanism, often automatic, is employed to keep the crystal diameter constant. At the end of the run the crystal is withdrawn from the melt and slowly cooled. A variant of the CZ process is called liquid encapsulated Czochralski (LEC) growth. In the L E C process, an inert substance, usually boric oxide ( B 2 O 3 ) , is floated on the surface of the melt to prevent the escape of volatile components from the melt. This is important in the growth of III-V compounds such as gallium arsenide and indium phosphide, where the stoichiometric ratio of the elements must be precisely controlled in order to produce a high-quality crystal. 8 2.2 Heat-Transfer Mathematical Modelling of the Czochralski Process Heat flow in the CZ furnace is complex, as illustrated schematically in Figure 2.2 [56]. In the L E C growth of compound semiconductors, an additional degree of complexity comes from the encapsulant layer. To do a quantitative heat-transfer analysis of the CZ and L E C growth, the thermal variables must be tightly controlled. The experimental difficul-ties encountered in the growth system make the CZ process conducive to mathematical modelling. Extensive mathematical models have been developed in an attempt to understand the heat transfer in the growth system. The models contain necessary assumptions to render the problem tractable, especially the L E C growth in which the encapsulant layer plays a crucial role in the heat transfer of the process. A relationship between crystal radius and pull rate was first derived by Bill ig [24] by solving an one-dimensional approximation for heat transfer in the crystal with no heat transfer from the melt. He postulated that the radius of a steadily growing crystal should be inversely proportional to the square of the steady-state growth rate of the crystal. A n approximate analytical representation of the shape of the meniscus in the CZ growth which takes into account the angle of contact at the three-phase boundary, the crystal radius and the magnitude of the surface tension was derived by Hurle [22]. Predicting the shape of the melt-crystal interface together with the temperature field in the CZ process is one of the major challenges in modelling. The shape of the in-terface influences thermal stress in the crystal during growth. The larger the interface deviates from planarity, the larger the thermal stress develops in the crystal. Such vari-ations of thermal stresses with the deviation from planarity have been established by numerical techniques [41]. Wilcox and Duty were the first to couple the calculation of the heat transfer and melt/solid interface shape for a steadily growing CZ crystal [85]. t J 1 J 9cond Figure 2.2: Schematic heat flow during the CZ crystal growth [56] 10 Their numerical approach solved for the temperature field and interface shape with the assumption that the heat transfer between the crystal and melt could be represented by a heat-transfer coefficient. Ramachandran et al. [59] developed a model based on an iterative finite element scheme for prediction of the temperature distribution in the crystal for an assumed melt temperature and position of the melt-crystal interface which accounts for the detailed radiative heat transfer between the various surfaces present in the crystal pulling apparatus. Srivastava et al. then extended the model to solve for the temperature field in the melt as well as in the crystal, assuming the heat transfer in the melt occurs by conduction only while simultaneously determining the shape of the melt-gas meniscus as well as the melt-crystal interface for a fixed crystal radius [69]. H . Kopetsch developed a finite difference algorithm for time-dependent incompressible viscous flow which includes the dynamics of the solid-liquid interface in the CZ crystal growth process [30]. A n additional degree of freedom to control the interface shape and the diameter of the crystal is by using jet cooling. The effects of cooling the crystal side surface by blowing a jet of an inert gas have been examined by Srivastava et al. [70]. It was found that the interface shape and the crystal diameter can be controlled by adjusting the gas flow rate through the jet. In order to quantitatively predict the temperature field in the Czochralski process, a precise knowledge of the heat transfer taking place in the entire furnace is necessary. System-oriented models for optimizing the design of particular growth systems which in-clude the effects of furnace heat shields and temperature distributions were developed by Williams and Reusser [29] and Matsumoto et al. [44]. These works assumed fixed inter-faces and required experimentally measured thermal boundary conditions. Heat-transfer modelling of the entire furnace was first done by Dupret et al. using external power in-put and convective losses on the enclosure as the boundary conditions [26]. Their model 11 allows the investigation, of the sensitivity of the temperature distribution with respect to the geometry of the furnace, the presence of a heat shield and the material properties of the constituents. The effects of geometrical parameters upon the L E C growth of GaAs crystals have been analyzed by Nicodeme et al. using a global heat-transfer model which entirely relies upon external control parameters such as the pulling rate and the power input [62,63]. The model takes full account of the geometrical configuration and the ma-terial properties of a L E C furnace. They were able to obtain a good agreement between the experimental and theoretical shapes of the interface during growth. To determine the interplay between the geometry of the crystal growth system and heat transfer, several workers have used complex moving-boundary mathematical mod-els in order to capture the coupled effects of heat transfer and capillarity. The most extensive modelling in this area has been the work of R. A . Brown and his co-workers at Massachusetts Institute of Technology. The dynamics of the heat transfer, the stability and control of the shapes of the interfaces in the CZ and L E C growth systems have been studied in great details with thermal-capillary models. They have used the finite ele-ment method in conjunction with the Newton-Raphson method to solve for the resulting system of non-linear equations [38,36,48,37,61]. Major parts of the growth process have been captured in integrated mathematical models [21]. Another area which involves extensive mathematical modelling is fluid flow in the melt. Flow and temperature measurements in the melt during crystal growth are dif-ficult because of the hostile high temperature and high pressure environment. Thus, besides using model experiments with liquids having similar properties as those of semi-conductors to provide better conditions, mathematical models have been very helpful for investigating the transport phenomena and the effects of externally applied mechanisms on transport processes during crystal growth. 12 Transport processes in the bulk liquid play an essential role in the quality of the re-sulting crystal. To study the transport phenomena in the moving bulk, combined heat and fluid flow modelling has been carried out by several researchers. One such researcher who has modelled extensively the Czochralski bulk flow is W. E. Langlois. He along with his co-workers has done simulations ranging from simple parameter sensitivity study [80] and melt convection in a magnetic field [81] to effects of Joule heating and finite electri-cal conductivity of the crystal on hydromagnetic Czochralski flow [83,84]. Crochet et al. [50] formulated steady-state and time-dependent finite element simulations to verify the basic mechanisms of Czochralski melt convection and to explore the combined effects of such mechanisms. Mihelcic et al. [53] attempted a three-dimensional simulation of the Czochralski bulk flow in the presence of a non-axisymmetrical temperature distribution at the crucible wall; they showed that crucible rotation is necessary in order to sym-metrize the temperature field in the vicinity of the solid-liquid interface. The transport phenomenon near the interface of a CZ grown crystal has also been theoretically analyzed by Balasubramaniam and Ostrach [67]. The study of Czochralski bulk flow has also been carried out by several other investi-gators. A theoretical study has been done by Tsukada et al. on the flow and temperature fields in CZ single crystal growth [77]. Numerical studies of heat transfer including free and forced convection in the L E C growth of GaAs have been carried out by Sabhapathy et al. [64]. Sackinger, Brown and Derby have also extended the thermal-capillary model to include fluid flow in the melt [60]. However, there are problems preventing extensive calculations of fluid flow in the melt [72]. One problem is the sensitivity of the results to the details of the heat-transfer environment, which may make the calculations meaningless. Another problem is that the database of the thermophysical properties of the materials is not accurate enough to support this level of analysis. Detailed simulation of the melt flow is presently at the 13 forefront of numerical modelling [73]. Many of the thermal models were utilized to calculate the corresponding resolved shear stresses in the crystal to estimate the dislocation density. A seminal paper by Jordan et al. established that the generation of dislocations in a GaAs crystal is related to the thermal stresses in the crystal by comparing the thermally-induced stress field obtained analytically with the experimental dislocation field [9]. Following this work, similar modelling techniques were applied to other semiconductors such as Si and InP [8,7,6]. Along with these modelling efforts, Jordan compiled and estimated the material properties and the transport coefficients of the growth system [3,4,5]. In a series of papers by Schvezov et al., the temperature distribution and the stress field within a GaAs crystal have been calculated using the finite element method [18,19]. They examined the effects of crystal radius fluctuations and different growth conditions on the shear stress distribution in the crystal [20,16] as well as the thermal and stress fields in the crystal after separation from the melt [17]. Heat-transfer modelling has played a profound role in understanding the growth pro-cess. Because of many simplifying assumptions that must be made to permit a closed-form analysis and the tremendous complexity in the growth system, accurate quantitative predictions are difficult to obtain. The limitations of today's modelling also He in the fact that the available computing power is not able to completely model the transport phenomena in the overall process. Even with the most powerful supercomputers and the most efficient algorithms, numerical modelling at this moment cannot capture all of the processes important in an industrial scale solidification system [72]. It is crucial to know the values for the thermophysical properties and transport co-efficients, as well as the dependence of these properties and coefficients on temperature and pressure [73]. However, the information for this database is limited and inconsis-tent. For example, the value of the coefficient of thermal expansion for molten silicon is 14 only available from two isolated measurements, which disagree by a factor of three [73]. The difference in these values is critical in flow simulations [80]. The values for GaAs and InP have been reported by Jordan [4], but they are now being questioned as recent measurements of the viscosity of molten GaAs over a range of temperatures show an order-of-magnitude increase as the melting point is approached [43]. Because of the paucity of the heat-transfer properties, most of the values have been assumed in the mathematical models. As a consequence, the results of all the modelling work have been qualitative. Useful quantitative predictions have been limited by the availability and the accuracy of the various properties essential in modelling. Therefore, characterization of these properties is a prerequisite for accurate understanding of the processing conditions for the melt crystal growth system [73]. 2.3 Effects of Magnetic Fields on the Czochralski Growth Process There is no information found in the literature on the effect of a magnetic field on the effective thermal conductivity of a CZ melt. However, there have been numerous studies concerning the effects of magnetic fields on the CZ process, especially on melt flow. A literature review of these studies is given here as these studies are indirectly related to the influence of a magnetic field on the effective thermal conductivity of a melt within the CZ configuration. A well established effect of an applied magnetic field on the CZ process is the sup-pression of the melt flow. The investigation of the magnetic effects on the melt flow has involved both experimental work and numerical modelling. The forces driving fluid flow in a CZ melt are buoyant forces due to density differences, thermocapillary forces along the free surface due to surface tension gradient, centrifugal forces due to crystal rotation, and rotational forces due to crucible rotation [71]. The 15 melt convection affects the radial inhomogeneity of the crystal due to the radial flow field, the shape of the interface due to the isotherms or iso-concentration lines which are flow determined, the interface transition due to the flow transition near the interface, and the irregular striations caused by the temperature-flow fluctuations [71]. The thermal convection and the associated temperature fluctuations in the melt can be suppressed by the application of a magnetic field because the electrical conductivi-ties of molten semiconductors are comparable to those of liquid metals [65]. The flow suppression is due to the electromagnetic Lorentz force which opposes the motion of the electrically conducting melt across the magnetic field lines. The magnetic field can be applied either vertically or horizontally. The vertical magnetic field, generated normally by a coaxial solenoid, maintains the axisymmetry of the system and thus it is more widely used. The flow suppression by a vertical magnetic field has been observed during the CZ growth of silicon [42] and during the L E C growth of GaAs [40,76]. A horizontal magnetic field of 0.125 T has been found to suppress the temperature fluctuations at 4 mm below the melt-encapsulant interface from 18 to 0.1 °C for a 2-inch diameter GaAs crystal [45]. A horizontal magnetic field of 0.32 T eliminated the temperatue fluctuations in the melt for a 3-inch diameter GaAs crystal growth [46]. Analytical and numerical studies on the effects of a vertical magnetic field during the growth of the crystals have also shown the flow suppression in the melt [49,82,65]. Mihelcic et al. studied the effects of an axial magnetic field on the flow and temperature oscillations in the melt, and found that the axial magnetic field damped the flow velocity vectors and ehminated the temperature oscillations [51]. Numerical studies of the flow and heat transfer in the L E C growth of GaAs with an axial magnetic field have also been carried out by Sabhapathy and Salcudean [65]. It was found that thermal convection, especially near the crucible side wall, was significant even in the presence of a strong 16 magnetic field of 0.75 T. The quantitative difference between the influence of a stationary transverse field and a vertical magnetic field on the flow and temperature distribution which were forced to be asymmetrical by applying an asymmetrical temperature boundary condition at the crucible wall in a Si melt was calculated [52]. The vertical magnetic field was found to be more effective in reducing the asymmetry in the system than the stationary transverse magnetic field, particularly in the vicinity of the growth interface. The effect of a horizontal magnetic field on residual impurity concentrations in GaAs crystals has been studied by Terashima et al. [46]. They observed that the carbon and boron concentrations were decreased as the magnetic field strength was increased. As a result of the beneficial effects of the magnetic field, better quality crystals have been successfully obtained. A n axial magnetic field stronger than 0.1 T eliminated the irregular striations in a GaAs crystal having a diameter up to 3 inches [40]. With an applied magnetic field of about 0.3 T, a completely striation-free In-doped 3-inch diameter crystal has been grown without the need to optimize the crystal rotational rate [76]. Dislocation-free and striation-free GaAs crystals of 50-mm diameter have been reported [33]. 2.4 Pyrolytic Boron Nitride as a Crucible Material for the Processing of Semiconductors 2.4.1 Production, Crystal Structure, and Properties of Pyrolytic Boron Ni-tride Boron nitride (BN) is manufactured by sintering hexagonal B N powder with added sin-tering agents such as boric oxide, borate and high silicate glass [79]. Undesirable pores 17 and impurities are found in the sintered product. Pyrolytic boron nitride (PBN) is pro-duced by reacting N H 3 and BCI3 gases at temperatures in the range of 1800-1900 °C and pressures below 1 Torr to form B N and HC1. The B N is deposited on a graphite mandrel of desired shape [10]. The deposition process is shown schematically in Figure 2.3. When the vapor deposition is complete, the finished P B N is slipped off the mandrel. Articles of P B N are displayed in Figure 2.4 [55]. The deposition process enables the formation of high purity and high density products. As a result of the deposition process, P B N products consist of laminar sheets. Boron nitride has a hexagonal structure as shown in Figure 2.5. During deposition, B N basal planes are oriented parallel to the deposition surface [55]. As a result, P B N normal to the deposition surface is predominantly parallel to the hexagonal c-axis while P B N parallel to the deposition surface is mostly parallel to the a-axis [10,25]. The crystal structure of the deposits has been found to be affected by the total pressure and the deposition temperature of the process [79]. Precipitates having five-fold symmetry are present in P B N and have been studied using transmission electron microscopy [47]. The results show that they are multiply-twinned hexagonal B N with (112) twin planes, and the five-fold symmetry axis is [201]. P B N is a highly anisotropic material. The properties in the "a" direction are markedly different from those in the "c" direction, as can be seen from Table 2.1 [10,2]. 2.4.2 Applications of Pyrolytic Boron Nitride Molten semiconductor melts react with crucible materials. This is due to the high tem-perature and the high dissolving power of the melts. Thus, a suitable crucible material must be one that reacts sufficiently slowly with the melt and does not contribute detri-mental elements to the crystal with respect to its electronic properties [13]. The crucible material affects the quality as well as the properties of the crystal. The Figure 2.3: Schematic diagram of the P B N deposition process [10] Figure 2.4: Articles of P B N [55] 20 • Boron o Nitrogen L 4 4 6 A Figure 2.5: Hexagonal crystal structure of boron nitride [25] Property Units a Direction c Direction (a) Thermal expansion (RT to 1500°C) Thermal conductivity at RT at 800°C Tensile strength at RT Tensile strength at 2200°C Flexural strength a tRT Young's modulus at RT W / m - K W / n v K MPa psi MPa psi MPa psi GPa 106 psi 0.20-0.65 63 63 41 6000 103 15000 83 12000 22 5.0 1.5 2.9 2.8 410 No data No data No data Summary of Thermal Properties Summary of Mechanical Properties Property Thermal Conductivity "a" Direction i 25°C 500°C 1000°C "c" Direction (&> 25°C 500°C 1000°C Thermal Expansion "a" Direction @ 500°C 1000°C "c" Direction & 500°C 1000°C Coefficient of Thermal Expansion "a" Direction above500°C "c" Direction Resistance to Thermal Shock: 1200°C into Liquid Nitrogen Specific Heat @ 25°C 500°C 1000°C Unit caU cm-sec 0 C mm/mm mm/mm-0C cal/gm-°C Value .25 .17 .15 .004 .005 .006 .001 .0025 .013 .027 3 x 1 0 ' 30 x 10* no damage .2 .4 .47 Property Apparent Density Gas Permeability (Helium) Compression Strength "a" Oirection @ 2S°C 1200"C "c" Direction ® 25°C 1200°C Tensile Strength "a" Direction @ 25°C Flexural Strength @ 25°C 1200°C Torsional Shear Strength @ 25°C Young's Modulus "a" Direction @ 25°C Poisson's Ratio "a" Direction @ 2S°C Flexural Modulus > 25°C 1200°C Hardness—Taken on surfaced "a" plane Unit g/cc cm'/sec Psi Psi Psi Psi Psi Psi Knoop Hardness # Table 2.1: Properties of P B N : (a) from reference [10]; (b) from reference [2] 22 effects of crucible material on GaAs crystals have been investigated [78]. The results showed that the amount of arsenic lost by vaporization during growth was dependent on the crucible material. The electrical properties of the crystal such as n- or p-type conductivity, resistivity, carrier concentration, and Hall mobility were different using different crucible materials. The impurity and carbon concentrations were also affected. There are several materials that can be used for the processing of semiconductors. The most commonly used crucible material is fused silica (quartz). It provides the best compromise with respect to chemical stability, cost, ease of fabrication and mechanical strength [13]. Pyrolyric boron nitride (PBN) is increasingly being utilized particularly to grow semi-insulating (SI) GaAs because SI GaAs can be obtained using this crucible material without intentional doping of impurities [75]. The use of P B N crucibles instead of conventional fused silica crucibles minimizes melt contamination by Si through the dissolution of the quartz crucible, eliminating the need for compensatory doping with chromium to achieve SI quality. Another advantage is its lower cost per GaAs crystal, because each P B N crucible can be used many times, whereas a quartz crucible can be used only once [10]. Other new crucible materials which are also actively being investigated are aluminium nitride (A1N), silicon nitride (Si3N4) , silicon carbide (SiC) and hexagonal boron nitride (hBN) [78]. These materials have high melting points and hardly react at all with other materials at high temperatures in either oxidizing or reducing atmospheres. P B N is very attractive as a container material for elemental purification, compound-ing, and growth of compound semiconductor crystals because of its structure, properties (thermal and mechanical), purity, and chemical inertness [10]. Not only has it been used as a crucible material for the L E C process, but it has also been used as crucibles for the vertical gradient freeze (VGF) growth of crystals, evaporation crucibles for deposi-tion of metals and dopants at high temperatures and ultra-high vacuum by molecular beam epitaxy (MBE) . A P B N - P G resistance heater has been developed which enables 23 the growth of compound semiconductor crystal layers by rapid thermal processing in metal organic chemical vapor deposition ( M O C V D ) and other epitaxial processes. It has also been used as high-temperature jigs and insulators [47]. However, P B N material in these applications still needs to be improved. In order to optimize P B N performance, an accurate knowledge of the structure, thermal and mechanical properties of the material is required [10]. 2.5 Experimental Determination of Thermal Conductivity The general principles of thermal conductivity measurements will be reviewed based on reference [88]. The various standard methods of measurement will be briefly described. Heat is conducted through solids by various carriers: electrons, lattice waves or phonons, magnetic excitations, and electromagnetic radiation. The total thermal con-ductivity is the sum of contributions from each type of carrier which is given by k = IE Ct vi u (2.i) where i denotes the type of carrier, C, is the contribution of each carrier to the specific heat per unit volume, i>, is the velocity of the carrier, and /,• is an appropriately defined mean free path. The method for measuring the thermal conductivity of a material depends markedly on whether the material is a solid, liquid or gas. It also depends on the temperature range required and whether the material has a high or low thermal conductivity. No one method is suitable for all the required conditions of measurement. The appropriateness of a method is further detected by such considerations as the physical nature of the material, the geometry of samples available, the required accuracy of results, and the time and funds entailed as well as the availability of experimental apparatus. 24 Two general categories of thermal conductivity measurements are used, namely steady-state measurements and transient methods. In the steady-state methods, the test spec-imen is subjected to a temperature profile which is time-invariant, and the thermal conductivity is determined directly by measuring the rate of heat flow per unit area and the temperature gradient after equilibrium has been reached. In the transient methods, the temperature field in the specimen varies with time, and measurements of the rate of temperature change replace the measurement of the rate of heat flow. The transient measurement normally determines the thermal diffusivity instead of the thermal con-ductivity. The thermal conductivity is calculated from the thermal diffusivity using the density and specific heat capacity of the test material. Obtaining a controlled heat flow in a prescribed direction such that the actual bound-ary conditions in the experiment agree with those assumed in the theory is the primary concern in most methods of measurement. In theory, the simplest method to accomplish this task is to use a specimen in the form of a hollow sphere with a heater in the cen-ter. The heat supplied by the heater passes through the specimen in a radial direction without loss. However, in practice it is very difficult to fabricate a spherical heater which produces uniform heat flux in all radial directions. Morover, it is difficult to fabricate spherical specimens and to measure the heat input and the temperature gradient in such an experimental arrangement. A method more commonly used to control the heat flow in the prescribed direction is the use of quard heaters combined with thermal insulation in most cases so adjusted that the temperature gradient is zero in all directions except in the direction of desired heat flow. In most methods of measuring thermal conductivity, a cylindrical specimen geometry ranging from a long rod to short disk is utilized, and the heat flow is controlled to be in either the longitudinal (axial) or the radial direction. 25 Detailed description, of each of the standard methods for obtaining the thermal con-ductivity of a material would take pages to complete, and thus only the main ideas are mentioned here. The references listed in reference [88] should be consulted for detailed discussions of the methods. 2.5.1 Steady-State Methods Steady-state procedures for measuring the thermal conductivity of solids include the following methods: (a) longitudinal heat flow, (b) Forbes' bar, (c) radial heat flow, (d) direct electrical heating, (e) thermoelectrical, and (f) thermal comparator. Depending on the means of measuring the heat flow, the longitudinal and radial heat flow methods are divided into absolute and comparative methods. In the absolute method, the rate of heat flow into a specimen is directly determined, by measuring the electrical power input to a heater at one end of the specimen, and the rate of heat flow out of a specimen is measured with a flow or boil-off calorimeter. In the comparative method the rate of heat flow is calculated from the temperature gradient in a reference sample of known thermal conductivity which is positioned in series with the specimen. It is assumed that the same heat flow occurs in both the reference and test samples. In the longitudinal heat flow method, the experimental arrangement is designed such that the flow of heat is only in the axial direction of a rod or disk specimen. The radial heat loss or gain of the specimen is prevented or minimized and evaluated. Assuming no radial heat loss or gain, under steady-state conditions the thermal conductivity is deter-mined by the following expression from the one-dimensional Fourier-Biot heat-conduction equation: where k is the average thermal conductivity corresponding to the average temperature 26 | (Ti + T 2 ) , A T is the temperature difference between T\ and T 2 , g is the rate of heat flow, A is the cross-sectional area, and Ax is the distance between points of temperature measurements for T\ and T 2 . The same idea is used in the radial heat flow method in which the axial heat flow is minimized. Other steady-state methods are variants of the longitudinal and radial heat flow methods with minor differences in the equations and experimental arrangements employed to obtain the thermal conductivity. 2.5.2 Nonsteady-State Methods The category of nonsteady-state methods includes the periodic and the transient heat flow methods. Each of them is also divided into longitudinal and radial heat flow methods. The transient heat flow methods incorporate the flash method, the line heat source and probe methods, the moving heat source method, and several comparative methods. In this method, the temperature field in the specimen varies with time. The rate of temperature change at- certain positions along the specimen is measured and no mea-surement of the rate of heat flow is required. These methods usually yield the thermal diffusivity from which the thermal conductivity can be calculated by knowing the density and specific heat capacity of the test material. In the periodic heat flow methods, the heat supplied to the specimen is modulated to have a fixed period, and the resulting temperature wave propagating through the speci-men with the same period is attenuated as it moves along. The thermal diffusivity can be determined from measurements of the amplitude decrement and/or phase difference of the temperature waves between certain positions in the specimen. The transient heat flow methods involve heating the specimen to a desired temper-ature and then letting it cool. The temperatures at two positions of the specimen are measured as a function of time. The thermal diffusivity can then be obtained from these measurements. Many varieties of the transient heat flow methods have been developed 27 employing similar ideas with various different heat sources and experimental arrange-ments. 2.5.3 Methods for Plate Specimens Plate methods are suitable for poor thermal conductors and insulators. This specimen geometry minimizes the ratio of the lateral heat losses to the heat flow through the specimen [88]. Plate methods can be categorized into absolute and comparative methods [88]. Two main kinds of experimental arrangements are used for the absolute method: the single-plate system and the twin-plate system. The single-plate system requires only one speci-men which is placed between a hot plate and a cold plate. The twin-plate system requires two similar specimens to be sandwiched between a hot plate in the middle and two cold plates on the outside. In the comparative method, a reference sample of known thermal conductivity is placed in series with the unknown specimen in which it is assumed the same rate of heat flow through both the reference and the specimen. The specimen and the reference sample are sandwiched between a hot plate and a cold plate. Comparative methods have the advantages of simpler apparatus, easier specimen fab-rication, and easier operation. Their disadvantages include the additional measurement errors due to the required additional measurements of temperatures and thermocouple separations, the difficulty in matched guarding, and the lower accuracy due to the addi-tional uncertainty in the thermal conductivity of the reference sample, the conductivity mismatch between the specimen and reference sample, and the interfacial thermal contact resistance [88]. Numerous techniques have been developed to obtain either the thermal conductivity (steady-state) or the thermal diffusivity (transient) of a plate specimen. Steady-state measurements include the de Senarmont method and the guarded hot plate method. 28 Transient measurements cover the calorimetric method, the flash method, and several other methods. In the de Senarmont method, the sample is coated with a thin film of wax. Heat is applied at a central point by means of a hot, thin silver tube tightly fitted in a hole at the center of the plate. The wax melts around the region where the heat is supplied and the bounding line of the melted wax is the visible isotherm. This method is used to determine anisotropy in thermal conductivity. However, it does not yield absolute values of thermal conductivity, and axial heat loss is not prevented. The guarded hot plate method has been adopted by the American Society for Testing and Materials (ASTM) as a standard method for thermal conductivity measurements. This method is described in the A S T M handbook under the fixed designation C177 [1]. However, it requires a sophisticated experimental apparatus. In the calorimetric method, a part of a thin sample is shadowed by a mask from chopped light irradiation and then periodic thermal energy is imparted to the remaining part of the sample [35]. The temperature at a position of the sample lying under the mask is detected by a fine thermocouple attached to the sample. To obtain the thermal diffusivity parallel to the plane surface of the sample, the temperature T is measured as a function of the distance x between the positions of the thermocouple and the edge of the mask. The thermal diffusivity can be calculated from the following equations: k = (u>/2a)2 where Q is the amplitude of the applied thermal-energy flux per unit area, u> is the angular frequency of the periodic heating, C is the heat capacity of the sample per unit volume, d is the thickness of the sample, k is the thermal decay constant, and a is the thermal diffusivity of the sample. 29 In the flash method, thermal energy from a flash tube or laser is imparted instan-taneously to the front surface of the plate-like sample by a light pulse [35]. The time interval of the flash is short compared to the time required for the resulting transient flow of heat to propagate through the specimen. The temperature of the rear surface is measured as a function of time. Assuming that the heat pulse is instantaneous and uniform over the entire front surface, the surface heat losses are negligible, and the heat flow is one-dimensional, the thermal diffusivity can be calculated from the relation 0.139 a? a = *l/2 where ti/2 is the time corresponding to a rise in the temperature to one-half of its maxi-mum value. There are several other methods more or less related to the transient methods out-lined above. These include the mirage-effect method, the picosecond time-resolved ther-moreflectance method, and the photoacoustic method. The procedures are described in reference [35]. 2.6 Reported Thermal-Conductivity Values for Pyrolytic Boron Nitride The thermal-conductivity values reported for P B N in the "a" and "c" directions which have been obtained from two different references are shown in Table 2.2 [2] and Table 2.3 [10]. Other sets of values for slightly different P B N materials are presented in Figure 2.6 and Figure 2.7 where Ka denotes the thermal conductivity in the basal plane (a direction) [25]. In Figure 2.6, specimen C has a higher Ka because it has a higher density than specimen B . Comparing the numbers in the two tables, it is evident that they are not in agreement with each other. The method used to obtain these values was not given. Thus, the values are suspicious in terms of their reliability and accuracy. The values shown in the figures 30 1 go o « ° * < > o o o 0 1 o 1 o - o — * o * • c Oo ° ° o ° o A A A B A A A A * * A — A * / 1 1 1 100 200 T(K) 300 400 Figure 2.6: Thermal conductivity of P B N reported in reference [25] for specimen B and C 31 10.0 .0 I £ o r -r -< O 0.1 0.01 0.001 1 1 1 I 1 1 111 1 1 1 1 1 1 1 1 1 1 I I -_ <pOXC«DC» -— o <£r> c — ocP o - o -- o -o o - CO — o — o - o — z o - <? 8 -6 8 o o cP o Mill -( 1 1 1 1 1 1 1 1 1 1 1 1 (Mil 1 I I I 10 100 T ( K ) 5 0 0 Figure 2.7: Thermal conductivity of P B N reported in reference [25] for specimen C 32 Temperature "a" direction "c" direction Unit 25 °C 104.5 1.672 Wjm.°C 500 °C 71.6 2.09 W/m.°C 1000 °C 62.7 2.51 W/m.°C Table 2.2: Thermal conductivity of P B N reported in reference [2] Temperature a direction c direction Unit room temperature 63 1.5 W/m.°K 800 °C 63 2.9 W/m.°K Table 2.3: Thermal conductivity of P B N reported in reference [10] are also different from the ones shown in the tables, but this is due to the different materials used. The values obtained were rather high, as can be seen from the figures. For instance, at room temperature, specimen C has a value of 390 W/m.°C which is comparable to that of copper (about 400 W/m.°C). Such high values are surprising for this type of material which is regarded as a thermal insulator (ceramic). Clearly, discrepancies in the reported thermal conductivity values exist for P B N ma-terials that are supposedly the same. Moreover, these values are not representative of the thermal conductivity of a P B N crucible. In a P B N crucible, the direction along the wall is parallel to the deposition surface which consists of predominantly the "a" direction, and the direction through thickness is perpendicular to the deposition surface comprising of predominantly the "c" direction. Chapter 3 Finite Element Solutions for the Two-Dimensional Heat-Conduction Equation The differential equation of heat conduction for a stationary, homogeneous, isotropic solid with heat generation or consumption within the body has the form [54,11,27,31] V.[kVT] + Q = pCp^ (3.3) where V is the differential operator, k is the thermal conductivity, T denotes tempera-ture, Q is the rate of heat generation or consumption, p is the density, Cp is the specific heat capacity, and t denotes time. In two-dimensional cylindrical coordinates, i.e. ax-isymmetrical system, „ d , d , V = — ir + — iz or oz where r and z are the radial and axial coordinates respectively, and ir and i2 are the corresponding unit vectors. For two-dimensional cartesian coordinate system, „ d , d , ox oy with x and y as the coordinate variables, and i and j the corresponding unit vectors. When the thermal conductivity is assumed to be constant, i.e. independent of position and temperature, and if there is no heat generation and consumption in the solid, the two-dimensional heat-conduction equation reduces to ld_( dT\ d2T _1 dT r dr 1 dr I dz2 a dt ' 33 34 in the cylindrical coordinates, and to d2T d2T 1_ dT_ a ~dt (3.5) dx2 dy2 in the cartesian coordinates, where a is the thermal diffusivity defined as In the special case of the Czochralski crystal growth, a quasi-steady-state assumption is usually used to calculate the temperature distribution in the crystal. The term |jp in equation 3.4 is transformed into ^ + v [9] to take into account the transient and convective heat transfer of the crystal respectively due to its upward motion, where v is the velocity of the moving crystal. However, according to Rosenthal's theory for moving heat sources, an observer moving with the coordinate system fails to detect any change in temperature with time in his surroundings [9]; hence, ^ is zero leaving only the convective term v in the transformed equation Equation 3.7 is expected to be valid after a crystal of some length has grown [9]. Equation 3.7 is no longer transient, and neither is it steady-state since the right-hand side still contains a non-zero term. Accordingly, the transformed equation is said to be quasi-steady-state. Solution to the heat conduction equation is often impossible analytically; only with very simple boundary conditions is the analytical solution obtainable. Consequently, nu-merical techniques are required for solving the temperature field in most cases. The most common numerical techniques employed are finite difference method and finite element method. The finite element method is used to obtain the solution as it can easily handle an irregular domain geometry and hence more complex boundary conditions. Moreover, dz (3.7) 35 the finite element method generates numerical approximations that satisfy certain global conservation laws, independent of the boundary conditions and geometry, which lead to computationally well-behaved results [23]. It also gives solutions for the entire domain instead of only at the nodal points. Generally, it yields a better approximation to the exact solution compared to the finite difference solution. 3.1 Finite Element Equations Finite element equations can be obtained by means of weighted residual methods. The most popular weighted residual method is the Galerkin method. In the finite element method, the domain of interest is discretized into a finite number of elements. The derivation of the equations is performed on each individual element. The elemental matrices and vectors are then assembled by adding their contributions at the elemental nodes to the corresponding global node locations in the domain to form the domain stiffness matrix and load (force) vector. The finite element formulation of the unsteady-state heat-conduction equation gen-erally yields a set of simultaneous algebraic equations in the form [K] {T} + [C] {T} = {/} (3.8) where [K] is the thermal stiffness matrix, [C] is the capacitance matrix, {/} is the thermal load (force) vector, {T} is the nodal temperature vector, and {T} is the vector containing time derivatives of nodal temperatures. The thermal stiffness matrix can be decomposed into [K] = [Kc] + [Kh] + [Kr] where [Kc] is the conduction matrix, [Kn] is the convection matrix, and [Kr] is the linearized radiation matrix. Similarly, the thermal load (force) vector can be subdivided 36 into {/} = {//»} + {fr} where {//,} is the convective vector and {/ r} is the linearized radiative vector. For the quasi-steady-state case (equation 3.7), the matrix equation is slightly different from equation 3.8, and is given by l[K} + [Cv}}{T}={f} (3.9) where [Cv] is the convective capacitance matrix. Equation 3.9 does not involve the time derivative of temperature but instead has the convective term due to crystal motion added to [K]. The finite element equations have been derived for both the cylindrical and the carte-sian coordinates using the Galerkin method. Full derivations of the equations are given in Appendix A . Good discussions on the finite element method can be found in references [57], [34], [68], [58] and [28]. After going through the derivation, the elemental matrices and vectors for the cylin-drical system are K Q = J r r h c N i N j d T JjfW = J r r hr Ni Nj dT C\f = JzJrrpCpNiNjdrdz ftf =JrrhcNiT00ldT ft? =JrrhrNiT002 dF and C $ =fsfrrvPCpNiNtjdrdz where N is the shape or trial function, hc is the convective heat-transfer coefficient, hr is the radiative heat-transfer coefficient, T 0 0 l and are the ambient temperatures, T 37 denotes domain boundary, and Nz is the derivative of N with respect to z. The subscripts i and j denote the number of rows and columns in the matrix, respectively. The equations in the cartesian coordinates are slightly different in that the r terms in the above equations are removed. Presented in a similar manner, the elemental matrices and vectors for the cartesian system are K $ = Jr hc ^ Nj dT r " r 3 (3.11) C\? = S J x P C p N i N j d x d y fif — Sr he Ni Too, dT = J r h r N i T 0 0 2 aT The next essential step in the finite element formulation is to choose the type of element and its associated shape functions. The type of element and shape functions chosen have an important implication on the accuracy of the solution. A very popular element type is the quadratic isoparametric element which has been used quite extensively and successfully. The quadratic isoparametric element considered here consists of four sides and eight nodes. The sides of the element can be curved due to the quadratic nature of the shape functions. The quadratic isoparametric element is shown in Figure 3.8a. The 38 Figure 3.8: Element type: (a) 8-noded quadratic isoparametric element; (b)Quadratic isoparametric boundary element 39 corresponding shape functions are JVi = - ( l - 0 ( l - i 7 ) ( l + f + i7)/4 N2 = -{l + Q(l-r,)(l-t + r,)/4 JVrs = - ( l + 0 ( l + ' 7 ) ( l - f - » / ) / 4 N4 = -{l-Q(l+r,)(l+t-V)/i AT5 = ( l - ^ ) ( l - r ? ) / 2 iV 6 = ( l - r ? 2 ) ( l + 0 / 2 Ar7 = ( l - ^ ) ( i + T ? ) / 2 JV8 = ( l - i f ) ( l - 0 / 2 where —1 < £, 77 < 1. At domain boundary, one-dimensional quadratic shape functions are used to account for convective/radiative boundary conditions. The boundary element is shown in Figure 3.8b and the shape functions are given by iV 1 = - 7 , ( 1 - ^ / 2 N2 = 77(1 + 17)/2 V " (3.13) N3 = 1 - 772 Since the shape functions are a function of £ and 77, the r-z or x-y space (global coor-dinates) must be mapped into the £-77 space (local natural coordinates). The coordinate transformation is shown in Figure 3.9, where r = Er,-JVi(e,i ?) z = J 2 z i N i ( t i V ) 1=1 i=i The mapping functions used to accomplish the transformation are the same as the shape functions as the element is isoparametric. To evaluate the integrals of the elemental matrices and vectors, numerical integra-tion must be used since the integrands are not simple functions that permit closed-form 40 Figure 3.9: Mapping of the quadratic isoparametric element into the (£, T/) coordinates: (a) domain element; (b) boundary element 41 integration [34]. A numerical integration technique most often used is the Legendre-Gauss method (Gaussian quadrature). The method uses a number of sampling points to evaluate the integrals. A larger number of sampling points results in better accuracy but is more computationally intensive. Using nine sampling points in the domain and three sampling points along the boundary proves to be sufficient. The equations in the cylindrical coordinates after the numerical integration are mm rp m rp [/#>] = £ WS rM hc {NM} {NM} A M «=1 [/#>] =^Wi r{rK) hr {NM} {NM} A M i=l mm rp |c(e)] = EE^^^,,)^p{%,,)}{%,,)} M-W] <3-14) i=l 3=1 m {fie)} = E Wi ' ( , ) K {N{v)} Too, A M m = E W> r M hr {NM} TOO2 A(,..) mm rp [^]=EE WiWmwvpc, {N^} {N,^)} dc*[j(ftfI8)] t = l 3=1 and in the cartesian coordinates m m rp t = l J = l m i = l . m ™ [ i t f >] = £ W- fc, { iV ( r K ) } {JV ( | B ) } A M (3.15) m m „ [ C ( e ) ] = E E Wi W ^ C P { % , „ ) } { % , , ) } <fef i = l J = l 42 where m {/„e)}=£ {*(„)} AM i=l t n {/r(e)} = E ^ ^ K*)} A(*) t = l m = number of sampling points in one dimension Wi and Wj = weighting constants £ and 77 = coordinates at sampling points The elemental matrices and vectors are calculated using computer programs, and placed in the corresponding locations in the global matrix and vector. The global matrix equations are solved using a suitable matrix solver to obtain the nodal temperatures. The equation to be solved in steady-state heat conduction [i<] m = {/} is algebraic in nature, and hence can be evaluated directly using matrix algebra. However, the transient heat-transfer equation [K] {T} + [C] {T} = {/} involves the time derivatives of temperature and is generally difficult to solve analytically. To overcome this problem, a time-stepping recurrence technique is employed in which temperature is approximated by {T} « £ i V , {Ti} (3.16) where {7}} is the temperature vector at time 7. The shape functions, iV}, are a function of time only. In this method, trial functions are used to discretize the time domain and 43 the solution at each time step is obtained using the temperatures from the previous time step(s). The number of trial functions used influences the accuracy and stability of the solution; in general the greater is the number of shape (trial) functions used the better will the accuracy and stability be [12]. The recurrence method which uses two shape functions is known as a two-point (lin-ear) recurrence scheme, while a three-point (quadratic) scheme uses three shape functions. The dimensionless time variable £ and shape functions for the two-point scheme are o < C < i C = £ Ni = l - C; Ni = -k (3-17) J V i + 1 = C; JVi+i = i and for the three-point scheme - i < c < i c = ih j V ^ - i c U - C ) ; A i - i = (-* + c) & JVi = ( i -C) ( i + 0 ; Ni = - ¥ t Ai+i = i C ( i + 0 ; = (! + <) & (3.18) where the subscripts I — 1, I and / + 1 denote previous time step, current time step and next time step, respectively. Two most widely used techniques are chosen for the two time-stepping schemes. The Crank-Nicholson method is employed in the two-point recurrence scheme, and the result-ing equation is given as [K] { ^ i ± L ± ^ } + [C] { ? » ! ^ - } « {/} (3.19) The Dupont three-level technique is used in the three-point scheme resulting in [ K ] 13r,+,H-r,_,j + [ c ] | ^ { / } ( 3 2 o ) 44 Both schemes are utilized in the following fashion. The two-point scheme which only requires the temperatures from previous time step is used to obtain the solution at the very first time step. The three-point scheme which requires the temperatures from two previous time steps is then used for all subsequent calculations. 3.2 Comparisons with Analytical Solutions In order to validate the finite element equations, computed results must be compared with analytical solutions. Checking the finite element solutions against the analytical solutions under different boundary conditions is necessary to ensure the validity of the mathematical models developed in chapter 4 and chapter 5. In chapter 4, a steady-state model is utilized to examine the effect of a magnetic field on heat transfer in a melt. In chapter 5, a transient model is required to determine the thermal conductivity of a P B N crucible plate. The models involve fixed, insulated and convective boundary conditions. Thus, in the following sections, thorough comparisons are made in one- as well as two-dimensional conduction, convection and transient solutions. The derivations of the analytical solutions are included in Appendix B . 3.2.1 Steady-State Conduction For pure steady-state conduction, the equation to be solved is [Kc] {T} = 0 (3.21) The boundary condition can be either fixed temperature or zero heat flux at the boundary. Before comparing with analytical solutions, the finite element calculations must satisfy two criteria. The first is that the computed temperatures must be zero at all nodes when no constraint is applied at the boundary. The second is that when a fixed temperature 45 case 1 case 2 case 3 case 4 r,- (m) r0 (m) 0.01 1.61 1000.01 1001.61 0.0001 1.6001 0.0001 3.2001 Table 3.4: Different pairs of values of r,- and ra used for comparisons of steady-state conduction finite element results with 1-dimensional analytical solutions is applied along the boundary the solution temperatures must be equal to the boundary temperature. 3.2.1.1 One-Dimensional Solutions The analytical solution is for a hollow tube with fixed temperatures at the inner and outer surface. The solution is [27] l n ( ^ ) T(r) =Ti + (T0 - Ti) — W - (3.22) ln0?) where T,- and Ta are the inner and outer surface temperatures, and r, and r0 are the inner and outer radii respectively. T is a function of r only. Since this is a steady-state con-duction problem with constant thermal conductivity and fixed boundary temperatures, the thermal-conductivity value does not affect the temperature profile. Ti is set to 0 °C and T0 is fixed at 1000 °C for the calculations. Four cases are examined for different combinations of r,- and r0 values as shown in Table 3.4. In all cases, 5 different numbers of elements of equal size are used to illustrate the convergence of the finite element results to the analytical solutions. The number of elements is doubled from 2 to 4, 8, 16 and 32. Figure 3.10 shows the results of the comparison for the first case. The finite element results clearly become closer to the analytical solutions as the element number increases. In this case, the solutions using 32 elements give essentially the same values as the 46 Figure 3.10: Comparison of finite element solutions with 1-dimensional analytical solu-tions for steady-state conduction using r,- = 0.01 m and r0 = 1.61 m 47 analytical. The error is larger at the region near the z axis where the temperature gradient changes rapidly. This is due to the inability of the quadratic element shape functions to fit a high curvature. Changing the shape functions to a higher-order polynomial or increasing the number of elements locally will improve the fit. In the second case, the inner radius is very large compared to the one in case 1 while the thickness of the tube remains unchanged. The temperature profile in such a tube is virtually linear, and therefore is easily fitted by the finite element solutions. Figure 3.11 illustrates the excellent agreement between the solutions. Even with only 2 elements, the finite element solutions are indistinguishable from the analytical solutions. The comparison for case 3 is displayed in Figure 3.12. As expected, the solutions converge as the element number is increased. However, the finite element solutions are unable to fit the region near the z axis well even with 32 elements. This case differs from case 1 in that the inner radius for this case is much smaller. As r*j gets smaller, the change in temperature gradient becomes more severe in the region close to the inner surface of the tube since the inner surface is nearer to the z axis. For this reason, the fit in this case is not as good as the fit in case 1. Increasing the thickness of the tube results in a larger change in the temperature gradient near the z axis. Figure 3.13 presents the results of comparison for case 4. This case uses the same r,- as case 3 but the thickness of the tube is twice the thickness in case 3. By comparing Figure 3.12 and Figure 3.13, it can be seen that the temperature profile in Figure 3.13 has a slightly higher change in the temperature gradient resulting in a larger error in the finite element solutions. 3.2.1.2 Two-Dimensional Solutions The finite element solutions are also compared with two-dimensional analytical solutions to further substantiate the convergence of the finite element solutions. The analytical 48 Figure 3.11: Comparison of finite element solutions with 1-dimensional analytical solu-tions for steady-state conduction using rt- = 1000.01 ra and r„ = 1001.61 ra 49 0.2 0.4 0.6 0.8 1 Radial Distance (m) 1.2 1.4 1.6 Figure 3.12: Comparison of finite element solutions with 1-dimensional analytical solu-tions for steady-state conduction using rf- = 0.0001 m and r 0 = 1.6001 m 50 I & O a -0.4 0.8 1.2 1.6 2 Radial Distance (m) 2.4 2.8 3.2 Figure 3.13: Comparison of finite element solutions with. 1-dimensional analytical solu-tions for steady-state conduction using r,- = 0.0001 m and r0 = 3.2001 m 51 profile 1 1'mid % for 0 < z < \ c 1'mid C Tmid Z for \c < Z < C profile 2 Tmid sin (s .) for 0 < z < c Table 3.5: Two different temperature profiles along the axial boundary for comparisons of steady-state conduction finite element results with 2-dimensional analytical solutions solution is for an axisymmetrical domain having a radial length b and an axial length c. The finite element domain is divided into n x m elements of equal size; n refers to the number of elements radially while m axially. The element numbers used are 1 x 2, 2 x 4, 3 x 6, 4 x 8 and 5 x 10. b and c are set to 1 m and 2 ra, respectively. The top and bottom boundaries are fixed at 0 °C. Two temperature profiles are imposed on the side of the domain as listed in Table 3.5. Results are obtained by setting the mid-height temperature Tmid to 1000 °C in both cases. Using profile 1, the resulting analytical solution is [54] rpt s 2Tmid ^ I0(r]mr) • J 1 . (rimc\ c (r)mc\\ T(r,z) = 2^ — — s i n ( 7 7 m 2 ) < — sin I— - cos I — - ) \ 0 < z < l-c ~ ~ 2 ^ s i n { V m C ) + ^ c o s ^ ( 3 . 2 3 ) — c<z<c 2 52 where rrnr Vm = c The solution for profile 2 is (3.24) where Jo() is the modified Bessel function of order 0 of the first kind. Figure 3.14 shows the comparison of the solutions along the z axis for profile 1. The convergence is clear with increasing number of elements. Except for 1 x 2 elements, the other solutions give a good agreement with the analytical. The corresponding comparison along the midline between the z axis and the boundary is shown in Figure 3.15. The finite element .solutions along this line are closer to the analytical than those along the z axis. It is expected that the solutions become more accurate in the region nearer to the boundary since the temperature along the boundary is fixed. The results of comparison for the second profile are presented in Figure 3.16 and Figure 3.17 along the z axis and the midline, respectively. Using this profile, the finite element solutions are generally better than those using profile 1. This can be seen by comparing Figure 3.16 with Figure 3.14 and Figure 3.17 with Figure 3.15. This is due to the discontinuity of the temperature slope in profile 1 which causes sharper changes in temperature gradient axially. To see the convergence of the whole temperature distribution, two-dimensional con-tour plots for profile 2 are compared in Figure 3.18. 3 finite element solution fields using 2 X 4, 4 X 8 and 8 x 16 elements are shown in the figure. Wi th 8 x 16 elements, the solution field is virtually identical to the analytical temperature distribution. 53 Figure 3.14: Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state conduction using temperature profile 1 along the z axis 54 Figure 3.15: Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state conduction using temperature profile 1 along the midline between the z axis and boundary 55 Figure 3.16: Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state conduction using temperature profile 2 along the z axis 56 700 I — i — i — i — i — i — i — i — i — i — i — i — i i i i i i i r Axial Distance Figure 3.17: Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state conduction using temperature profile 2 along the midHne between the z axis and boundary Figure 3.18: Comparison of finite element temperature distributions with 2-dimensional analytical solutions for steady-state conduction using temperature profile 2: (a) 2 x 4 elements; (b) 4 x 8 elements; (c) 8 x 16 elements; (d) analytical 58 3.2.2 Steady-State Convection For steady-state convective boundary conditions, the equation to be solved is [ [Kc] + [Kh]] { r} = {/„} (3.25) Equation 3.25 describes the combined conductive and convective steady-state heat trans-fer. Before comparing with analytical solutions, the equation must be tested against the following criterion. Using a constant heat-transfer coefficient h and a uniform ambient temperature Too, the resulting finite element solutions must be equal to the ambient temperature. The finite element calculations are compared with the analytical solutions of a two-dimensional axasymmetrical domain of radial length b and axial length c. Along the top and bottom boundaries fixed temperatures are imposed, while convection occurs along the side boundary to an ambient temperature of 0 °C with a constant h. With the top boundary Tt0 set to 0 ° C and the bottom boundary Tbo to a constant positive temperature T^,, the analytical solution is given by [54] r(r ~\ - 2Tb° V H s i n M£" ( c ~ *)) MM h s i n M / M J 0(/5m6) { 6 ' 2 b } where and /3m's are the positive roots of An Wmb) =H J0(flmb) and J 0 () and Ji() are the Bessel functions of order 0 and 1 of the first kind, respectively. The numbers of elements used to obtain the solutions are 1 x 2 , 2 x 4 , 4 x 8 and 8 x 16 elements. The size of the domain, the boundary conditions and the heat-transfer properties used are listed in Table 3.6. 59 b 1.0 m c 2.0 m Tu, 0.0 °C Tbo 1238.0 ° C T 0.0 °C k 1.0 W/m.°C h. 100.0 i y / m 2 . ° c Table 3.6: Size of domain, boundary conditions and heat-transfer properties used for comparison of steady-state convection finite element results with 2-dimensional analytical solutions The comparisons are shown in Figure 3.19, Figure 3.20 and Figure 3.21 along the z axis, the midline and the boundary, respectively. Obvious convergence is observed. The finite element results are the least accurate along the boundary. Again, this is due to a large change in the temperature gradient along the boundary. The corresponding two-dimensional contour plots are displayed in Figure 3.22. The finite element solutions are obtained using 2 x 4, 4 x 8 and 8 x 16 elements to show the convergence of the temperature field. The results of 8 x 16 elements resemble the analytical solutions very well. 3.2.3 Transient Conduction In transient heat-conduction problems, the equation to be solved is [Kc] {T} + [C] {f} = 0 (3.27) There is no heat exchange between system and ambient surrounding. Comparisons are made for 1-dimensional solutions in the cartesian coordinates. Two problems are considered. The first problem has insulated ends and a linear initial temperature profile with 0 °C 60 Figure 3.19: Comparison of finite element solutions with 2-dimensional analytical solu-tions for steady-state convection along the z axis 61 1300 i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — r Axial Distance (m) Figure 3.20: Comparison of finite element solutions with. 2-dimensional analytical solu-tions for steady-state convection along the midline between the z axis and the boundary 62 i — i — i — i — i — i — i — i — i — i — i — l i i i i i i r + 1x2 elements -0 2 x 4 elements -A 4 x 8 elements _ X 8 x 16 elements analytical r/R = 1.0 i i i i i i I I I 1 I I L 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Axial Distance (m) Figure 3.21: Comparison of finite element solutions with. 2-dimensional analytical solu-tions for steady-state convection along the boundary 63 Figure 3.22: Comparison of finite element temperature distributions with 2-dimensional analytical solutions for steady-state convection: (a) 2 x 4 elements; (b) 4 x 8 elements; (c) 8 x 16 elements; (d) analytical 64 at x = 0 and Tmax at x = L, where L is the length of domain. The analytical solution is [54] oo mt j.\ Tmax 2 Tmax —^\ —aQ2 t la \ T{x,t) = + L 2 }^ e m -cos (An*) m—1 r - in t\i cos(P™L) 1 L sm(BmL) + (3.28) An f% L is equal to 1 m in the calculations. The temperature at x = L is 1000 °C (i.e. T m a a ; = 1000 °C). The thermal diffusivity, a, used is 10" 4 m2/-s. The calculations for 2 different time-step sizes, 300 s and 150 s, are compared with the analytical solutions. A total of 5 equal-size elements is used; increasing the number of elements does not change the results significantly. The transient responses are plotted for x = 0 and x = L in Figure 3.23. The smaller time-step size clearly gives a better approximation. The finite element results are not so accurate during the first 2 time steps. The solutions become more accurate with time. The larger error during the first time step is probably due to the two-point time-stepping scheme used for obtaining the solution which is not as accurate as the three-point scheme used for the subsequent calculations. The second problem is for a semi-infinite region initially at a constant temperature To. At t > 0, a fixed temperature of 0 °C is imposed at x = 0. The analytical solution for this problem is given by [54] T(x,t) = T0erf(^=^ (3.29) where er/() is the error function. To is 1000 °C in the calculations. The thermal diffusivity used is 1 0 - 4 m2/s. 15 equal-size elements are used for a bar of 0.3 m in length in the finite element calculations. Figure 3.24 presents the transient response at x = 0.01 m. Two time-step sizes are tested, 1 s and 0.5 s. As expected, more accurate results are obtained with the Figure 3.23: Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient conduction at x = 0 and x = 1 m 66 Figure 3.24: Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient conduction at x = 0.01 m 67 smaller time-step size. Again, large errors are observed during the first few time steps and the accuracy of the solutions improves with time. 3.2.4 Transient Convection The equation to be solved is [[Ke] + {Kh}} {T} + [C] {f} = {/„} (3.30) The above equation combines the conductive and convective transient heat transfer. The analytical solution employed for comparison is for an one-dimensional slab of length L which is insulated at x = 0 and experiencing convective heat transfer with a medium of 0 °C at x = L. Initially the slab is at a constant temperature To. The solution is given by [54] T(x t) - 2Tn V e-a^t H cos(/?ms) ( ' ' ° h L(fPm+ IP) + H cos(M ( 3 - 3 1 ) where An's are the positive roots of An tan(A»£) = H The solutions are obtained for H = 100 m - 1 and a = 10~4 m?/s. The length of the slab is 1 m and T 0 is 1000 °C. The transient response at the convective boundary (i.e. x = 1.0 m) is shown in Figure 3.25. 3 equal-size elements are used for the finite element calculations. The two time-step sizes are 2 s and 1 s. The smaller time-step size yields better results. In the region of rapid change of slope, the finite element solutions fail to accurately represent the analytical solutions. This is because of the quadratic time-stepping shape functions used which are unable to fit a large change in slope. A smaller time-step size will improve the fit in this region. 68 1000 900 800 700 + At = 2 s o At = 1 5 — analytical £ 600 1=5 o O 500 (X ° £ 400 300 200 100 0 x/L = 1.0 20 40 Time 60 80 Figure 3.25: Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient convection at x = 1.0 m 69 A comparison is also done in the cylindrical coordinates. The problem considers a cylinder with negligible internal resistance, i.e. the resistance to conduction is much less than the resistance to convection from the surface. The Biot number (Bi = ^ ) must be much less than 1.0 for this to be true. The solution to this problem with T = T(r) and ^ = 0 is [27] where T 0 is the initial temperature, A is the surface area and V is the volume of the cylinder. Solutions are obtained for a cylinder having a radius (R) of 0.001 m and an axial length of 6 m. The thermal diffusivity value used is 1 0 - 9 m2/s and the Biot number is 0.01 in the finite element calculations. The initial temperature is 1000 °C. The results are shown in Figure 3.26 for the thermal history along the convective boundary (r = 0.001 m). Only 3 elements of equal size are used for the finite element calculations. The two time-step sizes, 1 s and 0.5 s, give essentially the same solutions as the analytical. 3.2.5 Discussion on the Finite Element Solutions A comprehensive validation of the finite element equations against the analytical solu-tions has been carried out. The finite element method gives excellent approximations to the analytical solutions for the different boundary conditions considered. The excellent agreements between the numerical and analytical solutions substantiate the capability of the finite element technique to yield accurate results for a series of problems which constitute some of the features of the models developed in chapters 4 and 5. When a sufficient number of elements is used, the approximation is accurate to within 5 % of the exact solution. In most cases, the error is less than 1 % even with a small number of Figure 3.26: Comparison of finite element solutions with 1-dimensional analytical solu-tions for transient convection at r = 0.001 m 71 elements if the change in temperature gradient is not large. For the transient solution, an appropriate time-step size is required to obtain the optimal finite element approxi-mation. A small time-step size is necessary to account for a large change in slope of the temperature-time curve. Chapter 4 Effect of a Magnetic Field on Heat Transfer in a Germanium Melt It is well established that heat transfer in a liquid is markedly influenced by fluid flow due to natural convection and, in some cases, forced convection. The effective thermal conductivity of the liquid, which accounts for the heat transfer due to fluid flow, can be several times the 'stagnant' thermal conductivity of the liquid, the amount depending on the extent of the flow in the liquid. The effective thermal conductivity in liquid tin under the influence of natural convection can be as high as ten times the 'stagnant' thermal conductivity [14]. A n applied magnetic field suppresses the fluid flow in liquid semiconductors. Thus, the magnetic field reduces the effective thermal conductivity of the melt. The effect of a vertically applied magnetic field has been investigated by utilizing temperature measurements within a germanium (Ge) melt in a Czochralski growth con-figuration and mathematical modelling of heat transfer. Specifically, the effective thermal conductivity of the melt was so adjusted that the calculated radial temperature gradients match with the measurements. 4.1 Temperature Measurements and Radial Temperature Gradients Temperature measurements in a Ge melt and related information have been reported in reference [87]. The temperature measurements were carried out in a melt contained in a graphite crucible inside a Czochralski crystal puller (Hamco crystal puller) using 72 73 chromel-alumel (type K) thermocouples. The mass of the melt was 5.499 kg, and the in-ner diameter of the crucible was 0.1512 m. The steady-state temperature distribution at a given melt depth was first measured without a magnetic field; then, the same measure-ments were repeated with a magnetic field. The magnet current was 500 amperes which produced an equivalent magnetic field strength (B) of 0.099 T (990 Gauss). The crucible rotational rate was 14.4 rpm. In this investigation, the depths to be considered from the experimental measurements are 5, 10, 30 and 40 mm. Measurements at the melt surface are also available, but they are unreliable because a small change in the thermocouple height would result in large temperature changes since the surface temperature gradients are much higher than those below the surface. The thermocouple readings showed that there were temperature fluctuations in the melt over time. It was observed that these fluctuations occurred whether the crucible was rotating or not. The extent of the fluctuations in the rotating crucible seemed to be different from that of the stationary crucible, but the amplitude of the temperature fluctuations was clearly affected by the presence of the magnetic field. Without the magnetic field, the amplitude of the temperature fluctuations was roughly 3 °C, but under the influence of the magnetic field it was approximately 2 ° C Thus, the magnetic field dampened the fluctuations of temperature in the melt. The temperature measurements were plotted against radial distance from the center of the crucible at each melt depth using the mean temperature obtained from the fluctuating temperature data. The points were linearly fitted. The deviation from the fitted line is within 10 °C without the magnetic field, and 30 °C with the magnetic field. The data points and the fitted lines are shown in Figure 4.27a to 4.27d for the case without the magnetic field, and Figure 4.28a to 4.28d for the case with the magnetic field. Radial temperature gradients at the four melt depths were taken from the slopes of the lines, as presented in Figure 4.29 for both cases. From Figure 4.29, the applied 74 1050 1040 1030 5 mm below melt surface B — 0 T data point linear fitted line 1020 •4—' y ^ § O 1010 ex ° 1000 crucible rotation = 14.4 rpm 990 980 970 0.02 0.04 Radial Distance from Center ( m ) 0.06 Figure 4.27a: Measured temperatures at 5 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 75 1050 1040 1030 1020 3 _ ° IOIO a . <D H 1000 990 980 crucible rotation = 14.4 rpm 10 mm below melt surface B = 0 T data point linear fitted line I • 970 0.02 0.04 Radial Distance from Center (m) 0.06 Figure 4.27b: Measured temperatures at 10 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 76 1050 1040 -1030 -1020 -t-i S3 O IOIO ex, o B ^ 1000 990 980 970 1 1 1 1 1 --— • — -I • - crucible rotation = 14.4 rpm -30 mm below melt surface B = 0 T -• data point linear fitted line -1 1 I 1 I 0.02 0.04 Radial Distance from Center (m) 0.06 Figure 4.27c: Measured temperatures at 30 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 77 1050 1040 1030 1020 -<D £ 1010 E-i 1000 crucible rotation = 14.4 rpm 990 980 40 mm below melt surface B = 0 T • data point linear fitted line 970 0 JL _L 0.02 0.04 Radial Distance from Center (m) 0.06 Figure 4.27d: Measured temperatures at 40 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under no applied magnetic field [87] 78 1050 1040 1030 1020 V-i S3 ^ IOIO <D 1000 5 mm below melt surface B = 0.099 T • data point linear fitted line crucible rotation = 14.4 rpm 990 980 970 0.02 0.04 Radial Distance from Center (m) 0.06 Figure 4.28a: Measured temperatures at 5 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] 79 1050 1040 1030 -10 mm below melt surface B = 0.099 T • data point linear fitted line 1020 crucible rotation = 14.4 rpm % t> 1010 B H 1000 990 980 970 0.02 0.04 Radial Distance from Center (m) 0.06 Figure 4.28b: Measured temperatures at 10 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] 80 1050 1040 -1030 " 1020 <L) H 1000 crucible rotation. = 14.4 rpm 990 980 30 mm below melt surface 1? = 0.099 T • data point linear fitted line 970 0.02 0.04 Radial Distance from Center (m) 0.06 Figure 4.28c: Measured temperatures at 30 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] 81 1050 1040 1030 1020 O 1010 H 1000 990 980 970 1 1 1 1 1 -• -- crucible rotation = 14.4 rpm 40 mm below melt surface B = 0.099 T • data point - linear fitted line -1 1 i i i 0 0.02 0.04 Radial Distance from Center (m) 0.06 Figure 4.28d: Measured temperatures at 40 mm below the Ge melt surface in a crucible rotating at 14.4 rpm under an applied magnetic field of 0.099 T [87] 82 700 •4—> c •3 <-( o CD t , S3 •4—' cd & O a , o S — <u H 600 500 400 300 200 100 B = 0T • + -•• B = 0.099 T crucible rotation = 14.4 rpm mcll depth 5 mm 10 mm ihcX , n coM 3U mm — 40 mm T,(i) Gemcb ho( * 1 « 0 A l U d u l Temperature Gradient = TiMg»iM 0 0.005 0.015 0.025 Melt Depth (m) 0.035 Figure 4.29: Radial temperature gradients in the Ge melt at the four melt depths in a crucible rotating at 14.4 rpm with and without an applied magnetic field of 0.099 T [87] 83 magnetic field clearly increased the radial temperature gradient by at least a factor of 2; in this case, the effect was the largest near the melt surface and became smaller towards the bottom of the melt. 4.2 Mathematical Model A steady-state conduction-dominated heat-transfer mathematical model was used to es-timate the influence of the applied magnetic field on the effective thermal conductivity of the melt. The model treats the melt as a stationary solid since the actual fluid motion in the melt is not taken into account. 4.2.1 Assumptions The model assumes an axisymmetrical temperature field within the melt. Such an as-sumption is often made in modelling two-dimensional heat flow in the Czochralski crystal growth system. The heat-transfer properties are taken to be temperature-independent. This is done because the heat-transfer properties of the system are not readily available. Assuming constant properties also simplifies the numerical procedure required in obtaining the solutions; if the properties were temperature-dependent, then the resulting system of equations would be non-linear and thus an additional numerical technique such as the Newton-Raphson method, which iterates the solution until convergence criterion is met, would be required in solving the equations; consequently, more time and computational cost would be necessary. On the other hand, such simplification limits the predictive capability of the model. 84 4.2.2 Domain and Boundary Conditions The domain of the model is shown in Figure 4.30. The dimensions are 0.0756 m in radius and 0.0575 m in height. The dimensions correspond to a melt of 5.499 kg contained in a crucible having an inner diameter of 0.1512 m used in the measurements. Within the domain, the temperature field is governed by u / dr\ <yr _ Q r dr \ dr J dz2 At the axis of symmetry, the radial temperature gradient is set to zero to ensure an axisymmetrical temperature distribution. The side surface of the melt is fixed at either a constant temperature or a linear temperature profile. The bottom surface of the melt has two boundary conditions: insulating and non-insulating. The top surface exchanges heat with the ambient surrounding by convection and radiation to an assumed constant ambient temperature. With the boundary conditions, the following equations hold. At the axis of symmetry (centre line), dT * = o («*> Along the side surface, T = Twall (4.35) where Twau is the fixed temperature on the side surface. At the top surface of the melt, dT - km — = hto (T - Ta) (4.36) where km is the effective thermal conductivity of the melt, and hto is the overall heat-transfer coefficient between the top surface and the ambient surrounding. Both the convectice and radiative heat-transfer coefficients are assumed to be accounted for in hto. Along the bottom surface, dT & = 0 <4-37> 85 axis of symmetry heat exchange with ambient surrounding t t t t t t t t . t germanium melt ( fc. ) insulation or U H U U ! heat exchange with ambient surrounding (ho) melt surface specified wall temperature _melt bottom scale 0.01 m Figure 4.30: IdeaHzed Domain of the mathematical model for investigating the influence of the applied magnetic field on the effective thermal conductivity of the Ge melt 86 for insulated bottom, and dT -km — = h b o ( T - T a ) (4.38) for non-insulated bottom, where hbo is the overall heat-transfer coefficient between the bottom surface and the ambient surrounding. In the model, the extent of heat transfer is determined by the effective thermal con-ductivity of the melt, the overall heat-transfer coefficients, and the differences between the melt surface temperatures and the ambient temperatures. The number of heat-transfer parameters is kept to a minimum because they are not known well in the real system. 4.2.3 Finite Element Mesh of the Domain Four different finite element meshes were tested to check the accuracy of the resulting temperature distribution, as shown in Figure 4.31. The finite element meshes are su-perimposed on the corresponding temperature fields. The numbers of elements used are 2 x 3, 4 x 6, 8 x 12 and 12 x 18. The results clearly illustrate the convergence of the temperature field with decreasing mesh size. As expected, a smaller mesh size yields smoother contour lines, indicative of more accurate solutions. Decreasing the mesh size further does not give any visible change in the contour lines. For most accurate results, the finest mesh (12 x 18 elements) should be utilized. However, the results obtained using 8 x 12 elements are not much different from those obtained using 12 x 18 elements. For higher efficiency in terms of computational time and cost, 8 x 12 elements will be used for all subsequent calculations. The mesh is shown in Figure 4.32. 87 Figure 4.31: Effects of finite element mesh on temperature field: (a) 2 x 3 elements; (b) 4 x 6 elements; (c) 8 x 12 elements; (d) 12 x 18 elements [unit in °C] 88 top surface axis of symmetry j j i i ;• ; • t ! j j i { s i i ! i | ! | | | i ermani am melt i 1 1 ! | i j i I S i j | j I t L i i : I | ] j side surface bottom surface scale 1 0.01 m Figure 4.32: Finite element mesh of the mathematical model for investigating the in-fluence of the applied magnetic field on the effective thermal conductivity of the Ge melt 89 Parameter Reference Value Twaii 1000 °C T 600 °C B i t o 0.0575 Bibo 0.0 Table 4.7: Reference parameter values used in the sensitivity analysis 4.2.4 Sensi t iv i ty Ana lys i s of the Parameters There are four essential parameters in the model which affect the temperature field: the temperature at the side surface of the melt Twau, the ambient temperature T„, the Biot number at the top surface B i t o , and the Biot number at the bottom surface Bibo- The Biot numbers are denned as Bito = h-~^ (4.39) and Bi^ = (4.40) where Lm is the vertical length of the melt. The set of reference values is listed in Table 4.7. The effects of the parameters are examined by varying the value of the parameter being analyzed and holding others at the reference values. The temperature field resulting from varying each of the values of the parameters is compared with the reference case in two-dimensional contour plots. The contour interval is 2 °C. 4.2.4.1 Side Surface Temperature, Twau Figures 4.33a and 4.33b show the effects of constant Twau. The value is changed from the reference value (1000 °C) to 1050 °C. It has a dramatic influence on the temperature field even though the overall appearance of the temperature field is similar. B y varying 91 Twau, the temperature distribution in the melt shifts toward the value of Twau. The absolute temperature increase in the melt is close'to the absolute temperature increase in Twan. The difference between the maximum and minimum temperature within the melt remains essentially unchanged. The radial temperature gradient is altered slightly, becoming steeper with higher Twau; this is because the higher surface temperature loses more heat to the ambient surrounding which has a lower temperature. Instead of using a constant Twau, a linear temperature profile is imposed on the side surface of the melt. The imposed temperature decreases linearly from 1025 to 975 °C as it goes up along the side surface. Figures 4.34a and 4.34b compare the resulting temperature fields for a constant and a linear Twau, respectively. The temperature fields are distincly different. The radial temperature gradient is higher than that of the reference case near the bottom, decreases towards the middle height, and becomes negative near the top surface of the melt. Such a temperature distribution would not be realistic in a melt heated from the side wall because the radial temperature gradient in such a system should always be positive. 4.2.4.2 Ambient Temperature, Ta The results for two different ambient temperatures, 600 °C (reference) and 800 °C, are shown in Figures 4.35a and 4.35b. The temperature distribution is not much affected by the change in Ta, as can be seen from the figures; increasing Ta by 200 °C only increases the melt temperature by a maximum of 8 °C, and the change in temperature gets smaller towards the bottom surface. Thus, the effect of ambient temperature is relatively insignificant compared to that of Twau. However, the change in the radial temperature gradient is more pronounced in this case compared to the case of varying Twall-Figure 4.35: Effects of Ta: (a) reference; (b) Ta = 800 [unit in °C] 94 4.2.4.3 Top Surface B i o t N u m b e r , B i t o The sensitivity of the temperature distribution to this parameter is illustrated in Figures 4.36a and 4.36b. Using a value twice the reference value changes the temperature field dramatically, reducing the top surface temperature by about 14 °C. The radial tem-perature gradient is markedly increased as B i t o becomes larger, especially near the top surface. Compared to the case of varying T a , doubling B i t o is more effective in changing the radial temperature gradient of the melt. 4.2.4.4 B o t t o m Surface B i o t N u m b e r , Bibo The effects of Bibo are displayed in Figures 4.37a and 4.37b. The value is changed from the reference value (insulated bottom) to 0.0575. The bottom ambient temperature used in the calculations is 600 ° C As can be seen from the figures, increasing the value of Bibo increases the overall radial temperature gradient in the melt, particularly near the bottom surface. This is because heat is also lost through the bottom surface in addition to the heat loss at the top. 4.2.4.5 Rela t ive Sensitivities of the Parameters To examine the relative effects of the parameters, the radial temperature profiles at the middle height (0.02875 m) of the melt from the sensitivity analysis results are plotted in Figure 4.38. As shown in the figure, increasing Twau essentially shifts the temperature profile upwards without significantly affecting the radial temperature gradient. On the other hand, Ta and the Biot numbers influence the radial temperature gradient more significantly; the Biot numbers have a greater overall effect than Ta. ure 4.37: Effects of Bibo'- (a) reference; (b) Bibo = 0.0575 [unit in ° 97 Figure 4.38: Relative effects of the model parameters (temperature profiles at the middle height of the melt) 98 4.3 Determination of the Magnetic Effect on the Effective Thermal Conduc-tivity of the Ge Melt In estimating the effect of the magnetic field on the effective thermal conductivity of the melt, the bottom of the melt is assumed to be insulated first. The case of non-insulated bottom is then employed to carry out the same analysis, and the results from this case are compared with those obtained assuming bottom insulation. 4.3.1 Estimation of the Biot Numbers In order to use realistic values of the Biot numbers in the fitting process, the values should be estimated from the information available in the measurements [87]. Heating the melt by turning the power to 100 % yielded a constant rate of temperature increase in the melt of 0.68 °C/s. Cooling it by turning the power off gave a constant rate of temperature decrease of 0.29 °C/s. The heating rate can be converted to the net rate of heat input to the melt as follows: jrp qin = rnGe C P G c — = 5.499 x 376 x 0.68 ss 1406 W at The net rate of heat input to the melt is far lower than the total power input to the puller (roughly 40000 W) , indicating that most of the power was lost to the surroundings. Similarly, the cooling rate can be used to obtain a rough estimate of the rate of heat loss from the melt: dT qout = rnGe C P G e — = 5.499 x 376 x 0.29 « 600 W at This can be taken as the total rate of heat loss through all the surfaces of the melt. To estimate the rate of heat loss from the top surface of the melt only, the following conditions are assumed. The heat is lost by radiation and convection, where the convec-tive rate of heat loss, qc, is assumed to be 10 % of the radiative rate of heat loss, qr. The 99 melt surface temperature, Ts, is fixed at 1000 °C, and the radiative ambient temperature, Tar, is given a value of 100 °C. Assuming an emissivity, e, of 0.1 which is reasonable for a shiny melt surface, the rate of heat loss from the top surface having an area Ato is q = qr + qc = Ato 1.1 {ea (Ts4 - Tar4)} « 292 W This rate of heat loss corresponds to a hto value of 54.2 W/m2 °C for an assumed ambient temperature of 700 °C in the fitting process. The 'stagnant' thermal conductivity of Ge is approximately 17 W/m.°C at 1000 °C. Therefore the top surface Biot number, Bito, is roughly 0.1833. This value is probably an overestimate because the effective thermal conductivity of the melt may be significantly higher than the 'stagnant' thermal conductivity used due to convection in the liquid. Thus, the Biot number used in the fitting process should be smaller. Since the total estimated rate of heat loss through all the surfaces of the melt (600 W) is about twice the estimated rate of heat loss through the top surface (292 W), the rate of heat loss through the bottom surface, which was mostly by conduction through the crucible, is presumably of the same order of magnitude as that through the top surface, if the rate of heat loss through the side surface is assumed to be much lower than that through the bottom surface of the melt. The side surface is nearer to the heater which tends to have a higher ambient temperature, resulting in a lower rate of heat dissipation through this surface. Hence, the value of the bottom surface Biot number, Bibo, should be comparable to the value of Bito-Such a large heat loss through the bottom surface from the calculations is supported by the experimental evidence. The radial temperature gradient increased rapidly towards the bottom of the melt, indicating a large amount of heat being lost through this surface. Thus, more realistic modelling should take the heat loss at the bottom surface into account by applying a non-insulating boundary condition. 100 4.3.2 Insulated Bottom Surface In this case, hbo is zero. For the purpose of simplifying the analysis, Twau is fixed at 1025 °C and Ta is set to 700 °C. The only parameter being adjusted to fit the radial temperature gradients is Bito. The side surface temperature and the ambient temperature are not available in the experimental data. In addition, the sensitivity analysis shows that the radial temperature gradient is more affected by the value of Bito. Since the objective is to investigate the effect of the applied magnetic field on the effective thermal conductivity of the melt, it is necessary that Bito be adjusted to determine the relative change in the value of the effective thermal conductivity. In calculating the measured temperature distribution under no magnetic field, a value of 0.0345 for Bito turns out to produce a reasonable fit to the temperature distribution in the melt and the corresponding radial temperature gradients under normal condition (no magnetic field), given the fixed values of Twau and Ta. It is noteworthy that this value is smaller than the estimated Biot number from the heat loss calculations, which is anticipated. Although the values of the parameters could be adjusted to obtain a slightly better fit, the improvement would not be significant; moreover, the uncertainties involved in the experimental measurements and the artificial boundary conditions used do not seem to justify the effort to gain such minor improvements. Besides, the assumptions made in the model set a limit to the extent of the overall improvement of the fit; improving one part of the fit is accompanied by degradation of the fit in another part. Figure 4.39 shows the calculated and measured (fitted) temperature profiles at 5, 10, 30 and 40 mm below the melt surface. There are several observations worth noting. First, since the side surface of the melt in the model is fixed at a constant temperature whereas it was not the case in the experimental system, the calculated lines are shifted to coincide with the fitted (experimental) lines at the side surface. The amounts the 101 Figure 4.39: Comparison of the measured and calculated temperature profiles at 5,10, 30 and 40 mm below the melt surface under no magnetic field assuming bottom insulation 102 calculated lines are shifted are indicated by the values (in °C) shown in the figure near the ends of the lines. A positive number indicates an upward shift of the calculated line, and a negative number a downward shift. Second, since the model does not take into account the fluid motion in the melt, a quantitative fit to the measurements is beyond the capability of the model. Even complex fluid flow models which encompass both natural convection (due to gravity) and forced convection (due to crucible rotation) in the melt do not guarantee 'good' quantitative predictions. Third, the experimental measurements were linearly fitted, and have a large scatter band, with a deviation up to 10 °C under the normal condition and up to 30 °C under the applied magnetic field. Thus, the seemingly large discrepancies between the calculated and measured lines are in fact not unreasonable since the absolute temperature difference between them is less than 30 °C, which is not too far off the scatter band. Lastly, the absolute temperature distribution in the experimental system is known to be dependent upon several factors such as ambient gas and ambient gas pressure [32], making quantitative predictions more difficult. However, for the present set of model boundary conditions, a comparison of the calculated and experimental values is not significant. Moreover, the absolute fit to the temperature distribution is not the main concern here. A more important condition in meeting the objective is to obtain a good fit to the temperature gradients perpendicular to the applied magnetic field as their changes are a direct consequence of the effect of the magnetic field on the temperature distribution of the melt. Since the magnetic field was applied vertically, the radial temperature gradients are a most important indication of the effect of the magnetic field on the effective thermal conductivity of the melt. The measured and calculated radial temperature gradients are shown in Figure 4.40. The radial temperature gradients are calculated by taking the difference between the maximum and minimum temperatures at each melt depth divided by the radius. Figure 4.40: Comparison of the measured and calculated radial temperature gradients as a function of melt depth under no magnetic field assuming bottom insulation 104 The calculated radial temperature gradient can be changed by varying B i t o . In fitting the measured radial temperature gradients under the applied magnetic field, B i t o is increased to 0.2415 which produces a good fit for the radial temperature gradients. Figure 4.41 gives the comparison of the temperature profiles at the four melt depths, and the corresponding radial temperature gradients are compared in Figure 4.42. Figure 4.43 combines the fits of the radial temperature gradients for both cases. The overall agreement is quite satisfactory except near the bottom of the melt. The calcu-lated lines continue to decrease while the measured gradients do not. This discrepancy is a consequence of the assumption of an insulated bottom which can only predict a con-tinuously decreasing radial temperature gradient towards the bottom surface since there is no heat loss through this boundary. Hence, to account for this drawback, heat must be allowed to dissipate at the bottom surface in the mathematical model. Such a model is more representative of the experimental system because the actual boundary is not insulated. 4.3.3 Non-Insulated Bottom Surface The results obtained assuming an insulated bottom surface cannot give an increase in the radial temperature gradient near the bottom surface. This limitation must be overcome by slightly modifying the model to incorporate heat loss at the bottom surface of the melt. In terms of model parameter, Bibo has a positive value. In obtaining the best fit to the radial temperature gradients, both B i t o and Bibo are adjusted holding Twau and Ta at the same values as before. Using the Biot numbers listed in Table 4.8 yields the temperature profiles at the four melt depths shown in Figure 4.44 for the normal case. For the magnetic case, the temperature profiles are displayed in Figure 4.45, and the values of the model parameters are shown in Table 4.9. The comparisons of the radial temperature gradients are presented in Figure 4.46, which 105 Figure 4.41: Comparison of the measured and calculated temperature profiles at 5,10, 30 and 40 mm below the melt surface under an applied magnetic field of 0.099 T assuming bottom insulation 106 Figure 4.42: Comparison of the measured and calculated radial temperature gradients as a function of melt depth under an applied magnetic field of 0.099 T assuming bottom insulation 107 700 c T3 o § — 400 1 £ & o a - 3 o o H \ X \ no magnetic field • measured 600 X x x ( B = 0 T ) calculated x •v X magnetic field o measured 500 « \ X N X ( J B = 0.099 T ) calculated T3 200 100 •km " 7 \ • 0.02 0.04 Melt Depth ( m ) 0.06 Figure 4.43: Comparison of the measured and calculated radial temperature gradients as a function of melt depth with and without an applied magnetic field of 0.099 T assuming bottom insulation 108 0 0.02 0.04 0.06 0.08 Radial Distance from Center (m) Figure 4.44: Comparison of the measured and calculated temperature profiles at 5, 10, 30 and 40 mm below the melt surface under no magnetic field without bottom insulation 109 0 0.02 0.04 0.06 0.08 Radial Distance from Center ( m ) Figure 4.45: Comparison of the measured and calculated temperature profiles at 5, 10, 30 and 40 mm below the melt surface under an applied magnetic field of 0.099 T without bottom insulation 110 Parameter Value Twaii 1025 °C T 700 °C Bito 0.0246 Bi^ 0.0164 Table 4.8: Parameter values used in the calculations for the normal case Parameter Value Ttuall 1025 °C Ta 700 °C B i t o 0.1725 Bibo 0.1150 Table 4.9: Parameter values used in the calculations for the magnetic case also includes the calculated radial temperature gradients from the model assuming an insulated bottom surface. Prom Figure 4.46, it can be seen that the radial temperature gradient increases near the bottom surface when heat is allowed to dissipate at this boundary. Since the increase of the radial temperature gradient was observed experimentally, this model is a more realistic representation of the actual experimental system. It is worth noting that the values of the Biot numbers used in fitting the measurements under the applied magnetic field are only slightly smaller than the values estimated from the heat loss calculations, done using the 'stagnant' thermal conductivity of the melt. This seems to indicate that the applied magnetic field reduced the effective thermal conductivity of the melt to the value approaching the thermal conductivity of a 'stagnant' melt. I l l a , ° H 700 600 K * 500 400 300 200 100 magnetic field ( B = 0.099 T ) A measured calculated with insulated bottom calculated with non-insulated bottom 7 no magnetic field( B = 0 T )" measured calculated with insulated bottom calculated with non-insulated bottom 0.02 0.04 0.06 Melt Depth ( m ) Figure 4.46: Comparison of the measured and calculated radial temperature gradients as a function of melt depth with and without an applied magnetic field of 0.099 T for both the insulating and non-insulating bottom boundary conditions 112 No Magnetic Field Magnetic Field of 0.099 T B i t o Bi^ 0.0345 0 0.2415 0 Table 4.10: Biot numbers used in fitting the measurements assuming bottom insulation No Magnetic Field Magnetic Field of 0.099 T Bito Bibo 0.0246 0.0164 0.1725 0.1150 Table 4.11: Biot numbers used in fitting the measurements accounting for bottom heat loss 4.3.4 Reduction Factor in the Effective Thermal Conductivity of the Ge Melt The values of the Biot numbers for both the normal and magnetic cases are compared in Table 4.10 for the insulated bottom surface and Table 4.11 for the non-insulated bottom surface. The Biot numbers have been defined previously as hto Lm B i t o = Bibo = km hbo Lm Assuming that the overall heat-transfer coefficients are the same for both the normal condition and the magnetic condition, the ratio of the Biot number for the magnetic case to that for the normal case is inversely proportional to the corresponding effective thermal conductivity ratio. The mathematical equation is Bi magnet k m normal Bir k„ (4.41) '•normal ^m magnet Using equation 4.41, the effective thermal conductivity ratio for the case of bottom 113 insulation is km normal _ 0-2415 km magnet 0.0345 For the case of non-insulated bottom, it is 7.0 The resulting ratio is 7.0 regardless of the different boundary conditions used at the bottom surface of the melt. Hence, the effective thermal conductivity of the melt under no magnetic field, kmnormah is 7x that under the applied magnetic field, kmmagnet. This corresponds to a reduction factor of 7 in the effective thermal conductivity of the Ge melt due to the effect of an applied magnetic field of 0.099 T. 4.4 Discussion on the Reduction Factor The results for both the insulated and non-insulated bottom cases turned out to be the same. The estimated reduction factor due to the applied magnetic field was 7. The number was derived by fitting the model results to the measurements. Consequently, depending on the accuracy of the measurements and the decision of what the 'best fit' was, the resulting number could be smaller or larger. It should also be noted that the resulting factor may be model-dependent. The factor 7 is a number obtained with a steady-state heat-transfer mathematical model using sim-plified boundary conditions and constant heat-transfer properties. More complex bound-ary conditions and temperature-dependent heat-transfer properties would give different temperature distributions and radial temperature gradients. Thus, a different reduction factor in the effective thermal conductivity would probably be obtained employing more complex models. It also depends on the set of values used for Twau and Ta. A different set of values might yield different results. Moreover, the reduction factor may not be the same everywhere 114 in the melt; it may be temperature-dependent; the non-uniformity of the magnetic force within the melt may also lead to positional dependence of the reduction factor. Nevertheless, since the magnetic field reduced the fluid flow in the melt, mixing of the melt was decreased and accordingly the amount of heat transferred by the fluid flow was reduced. Assuming that the overall heat-transfer coefficients characterizing the rate of heat loss from the melt to the surrounding under no magnetic field are the same as those under the applied magnetic field, the suppression of the melt flow by the magnetic field is thermally equivalent to reducing the effective thermal conductivity of the melt. Fluid flow in a CZ melt is very complicated and loosely defined; even a small density difference in the melt causes complex flow patterns which change appreciably with time. Extensive fluid flow modelling is often not justified because the model predictions are hardly representative of the flow in a real growth system. Hence, simpler steady-state conduction-dominated modelling can sometimes be a more efficient alternative to achieve the same results. It is in this purely conductive heat-flow modelling of electrically conduc-tive melts that the reduction factor of the effective melt thermal conductivity can prove to be useful as a first approximation to account for the effect of an applied magnetic field on heat transfer in the melt. A similar approach has been used in mathematical modelling of solidification of met-als. The liquid pool was assumed to be completely mixed everywhere and to act as a conducting solid. To account for the enhancement in heat transfer due to convection in the liquid, an effective thermal conductivity was used. The value of the effective thermal conductivity was taken to be seven to ten times the 'stagnant' thermal conductivity of the liquid [14]. As the effective thermal conductivity of liquid tin due to natural convection only can be as high as ten times the 'stagnant' thermal conductivity [14], the reduction factor of seven by the applied magnetic field of 0.099 T has not eliminated the melt flow especially 115 the flow due to the crucible rotation. A stronger magnetic field is necessary to further reduce the melt flow. In the limit, the melt flow would be completely eliminated, but such a case is practically impossible. Nevertheless, with a large magnetic field, the results obtained using a purely conductive heat-transfer analysis compare fairly well with the experimental results for silicon growth [66]; such a good agreement seems to indicate that the melt convection is not significant when a large magnetic field is present. 4.5 Relative Thermal Resistance Analysis of the Mathematical Model A relative extent of the different mechanisms of heat transfer can be analyzed using thermal resistances. The thermal resistance circuit of the model domain is shown in Figure 4.47. There are two, convective/radiative resistances: at the top surface, Rhtoi a n ( l at the bottom surface, Rhboi ^ o r * n e c a s e 0 1" non-insulated bottom. There is a conductive resistance through the Ge melt, RkGe-The resistances are defined as *»" i6l ( 4 4 2 ) = fclXT ( 4 4 3 ) "•bo -rt-m where Am is the horizontal cross-sectional area of the melt. Using the heat-transfer properties at 1000 ° C from the heat loss calculations, 116 convection/radiation from top surface of domain conduction through germanium melt Ge top surface <3e melt side surface bottom surface convection/radiation from bottom surface of domain scale 0.01 m Figure 4.47: Thermal resistance circuit in the model domain used for investigating the influence of the magnetic field on the effective thermal conductivity of the Ge melt 117 convection/radiation from domain boundary 1.025 °C/W 1.474 °C/W conduction through germanium melt •0.1879 ° C / W Figure 4.48: Thermal circuit with relative resistances of the mathematical model used for investigating the influence of the magnetic field on the effective thermal conductivity of the Ge melt 118 The thermal circuit as well as the resistance values is displayed in Figure 4.48. The total resistances are heat loss to surrounding: i?„ = J Z h t o = 0.6046 °C/W heat conducting through liquid germanium: i2jtGe = 0.1879 °C/W Thus, the conductive resistance is several times smaller than the convective/radiative resistance. Since the Biot numbers for the magnetic case are close to those obtained from the heat loss calculations, the above resistances are representative of the magnetic case. For the case of no magnetic field, the conductive resistance is 7x lower than that for the magnetic case because the effective thermal conductivity under no magnetic field is larger by a factor of 7 from the analysis of the magnetic effect on the effective thermal conductivity of the melt. Chapter 5 Determination of the Thermal Conductivity of a Pyrolytic Boron Nitride Crucible The method employed in the present investigation to obtain the thermal conductivity of a pyrolytic boron nitride crucible is not a standard procedure for determining the thermal conductivity of a material since suitable apparatus was not available. The apparatus used in this investigation was developed for this project, and is relatively simple. However, a mathematical model of the heat flow in the experimental system is necessary in order to analyze the experimental results, which is time-consuming. The procedure used has not been reported and therefore represents a novel technique which requires testing and characterizing for reliability and accuracy of the derived thermal conductivity. The thermal conductivity of a pyrolytic boron nitride (PBN) crucible in the thickness direction is determined utilizing temperature measurements and mathematical modelling. Reference specimens of known thermal conductivities are used to calibrate the model parameters. The calibrated model is employed to obtain the thermal conductivity of the P B N crucible. The determination is accomplished by matching the calculated (predicted) temperature-time curves with the experimental measurements. 5.1 Description of the Material The P B N material examined is from a used P B N crucible supplied by Johnson-Matthey Electronics. The crucible has been used to grow GaAs single crystals. Consequently, it has undergone heating and cooling cycles that caused expansion and contraction of the 119 120 crucible. Such dimensional changes normally give rise to differential stresses in a crucible which in turn shorten the crucible life. A n occasional failure of P B N crucibles is uncon-trolled delamination of layers. The delamination occurs as a result of the highly-oriented crystallographic texture and the highly-anisotropic properties of P B N [55]. Figure 5.49 shows a broken edge of the crucible. The different delaminated layers which make up the thickness of the crucible are clearly observed. Since the material examined is from a used crucible, the material properties may not be the same as those of a fresh crucible especially when dimensional changes have occurred over a long period of time. In practice, the material in use is changing continually. Hence, the thermal conductivity of this material is probably different to some extent from that of an unused P B N crucible. 5.2 Methods of Measurement A n experimental apparatus was designed to measure the rate of heat flow across the specimen to determine the thermal conductivity. Temperature measurements were taken at different positions in order to characterize the heat flow in the system. Initially, a steady-state procedure was attempted. Unfortunately, this method failed to produce adequate results. A transient method was then used by modifying the apparatus and procedure. 5.2.1 Steady-State Method In this method, the idea is to measure the steady-state rate of heat flow across the specimen. To a first approximation, knowing the heat-flow rate, the thermal conductivity of the specimen can be calculated using the Fourier's law of heat conduction. 121 Figure 5.49: A broken edge of the PBN crucible revealing the layers that make up its thickness (magnification = 16 x) 122 5.2.1.1 Experimental Apparatus The experimental apparatus consisted of a boat made of boron nitride, as shown in Figure 5.50. Two specimens were inserted into the two grooves in the boat separating it into 3 compartments, the middle being smaller than the two equal-size outer compartments. A resistance heater made of platinum wire was placed inside the middle compartment to provide the power input to the system as shown in Figure 5.51. A liquid conducting medium was used to give a better heat transfer from the heater to the specimens. Liquid gallium (Ga) was chosen for the purpose as it has a low melting point and thus is easier to handle. Three 0.02-inch chromel-alumel thermocouples (type K ) sleeved in alumina tubes were placed in the liquid gallium at positions shown in Figure 5.51. As liquid gallium reacts readily with metals at high temperatures, the tips of the thermocouples were coated with a ceramic paste (ceramabond 571) for protection. Fibrefrax was used to insulate the boat. The whole arrangement was contained in a horizontal tube furnace. Argon gas at a rate of roughly 1.33 x 1 0 - 5 m3/s was used to replace the air inside the tube primarily to prevent oxidation of the liquid gallium. 5.2.1.2 Temperature Measurement Procedure The system was thermally equilibrated at a given furnace temperature. Power was then supplied to the heater to raise the temperature of the liquid gallium in the middle com-partment. Heat flow occurred from this compartment to the two adjacent compartments within the insulation layer. The temperatures in the three thermocouple locations shown in Figure 5.51 were recorded using a potentiometer at discrete time intervals. The tem-peratures were measured relative to an ice/water bath reference cold junction. 123 A - A Cross Section A 4 A A scale 0.01 rn top view Figure 5.50: Boron nitride boat used in the experimental apparatus 124 jNianztube cortnpctjTg ter ; sl( eved :na Figure 5.51: Schematic diagram of the experimental apparatus for the steady-state mea-surements 125 5.2.1.3 Preliminary Experimental Results and Discussion Only a limited number of measurements was obtained using this method because of the difficulties encountered in the measurements. The problems associated with the experimental system are addressed in the next section. The material plates used in this preliminary experiment were the same as the boat material (BN). This experiment was intended to calibrate the whole system and to find out how reliable the experiment was. Since the thermal conductivity of B N is known, the calibration and reliability can be estimated by comparing the resulting experimental thermal conductivity with the known value. Figure 5.52 gives the temperatures of the three thermocouple locations as a function of time. The temperatures rise initially and then reach steady-state values. The increase in temperature by the heater is roughly 35 °C. Figure 5.53 shows the corresponding rates of temperature increase versus time. T C I has the highest rate as expected, followed by TC2 and then TC3 . After approximately 25 minutes of heating, the rate of temperature increase at T C I drops below those of TC2 and TC3, and finally goes to zero. This occurs because most of the additional heat introduced into the liquid gallium in the middle compartment is lost to the surroundings. The rate of temperature rise in the middle compartment can be estimated assuming that all the heat generated by the heater stays inside the liquid gallium in the middle compartment. The estimated value for dT/dt is 12.75 °C/minute (Appendix C). Com-paring this value of dT/dt to the values in Figure 5.53, it is clear that almost all the heat generated by the heater is lost to the surroundings. If the heat generated by the heater flowed only through the material plates, the temperature rise in the middle compartment would be much higher than the maximum value shown in Figure 5.53. For the rate of temperature increase in the middle compartment shown in Figure 126 Figure 5.52: Temperatures as a function of time for T C I , TC2 and TC3 positioned in the three cells as shown in Figure 5.51 127 i 1 1 1 1 1 1 1 1 1 1 1 1 r 10 14 18 22 26 30 34 38 Time ( minutes) Figure 5.53: Rates of temperature increase as a function of time for T C I , T C 2 and TC3 positioned in the three cells as shown in Figure 5.51 128 5.53 (roughly 1.5 °C/minute) , the corresponding rate of temperature increase in the side compartment can be calculated assuming uni-directional heat flow through the plates. The calculation is given in Appendix D. The rate of temperature increase (dT/dt) for the side compartment is calculated to be approximately 0.461 °C/minute . Compared to the measured value shown in Figure 5.53 which is roughly 1.5 °C/minute , the calculated value is low. It is unclear why the observed rate of temperature increase is higher than the estimated value. 5.2.1.4 Problems Associated with the Steady-State Method There were uncertainties and problems encountered in the steady-state measurements. A major problem was due to insulation. In order to obtain an accurate thermal conductivity value, heat flow, other than across the material plates being tested, must be kept to a minimum. The amount of heat that flows through the plates must be known accurately, which means that the heat generated by the resistance heater in the middle compartment must flow only through the test plates so that the thermal conductivity of the test material can be estimated accurately. Q For the above experimental system, it was difficult to obtain good insulation. The diameter of the quartz tube was only slightly larger than the diameter of the boat. This limited the amount of insulating material that could be used to wrap around the boat. Moreover, it was difficult to know how thick the insulating layer must be to reduce the heat gain from the surrounding to a negligible value. When an electrical current was passed through the Pt wire, not only was the part of the Pt wire which was immersed in the liquid gallium producing heat, but the extended Pt wires which were connected to the power supply via tungsten connecting wires were generating heat also. Thus, the additional heat coming from the extended Pt wires complicates the heat balance of the experimental system. 129 The B N boat used in the experiment has a thermal conductivity which is comparable to the thermal conductivity of the test plate. As a consequence, much of the heat generated by the heater inside the middle compartment was conducted through the B N boat walls. Moreover, since the absolute amount of heat being generated was relatively small, the heat balance of the system was very sensitive to any fraction of heat being lost through the B N boat walls. The heat conducted through the B N boat walls would contribute to the temperature rise in the other compartments which would result in inaccurate thermal conductivity values. The properties of the conducting liquid, such as density and heat capacity, must be known precisely since these properties affect the thermal conductivity values calculated from the temperature measurements. Thus, the liquid gallium must be kept as clean and pure as possible. At high temperatures, the liquid gallium is easily oxidized and therefore the properties are easily altered. Even with the continuous flow of argon, it was difficult to prevent oxidation of the liquid gallium. The extent of oxidation is hard to estimate, and the oxidation effect on the accuracy of measurement is difficult to characterize. A n additional problem associated with liquid gallium is that it reacts readily with metals at high temperatures, dissolving any metal that comes in contact with it. Since the heater and the thermocouples must be dipped in the liquid gallium, they must be protected with a ceramic coating to prevent any direct contact with the liquid gallium. While the ceramic coating protected the heater and the thermocouples from dissolving in the liquid gallium, the influence of the coating on the accuracy of the experiment was not known. The effect of the continuous flush of argon on the heat transfer of the system was also difficult to characterize. It presumably enhanced heat loss, and hence reducing insulation. Another major problem arose from the resistive heating system. The resistance heater was made of platinum wire. The heat was generated by passing a current through the 130 platinum wire. For a poorly insulated system, there would be competition between the heat entering the middle compartment and heat being lost to the surroundings. Thus, the amount of current needed to raise the middle compartment to a certain temperature was dependent on the insulation. Wi th better insulation, less current would be needed to raise the temperature. Since the insulation of the system was poor, it required much more current than the amount needed theoretically. Experimentally a current of 5 amperes was required to produce a temperature rise of 1 °C/minute . Unfortunately, the platinum wire tended to burn out and break with this input current. Better insulation of the system would have increased the life of the heating wire since a lower current could have been used, but this was not easy to achieve with a B N boat and the existing apparatus. Wi th all the uncertainties and problems involved in the experimental system, it was difficult to have confidence in the thermal conductivity values calculated from the mea-surements. But more unfortunately, the frequent failure of the Pt resistance heater prevented any meaningful measurements to be made. It was decided to modify the measuring system and use the transient method to obtain the thermal conductivity. 5.2.2 Transient M e t h o d In this method, the transient rate of heat flow across the test sample is measured. This can be done by measuring the temperature changes in the system as a function of time. Knowing the thermal responses, the thermal diffusivity of the specimen can be deter-mined. The thermal conductivity can be calculated using the density and specific heat capacity of the test sample. 5.2.2.1 Expe r imen ta l Appara tus The experimental apparatus used for the transient measurements is shown in Figure 5.54. In this case, one end of the boat was cooled using argon gas which flows through 131 quartz tube test specimen quartz tube stainless-steel tube chromel-alumel thermocouples sleeved in alumina tubes alumina tube argon gas 1 1 BN boat I | liquid gallium alumina cover fibrefrax insulation ceramic coating cold reference junction data acquisition system scale 0.05 m Figure 5.54: Schematic diagram of the modified experimental apparatus for the transient measurements 132 the stainless-steel tube as shown in the figure. One test sample was used which divides the boat into two compartments. When the liquid gallium in the end compartment is cooled, there is a corresponding temperature change in the liquid gallium in the other compartment. The temperature response depends on the heat transfer in the system, primarily determined by the thermal conductivity of the test sample. Four chromel-alumel thermocouples (type K ) were used to measure the temperature responses at 4 locations in the liquid gallium, as shown in Figure 5.54. The estimated distances between the specimen surfaces and the locations of T C I , TC2 , TC3 and TC4 were 27.5 mm, 2 mm, 2 mm and 31.5 mm, respectively. The whole arrangement was insulated using fibrefrax. The assembly was inserted into a horizontal tube furnace. Argon gas was flushed through the tube furnace at a rate of roughly 1.33 x l O - 5 m 3 / s to provide an inert atmosphere. To reduce the effect of the high thermal conductivity of the B N boat walls, the inside surface of the boat was coated with a ceramic paste to minimize the heat flow through the walls. The tips of the thermocouples and the end of the stainless-steel tube in contact with the liquid gallium were protected with a ceramic coating. 5.2.2.2 Temperature Measurement Procedure The system was initially equilibrated at a given furnace temperature. Once the system was at steady-state thermally, the flow of argon was initiated at a rate of approximately 1.333 x l O - 4 m 3 / s , cooling the liquid gallium in the end cell. After cooling for about 300 seconds, the argon flow was stopped following which the temperature of the system rose as a result of the higher furnace temperature. Each cycle took 30 minutes. During each cycle, the temperatures at the four thermocouple locations were recorded continuously every 4 seconds using a data acquisition system and the readings were stored in disks for subsequent analyses. For each test sample, several experiments at different furnace 133 Test Sample Thickness ( m ) Length ( m ) Width ( m ) Pyrex B N The P B N Crucible 0.00099 0.00066 0.00071 0.0158 0.0158 0.0158 0.0156 0.0156 0.0156 Table 5.12: Dimensions of the test samples temperatures were carried out. A cold reference junction was provided by an ice/water bath. The temperature mea-surements were accurate to within 0.5 °C. Temperature fluctuations were observed due to the electrical noise and fluid motion. The amplitude of the fluctuations was less than 0.5 °C in most cases. 5.2.2.3 Test Materials and Specifications Three test materials were used: (a) borosilicate glass (Pyrex), (b) boron nitride (BN), and (c) the pyrolytic boron nitride (PBN) crucible. Pyrex and B N were tested to calibrate the heat transfer of the system since their thermal conductivities are known. The P B N crucible test material was the material of interest. The test samples were in the form of rectangular plates. The P B N crucible plate was cut from the bottom of the crucible. The B N and P B N crucible specimens were polished to 600 grit to flatten the surfaces. The dimensions of the specimens measured using digital calipers are shown in Table 5.12. 5.2.2.4 Experimental Results and Discussion Temperature-time curves were obtained for a series of experiments for the temperature ranges listed in Table 5.13, Table 5.14 and Table 5.15. The furnace (ambient) tempera-tures shown in the tables have been adjusted using a procedure which will be described Experiment no. Adjusted Furnace (Ambient) Temperature ( °C ) 1 65.9 2 102.8 3 142.0 4 184.1 5 251.8 6 303.1 7 356.8 8 404.5 9 458.3 Table 5.13: Adjusted furnace (ambient) temperatures for the Pyrex test sample Experiment no. Adjusted Furnace (Ambient) Temperature ( °C ) 1 431.4 2 532.9 3 630.2 4 728.7 5 820.1 Table 5.14: Adjusted furnace (ambient) temperatures for the B N test sample 135 Experiment no. Adjusted Furnace (Ambient) Temperature ( °C ) 1 133.9 2 175.3 3 206.8 4 252.1 5 300.3 6 346.5 7 390.8 8 440.5 9 489.8 10 540.7 11 585.6 12 630.7 13 676.7 14 725.9 15 773.7 16 822.0 Table 5.15: Adjusted furnace (ambient) temperatures for the test sample of the P B N crucible in the thickness direction 136 later. Since the temperature responses behave in a similar manner for all the measurements, only typical results are presented and discussed. Results for the three different materi-als are compared to distinguish different features with reference to the relative thermal properties of the materials. 5.2.2.4.1 Measured Results Typical results, taken from experiment no. 9 of Pyrex, are given in Figure 5.55. Four curves are shown corresponding to the temperatures at the locations of thermocouple 1 to 4 for T C I to TC4, respectively. The curves were reproducible to within 5 °C. The differences in the initial temperatures of the four thermocouple readings were due to the furnace temperature gradient. Normally, the temperature at the position of T C I was the highest, followed by that at the position of TC2 , TC3 and then TC4. As expected, the temperature at position 4 dropped faster than those at the other 3 positions during cooling, with the least change at position 1. After cooling stopped, the temperature in the system returned to the initial steady-state distribution after a short period of time. The temperature at position 4 rose at a much faster rate than the rest. This is because the temperature difference was the largest between the ambient surrounding and position 4. Each curve consists of a cooling section, and a reheating section. Considering the cooling section, heat was extracted into the cooling block causing heat to flow along the axial length of the liquid gallium. The furnace temperature, which was higher than the temperature of the liquid gallium, was maintained during cooling, and heat flowed into the liquid gallium from the surrounding primarily in the radial direction. Within the liquid gallium, the radial heat flow was much lower than the axial flow, resulting in effectively uniaxial heat flow along the liquid associated with the heat extraction. Such a 137 340 ' 1 1 t i i i i i i i i i i i 0 200 400 600 800 1000 1200 1400 Time (seconds) Figure 5.55: Typical output of the four thermocouple readings (taken from experiment no. 9 of Pyrex) 138 urn-directional heat extraction is extremely important as it gives a good representation of the rate of heat flow through the specimen in the experimental system. For the reheating sections of the curves in Figure 5.55, the radial heat input from the surrounding became much more significant than the heat flow across the specimen, as observed from the rapid temperature rise of the liquid gallium towards the furnace temperature, which makes analysis of these sections uncertain. Therefore, analysis is confined to the cooling sections. For the experiments at higher furnace temperatures, the temperature drops during cooling at the four positions were larger. This was attributed to a larger difference be-tween the temperature of the argon gas used to cool the liquid gallium and the initial temperature of the liquid. A higher argon gas flow rate also resulted in a larger temper-ature decrease. 5.2.2.4.2 Comparison of the Results for the Three Different Materials The temperature-time curves for the three materials examined are very similar. How-ever, if examined carefully, there are differences with respect to the rate of heat flow through the specimens. Figures 5.56, 5.57 and 5.58 show the temperature-time curves for Pyrex, B N and the P B N crucible in the thickness direction around 450 °C7, for tests no. 9 of Pyrex, 2 of B N and 10 of the P B N crucible as listed in Tables 5.13, 5.14 and 5.15, respectively. The curves have been corrected for the thermocouple errors and furnace gradient, the procedure of which will be described in section 5.4.1. Since most of the heat flow across the test sample occurred during cooling, and not during heating, only the cooling portions of the curves are used. In comparing the three figures, it is assumed that the heat exchange with the ambient surrounding is the same for all cases. It is possible that the insulation was slightly different 139 Figure 5.56: Temperature responses for Pyrex (taken from experiment no. 9 of Pyrex) 140 i 1 1 1 1 1 1 1 1 1 r 0 200 400 600 800 1000 1200 1400 Time (seconds) Figure 5.57: Temperature responses for B N (taken from experiment no. 2 of BN) 141 Figure 5.58: Temperature responses for the P B N crucible in the thickness direction (taken from experiment no. 10 of the P B N crucible) 142 from one case to another, but the differences were likely too small to alter the curves significantly. From the thermal responses shown in the figures, the temperature difference between the positions of T C I and TC2 is small. This difference is mainly controlled by the effective thermal conductivity of the liquid gallium. The significant information related to the thermal property of the test sample is the temperature difference between the positions of TC2 and TC3. The temperature difference between the positions of TC2 and TC3 is determined by two major variables: the thermal diffusivity and the thickness of the test sample. It is apparent that B N has the smallest temperature difference while the P B N crucible has the largest. For the same thickness, it appears that B N has the largest thermal diffusivity among the three materials. But the B N specimen used in the experiments has the smallest thickness which tends to make the temperature difference smaller. One definite observation is that the P B N crucible has a lower thermal diffusivity compared to Pyrex because the P B N crucible has a larger temperature difference and a smaller specimen thickness. Since the thicknesses of the materials are different, no thorough conclusion can be made about the relative thermal diffusivities of the three materials from the temperature difference. A proper analysis of these curves requires a mathematical model which incorporates all the essential variables in order to determine the thermal properties of the unknown material (the P B N crucible in the thickness direction). 5.3 Mathematical Model A two-dimensional transient model was formulated to simulate the experimental process for the transient method. The thermal conductivity was obtained by matching the calcu-lated temperature-time curves with the measured ones. The model simplifies the actual 143 experimental system in the sense that it is two-dimensional, but it takes into account the heat exchange between the system and the surrounding due to imperfect insulation. 5.3.1 Assumptions The model assumes that the system is two-dimensional. The ambient temperature is taken to be constant which is reasonable since the furnace temperature was not affected by the temperature change in the liquid gallium. In addition, the model assumes a domain which consists of the specimen and liquid gallium only. A n one-dimensional model would probably be sufficient for the purpose; in fact, for simpler boundary conditions and constant heat-transfer properties, analytical solutions could be employed. The reasons in using a two-dimensional numerical model instead of an one-dimensional model or analytical solutions are that it allows more complex and re-alistic applied boundary conditions and temperature-dependent heat-transfer properties to be incorporated. The temperature-dependent heat-transfer parameters are calculated using the tem-perature field from the previous time step. This is a justified approach since the time-step size is relatively small and the temperature distribution does not change appreciably within the time interval. The usage of the time-varying input temperature profile warrants an explanation. Since this temperature profile belongs to only one location in the experimental system, applying it along the whole length of the specimen surface as the boundary condition assumes that there exists no temperature gradient along this boundary. This is a plausible assumption because the radial temperature gradient in the high-conducting liquid gallium is small. 144 5.3.2 Domain and Boundary Conditions The domain of the model is shown in Figure 5.59. It represents a horizontal cross-sectional area of the core of the experimental system. Only half of the area is modelled since the actual area is assumed to be symmetrical about the centre line. It comprises of two distinct regions representing the specimen and liquid gallium. The dimensions used in the calculations are 0.005 m x t for the specimen and 0.005 m x 0.028 m for the liquid gallium, which correspond to the actual size of the experimental system, t is the thickness of the specimen which may vary from one specimen to another. Within the domain, the temperature distribution is governed by Bx* + dy*~ cc dt ( 5 - 4 5 ) The boundary conditions are as follows: 1. at the axis of symmetry (centre line), thermal insulation is used as required by the model domain, i.e. & = 0 <5-4 6> 2. along the surface of the specimen, a time-varying input temperature profile is im-posed, i.e. T = f(time) (5.47) 3. along the other boundaries, the system exchanges heat with the ambient surround-ing by - f c ( n . V T ) = h(T-Ta) (5.48) where h is an overall heat-transfer coefficient. The last boundary condition can be changed to an insulating boundary condition by setting the overall heat-transfer coefficient to zero. Such a boundary condition would only be applicable if the experimental system were perfectly insulated. 145 heat exchange with ambient surrounding t t t t liquid gallium test specimen scale 0.005 m axis of symmetry specimen thickness = t O TCI TC2 O heat exchange with . ambient surrounding x time-varying input temperature Figure 5.59: Idealized Domain of the Mathematical Model for the determination of the thermal conductivity of the P B N crucible 146 Parameter Reference Value Unit number of finite elements 1 x 16 elements time-step size 4.0 s thermal diffusivity of specimen, oca io- 7 m2/s effective thermal diffusivity of liquid gallium, ctGa 10" 5 m?/s overall heat-transfer coefficient, h 0.0 W/m2.°C specimen thickness, t 0.001 m Table 5.16: Reference values used in the sensitivity analysis of the model parameters 5.3.3 Finite Element Mesh of the Domain The finite element mesh employed for all the calculations is shown in Figure 5.60. The whole domain is discretized into 16 elements: 2 in the specimen and 14 in the liquid gallium. Since the temperature gradient across the specimen thickness is higher than that across the liquid gallium, smaller-sized elements are utilized in the specimen domain. 5.3.4 Results of Sensitivity Analysis and Discussion Sensitivity analysis is helpful in understanding the relative effects of the model parameters on the calculated temperature-time curves. A set of reference values is used for the sensitivity analysis of the model parameters. The values are listed in Table 5.16. The time-varying input temperature used for the model is given in Figure 5.61. The number of finite elements and the time-step size shown in the table are optimal. Increasing the number of elements 2 times does not alter the calculated results to any significant extent, and neither does a smaller time-step size improve the results. The temperature-time curves shown below are at 2 mm from specimen/liquid gallium interface (corresponding to the location of TC2) and 2 mm from the other end of the liquid (corresponding to the location of T C I ) at the axis of the domain. The calculated results show a temperature difference of less than 0.5 °C between the chosen location 147 TCl axis of symmetry y T C 2 e specimen thickness = t liquid gallium x test specimen Figure 5.60: Finite element mesh of the model domain for the determination of the thermal conductivity of the P B N crucible 148 Figure 5.61: Time-varying input temperature profile used for the sensitivity analysis (taken from experiment no. 14 of the P B N crucible) 149 and any location within a radius of 2 mm. Thus, the location chosen is not as critical as it might be anticipated. In Figure 5.62, the thermal diffusivity of specimen is reduced and increased by 10 times the reference value to 1 0 - 8 and 1 0 - 6 m2/s, respectively. As can be seen from the figure, this parameter significantly impacts the thermal responses at the two locations. Reducing as inhibits the rate of heat flow through the specimen. This in turn results in a smaller temperature decrease in the liquid and consequently the temperature distribution is more uniform as indicated by the smaller difference between the two locations. On the other hand, the rate of heat flow across the specimen increases when cta is increased. This means that the temperature in the liquid is more readily affected by the temperature change on the other side of the specimen. While the temperature in the liquid stays relatively unchanged irrespective of the changing input temperature in the other case where cts is decreased, the temperature in this case follows the input temperature very closely indicating a faster rate of heat flow across the specimen. The effects of the effective thermal diffusivity of liquid gallium are displayed in Fig-ure 5.63. ctGa influences the uniformity of the temperature distribution in the liquid. As shown in the figure, a lower « G a value increases the temperature gradient in the liquid, while a higher value decreases it. When the value is decreased by 10 times, the temper-ature gradient is increased by several times; reducing it by 10 times virtually eleminates the temperature gradient in the liquid. A larger ctGa value means a faster transfer of heat in the liquid which tends to homogenize the temperature distribution while a smaller ctGa gives the opposite effect. The overall heat-transfer coefficient has a significant effect on the temperature field of the liquid. Figure 5.64 shows the temperature change as h is changed to a positive value. The reference case assumes an insulated boundary, i.e. h = 0. When h is increased to some positive value, there is an exchange of heat between the liquid and the ambient 150 T I I 1 ^(21 ' ' ' ' ' ' ' ' I I I I I ' ' ' I I I I I < 1 _ 0 200 400 600 800 1000 1200 1400 Time (seconds) Figure 5.62: Effects of the thermal diffusivity of specimen 151 Figure 5.63: Effects of the thermal diffusivity of liquid gallium 152 Figure 5.64: Effects of the overall heat-transfer coefficient 153 surrounding. Since the ambient temperature is higher than the liquid temperature, the liquid receives heat from the ambient surrounding resulting in a temperature increase as indicated by the upward shift of the curve. A larger h value also increases the temperature difference between the two locations in the liquid; this increase is expected since the rate of heat accumulation at the location near the end of the liquid pool becomes higher than that at the location near the specimen/liquid gallium interface as h is increased (i.e. the rate of heat loss is lower for the location farther away from the boundary where heat is being extracted). The specimen thickness has a similar effect on the resulting temperature distribution as the specimen thermal diffusivity. As indicated by Figure 5.65, an increase in the spec-imen thickness reduces the amount of heat extraction from the liquid. The temperature change at the other side of a thicker specimen has less influence on the liquid tempera-ture. Decreasing the thickness results in the opposite effect. The behaviour exhibited in Figure 5.65 is similar to that in Figure 5.62. 5.4 M o d e l Ana lys i s of the Exper imen ta l Results The analysis was done by using the finite element model. As mentioned before, the model takes experimental input temperatures at one boundary to predict the corresponding temperature changes at other locations in the liquid gallium. The calculated temperature responses were adjusted to fit the measurements by altering the values of the parameters in the model. As was done in the sensitivity analysis, the locations chosen for predicting the experi-mental temperature-time curves are at 2 mm from the specimen/liquid gallium interface and 2 mm from the other end of the liquid at the vertical axis of the domain for the 154 thickness of specimen ( m ) 0.002 0.001 ( reference ) 0.0005 630 i i i i I I I L J L J L 200 400 600 800 Time ( seconds) 1000 1200 1400 Figure 5.65: Effects of the specimen thickness 155 locations of TC2 and T C I , respectively. Since the calculated temperature does not crit-ically depend on the location in the domain within a radius of 2 mm (a change of 2 mm results in a temperature change of less than 0.5 °C), the positions of TC2 and T C I in the experiments need not be known to a great precision; this fact alleviates the fitting process. 5.4.1 Corrections and Utilization of the Temperature Measurements Minor corrections were made to the experimental results, and a method was deviced to utilize them so that they could be analyzed by the mathematical model. First, the measurements were corrected for the thermocouple errors and furnace temperature gradi-ent. The time-varying input temperature profile required by one of the model boundary conditions was then obtained from the corrected measurements. 5.4.1.1 Corrections for the Thermocouple Errors and Furnace Temperature Gradient The thermocouples were calibrated against the melting point of tin. The readings were within 1.5 °C of the melting temperature of tin. The calibration corrections were applied to the measured temperatures. The model assumes a homogeneous temperature distribution in the liquid gallium. Experimentally, temperature gradients in the liquid were observed. In accounting for the furnace temperature gradient, the average temperature at the positions of T C I and TC2 was used as the starting temperature before cooling started. The temperature readings at time equal to 0 second were adjusted to the value of the average temperature; all subsequent readings were shifted by an amount equal to the differences between the average temperature and the initial measured temperatures. Hence, Figure 5.55 now looks like Figure 5.66 without the initial temperature differences. In so doing, the steady-156 Figure 5.66: Corrected thermocouple readings 157 state temperature gradient in the furnace was eliminated and the initial temperature distribution was homogenized. 5.4.1.2 Utilization of the Thermocouple Readings The four thermocouple readings were utilized in the following way. The TC4 and TC3 readings were used to obtain the temperature-time curve at the surface of the specimen. This was done by extrapolating these two temperature-time profiles to the specimen surface assuming a constant temperature gradient along the length of the liquid gallium. This surface time-varying temperature profile is necessary to provide one of the boundary conditions for the model. The model then calculates the corresponding temperature-time curves at the other two thermocouple positions, i.e. TC2 and T C I locations. Depending on the thermal diffusivity of the specimen used in the model, the calcu-lated TC2 and T C I curves vary. Hence, by matching the calculated curves with the experimental curves using a set of calibrated parameters, the model should give unique specimen thermal diffusivity from which the thermal conductivity can be calculated. A flowchart describing the procedure for the calculation is shown in Figure 5.67. 5.4.2 Calibrations of the Overall Heat-Transfer Coefficient In the model, there are 3 unknown parameters which must be determined in order to reproduce the experimental curves. They are the thermal diffusivity of the specimen as, the effective thermal diffusivity of the liquid gallium aaa which incorporates heat transfer by convection in the liquid, and the overall heat-transfer coefficient h. The data points [88] of as as a function of temperature for Pyrex and B N are shown in Figure 5.68 and Figure 5.69, respectively. Included in the figures are the best fit lines 158 Temperature Measurements : I thermocouples specimen materials Mathematical Model known for Pyrex and BN, but not for PBN crucible known for all specimens needs to be calibrated known for all specimens parameters calibrated overall heat-transfer coefficient fit TCI andTC2 temperature-time curves density, specific heat capacity thermal diffusivity of the PBN crucible thermal conductivity of the PBN crucible Figure 5.67: Reference flowchart for the procedure used to obtain the thermal conduc-tivity of the P B N crucible 159 X <D o CO o " 5§ JL X • H Q 5 -4 -T 1 1 1— • data point — fitted line — extended fitted line i 1 1 1 r i i i i _ - g - * - W M - « - r J 1 I I l ' i 200 400 600 Temperature (°C) 800 1000 Figure 5.68: Thermal diffusivity of Pyrex [88] used for model calculations in the fitting process (equation 5.49) 160 I i i I i i i i i i I 0 200 400 600 800 1000 Temperature Figure 5.69: Thermal diffusivity of B N [88] used for model calculations in the fitting process (equation 5.50) 161 to the data points, which show a large scatter. The fitted equations are apyre* = 6.9 x l f r 7 + 4.3 x 10~9 T - 2.4766 x 10" 1 1 T 2 + 8.9048 x 10~ 1 4 T 3 (5.49) and aBN = 2 x 10~5 - 5.01 x 10~8 T + 6.7458 x 1 0 - 1 1 T2 - 3.189 x 1 0 - 1 4 T 3 (5.50) where a is in m2/s and T is in °C. ocpyrex was obtained by directly curve-fitting the available thermal diffusivity data, while &BN was calculated from the fitted equations of kBN and CPBN given below kBN = 32.82 - 6.547 x 10" 3 T - 3.757 x 10" 7 T2 (5.51) CPBN = 688.13 + 2.688 T - 1.833 x 10" 3 T2 + 3.984 x 10~7 T 3 (5.52) and a density value of 2250 kg/m3 [88]. kBN and CPBN are in W/m.°C and J/kg.°C, respectively. The whole equation for a was incorporated into the model when fitting. To simplify the calibration process, the extent of liquid convection was assumed to be temperature-independent. This is a reasonable assumption since the temperature gradi-ent in the liquid was small for all the experiments. To account for the liquid convection, an arbitrarily chosen factor of 5 was used for OLGO- Figure 5.70 shows the 'stagnant' and the magnified (effective) a a a as a function of temperature [88]. The data points were fitted to obtain aoa = 1.03221 x 10~5 + 4.25 x 1 0 - 8 T + 4.457125 x 10" 1 2 T2 (5.53) With otGa which accounts for the liquid convection held at 5 times the 'stagnant' value, the only unknown parameter in the model is h. h is determined by fitting the model-calculated curves to the experimental curves of Pyrex and B N for which a„ is known. 162 30 T r ~i 1 1 1 1 r • data point — fitted line extended fitted line J L 200 400 600 Temperature C O -800 1000 Figure 5.70: 'Stagnant' and effective thermal diffusivity of liquid gallium [88] used for model calculations in the fitting process (equation 5.53) 163 5.4.2.1 Assumed Function of the Overall Heat-Transfer Coefficient The overall heat-transfer coefficient between the system and the surrounding, h, can be a complex function of many variables such as the temperature difference between the surface temperature and ambient temperature, the ambient temperature itself, the cooling or heating rate, the ambient gas, the gas pressure, etc. Here, h is simplified to be dependent upon the ambient temperature and the temperature difference between the domain boundary and surrounding only. h is assumed to follow an expression given by h = Clexp {C2 (T - Ta)} (5.54) where Ta is the ambient temperature in °C. The unit of h is W/m2 °C. C\ and C2 are constants which are assumed to vary with the ambient temperature. The values of these two constants are determined by fitting the experimental results. 5.4.2.2 Discussion on the Discrepancy between the Experimental and Cal-culated Curves A typical fit for the locations of T C I and TC2, taken from the fit for experiment no. 9 of Pyrex in Table 5.13, looks like Figure 5.71. The calculated curves agree with the experimental curves quite well during cooling. However, the calculated curves tend to be higher during heating. Such a discrepancy can be understood by examining the actual experimental system and the model. In the actual experimental system, the temperature at the location of TC3 reached that of TC2 in a short period of time after cooling. This indicates that the rate of heat input into location 3 is higher than that into location 2. The reason for this is that the rate of heat input from the ambient surrounding is theoretically higher for location 3 since it has a larger temperature difference with respect to the surrounding. 164 410 0 200 400 600 800 1000 1200 1400 Time (seconds) Figure 5.71: Typical fit between tlie experimental and calculated temperature-time curves (taken from experiment no. 9 of Pyrex) 165 In the model, the location at the surface of the specimen is at a fixed temperature as required by the boundary condition. Since the specimen surface temperature is fixed, the temperatures at the locations of TC2 and T C I have to reach that of the surface instead of the reverse which was observed experimentally. However, this condition is not possible for the model to satisfy unless the heat-transfer coefficient h is extremely small. For h values that fit the cooling sections of the curves, the calculated temperatures at the locations of TC2 and T C I during heating are higher than those of the experimental. It is possible to force the calculated curves to more accurately match with the heating sections of the experimental curves by adjusting the values of as and h as can be seen from the sensitivity analysis. In light of the many assumptions made in the model and the uncertainties of the thermal-physical properties used, an excellent match of the calculated and experimental curves can not be expected. The calculated curves shown in the figure can be considered as the best representation of the experimental curves under the simplifications and uncertainties involved. Consequently, this method is both relative and approximate. As the objective is to know the rate of heat flow across the specimen, it is much more important to match the cooling sections of the curves since the cooling process represents a uni-directional heat flow through the specimen better. This can be seen from the experimental curves which have a large temperature difference between the locations of TC3 and TC2 during cooling but only a small difference during heating. Thus, the effort is directed in fitting the cooling sections of the curves. 5.4.2.3 M o d e l Predic t ions A s s u m i n g Zero Overa l l Heat-Transfer Coefficient Calculations are performed to predict the experimental curves assuming that h is zero. This corresponds to the case of an insulated boundary. Using the time-varying input temperature profile taken from experiment 9 of Pyrex in Table 5.13, as shown in Figure 166 5.72, the predicted curves are obtained and compared with the experimental in the same figure. It can be seen that the predicted curves are lower than the experimental. Realis-tically, this is anticipated since there is no heat gain from the surrounding in the model calculations while there is heat gain in the actual experimental system. Similar predictions are observed for all the other curves. This observation indicates that the model is able to simulate the basic features of the experimental process within the limitations involved. It also shows that the experimental system was not perfectly insulated. Thus, h has a positive value which depends primarily on the extent of insula-tion. 5.4.2.4 F i t t i n g the Temperature-Time Curves of P y r e x The temperature for this set of experimental measurements ranges from 60 to 460 °C approximately. The time-varying input temperature profiles are shown in Figures 5.73a and 5.73b. The comparisons of the calculated curves with the experimental curves at location 1 (TCI) and 2 (TC2) are shown in Figures 5.74a and 5.74b. The corresponding C\ and C2 values obtained in the fitting process are shown in Table 5.17 for the different ambient temperatures. 5.4.2.5 Consis tency of C\ and C2 Constants Before fitting the experimental results of B N to find the values for C\ and C2, the values obtained from the experimental results of Pyrex were used to predict the experimental curves of B N at about 430 °C. The comparison is shown in Figure 5.75. It can be seen that the fit is reasonably well but the agreement is not as good as that for Pyrex. In order to get a better fit, the thermal conductivity of B N was reduced by half as shown in Figure 5.76. The curves calculated using the reduced ksN is included in Figure 5.75. Such a reduction of ksN necessary to fit the experimental results of B N suggests that 167 460 450 440 V-( £ 430 a H 420 410 — — i 1 1 1 1 1 1 1 1 1 1 1 1 r adiabatic boundary conditions experimental calculated input profile 400 J I I I I I I I I I I I L 0 40 80 120 160 Time (seconds) 200 240 280 Figure 5.72: Comparison of the model predictions with experimental curves assuming zero overall heat-transfer coefficient (insulated boundary) 168 Figure 5.73a: Time-varying input temperature profiles of Pyrex used for model calcula-tions (experiment no. 1 to experiment no. 5) 169 Time (seconds) Figure 5.73b: Time-varying input temperature profiles of Pyrex used for model calcula-tions (experiment no. 6 to experiment no. 9) 170 T 1 i 1 1 1 r ~ — i 1 1 1 1 1 r 0 40 80 120 160 200 240 280 Time (seconds) Figure 5.74a: Comparisons of the experimental and calculated curves at the locations of T C I and T C 2 for Pyrex (experiment no. 1 to experiment no. 5) 171 Figure 5.74b: Comparisons of the experimental and calculated curves at the locations of T C I and T C 2 for Pyrex (experiment no. 6 to experiment no. 9) 172 Figure 5.75: Comparisons of the experimental and calculated curves at the locations of T C I and T C 2 for B N using C\ and C2 values from Pyrex for normal kBN and reduced kfjN (0.5 x normal ksN) 173 Figure 5.76: Normal and reduced (0.5 x normal) thermal conductivity of B N [88] (equa-tion 5.51) 174 Ambient Temperature, Ta (°C) C i x 106 C2 x 104 65.9 65 140 102.8 70 145 142.0 80 150 184.1 90 155 251.8 100 160 303.1 200 170 356.8 300 180 404.5 400 190 458.3 500 200 Table 5.17: C\ and C2 values for Pyrex at different ambient temperatures Ambient Temperature, Ta (°C) C i x 104 C2 x 103 431.4 6 16 532.9 9 17 630.2 11 18 728.7 13 19 820.1 15 20 Table 5.18: C\ and C2 values for B N at different ambient temperatures there might be a 50 % error in the values of the derived thermal conductivity using this fitting method. But the large scatter in the reported thermal-conductivity values for B N (see Figure 5.76) might well be a major reason for this discrepancy. 5.4.2.6 Fitting the Temperature-Time Curves of B N The approximate temperature range is from 390 to 820 °C. Using the same fitting pro-cedure as was done for Pyrex, a set of values for C\ and C2 was obtained for the B N results employing the experimental input profiles shown in Figures 5.77a and 5.77b. The values of C\ and C2 are listed in Table 5.18. The calculated results are compared 175 Figure 5.77a: Time-varying input temperature profiles of B N used for model calculations (experiment no. 1 to experiment no. 3) Figure 5.77b: Time-varying input temperature profiles of B N used for model calculations (experiment no. 4 to experiment no. 5) 177 with the experimental curves in Figures 5.78a and 5.78b. 5.4.2.7 F i t t e d Equations for C\ and C2 The values of C\ and C2 for the two temperature ranges listed in Tables 5.17 (Pyrex) and 5.18 (BN) are combined and fitted to obtain two equations to be used in the expression for h. The points and the polynomial fitted lines are displayed in Figure 5.79 and Figure 5.80 as a function of ambient temperature for C\ and C2, respectively. The corresponding equations are Ci = 2.3656 x 1 0 - 4 - 2.76 x 1 0 - 6 Ta + 1.09 x 10~8 Ta2 - 6.9 x 1 0 - 1 2 Ta3 (5.55) and C2 = 1.43212 x 1 0 - 2 + 7.113 x 10" 6 Ta (5.56) The C\ and C2 values obtained for B N differ from those obtained for Pyrex. The discrepancy, particularly in the C2 values, can be attributed to the difference of insulation in the Pyrex and B N experiments. It was difficult to maintain the same insulation from one set of experiments to another. In fitting the temperature-time curves of the P B N crucible, the fitted equations of C\ and C2 are assumed to be applicable since this set of experiments presumably has an average insulation. 5.4.3 T h e r m a l Diffusivi ty F i t t i n g of the Tempera ture-Time Curves of the P B N Cruc ib l e The temperature range covers approximately from 110 to 820 °C. For this set of ex-perimental results, all the parameters involved in the model were known except for the thermal diffusivity of the P B N crucible in the thickness direction, OCPBN- A S before, ctGa was assumed to be 5 times the 'stagnant' value. Equation 5.54 along with equations 5.55 178 Figure 5.78a: Comparisons of the experimental and calculated curves at the locations of TCI and TC2 for BN (experiment no. 1 to experiment no. 3) 179 830 820 810 800 t-790 780 770 760 i 1 1 1 1 1 1 1 1 1 1 1 1 r experimental calculated 0 40 80 120 160 Time (seconds) 200 240 280 Figure 5.78b: Comparisons of the experimental and calculated curves at the locations of TCI and TC2 for BN (experiment no. 4 to experiment no. 5) 180 0 I 1 I 1 I I I 1 l _ 0 200 400 600 800 Ambient Temperature, Ta (°C) Figure 5.79: C\ as a function of ambient temperature (combined Pyrex and B N values) (equation 5.55) 181 Figure 5.80: C2 as a function of ambient temperature (combined Pyrex and B N values) (equation 5.56) 182 and 5.56 was employed to provide the overall heat-transfer coefficients in determining &PBN for this set of experimental curves. The fitting procedure carried out here is anal-ogous to the procedure used in fitting the Pyrex and B N curves. The difference is that as was adjusted to match the calculated and experimental curves in this case instead of h in the other two cases (Pyrex and BN) . Using the time-varying input profiles shown in Figures 5.81a, 5.81b, 5.81c and 5.81d, an equation for QPBN was derived by trial and error until a reasonable fit between all the calculated curves and experimental curves was attained. The agreement of the fit can be seen in Figures 5.82a, 5.82b, 5.82c and 5.82d. The resulting OIPBN is presented in Figure 5.83 as a function of temperature. Two equations are used to describe the line shown in Figure 5.83 for two temperature ranges as given by aPBN = 2.52 x 1(T 7 + 6 x l f r 1 0 T + 2.8 x l f r 1 2 T 2 + 1.18 x 10~ 1 4 T 3 (5.57) 110 < T < 677 °C and ctPBN = -2.539 x l O - 6 + 2.22 x 10~8 T - 1.80 x 1 0 - 1 1 T 2 + 4.6 x l O - 1 5 T 3 (5.58) 677 < T < 820 °C T is in °C, and OCPBN is in m2/s. The temperature range in the figure has been extended to cover from 0 °C to 1000 °C. The decrease in the slope of the thermal diffusivity curve above 675 °C may be an artifact of the fitting process. However, among many other shapes of the curve, the curve presented in Figure 5.83 gives the best overall agreement to the measurements of Figure 5.82a to 5.82d. 183 Figure 5.81a: Time-varying input temperature profiles of the P B N crucible used for model calculations (experiment no. 1 to experiment no. 4) 184 450 i \ 1 1 1 1 1 1 i 1 1 1 1 i r Time (seconds) Figure 5.81b: Time-varying input temperature profiles of the PBN crucible used for model calculations (experiment no. 5 to experiment no. 8) 185 Figure 5.81c: Time-varying input temperature profiles of the P B N crucible used for model calculations (experiment no. 9 to experiment no. 12) 186 Figure 5.81d: Time-varying input temperature profiles of the P B N crucible used for model calculations (experiment no. 13 to experiment no. 16) 187 Figure 5.82a: Comparisons of the experimental and calculated curves at the locations of T C I and TC2 for the P B N crucible (experiment no. 1 to experiment no. 4) 188 ( seconds ) i >oA rurves at the locations . of the experimental and ^ ^ e r i m e n t no. 8) r ROW Comparisons of tae exy . t n o . 5 to exper T C 2 fo the P B N « u a U e 189 Figure 5.82c: Comparisons of the experimental and calculated curves at the locations of T C I and TC2 for the P B N crucible (experiment no. 9 to experiment no. 12) 191 Figure 5.83: Thermal diffusivity of the P B N crucible in the thickness direction obtained from the fitting process (equations 5.57 and 5.58) 192 5.4.4 Sensitivity of OLPBN, &Ga, h and t After obtaining CCPBN, a sensitivity analysis was done to see the effects of the parameters used in the fitting process. Unlike the sensitivity analysis done before in which the fitting procedure used constant values for all the parameters, the parameters in this case are in the form of equations as a function of temperature which were determined during fitting except for t. Accordingly, the fitted equations were used as the reference of comparison. The input profile is the same as the one shown in Figure 5.61. Each parameter is multiplied by 2 and 0.5. The calculated curves are compared with the reference curves. The comparisons are presented in Figures 5.84, 5.85, 5.86 and 5.87. Similar effects are observed as the sensitivity results using constant values. Since the parameters characterize the heat transfer of the experimental system, such a sensitivity analysis may be utilized to assess the error involved in the fitting process. A significant shift of the curves can be seen in the figures as a result of multiplying or dividing the parameters by 2. The factors of 2 and 0.5 seem to provide an upper bound and a lower bound of the fitting error, respectively. 5.4.5 Sources of Errors in the Derived Thermal Diffusivity There are three major sources of errors in deriving the thermal diffusivity using this method. The first one is from the experimental measurements. The second contribution of errors comes from the mathematical model. The last source of errors lies in the analysis of the experimental results by the mathematical model. The main experimental measurements are temperatures. The temperature measure-ments contain errors due to the inaccuracy of the thermocouples and the data acquisition system. Even though the thermocouples were calibrated against the melting point of tin, the errors in the measurements would still exist, and they were expected to be larger 193 Figure 5.84: Effects of a P B N 194 730 I 1 1 1 1 1 \ 1 1 1 1 1 1 1 r 6 5 0 I 1 1 1 1 1 1 1 1 1 1 1 1 1 1— 0 200 400 600 800 1000 1200 1400 Time (seconds) Figure 5.85: Effects of ctGa 195 Figure 5.86: Effects of h 196 Figure 5.87: Effects of t 197 at higher temperatures. The data acquisition system used to record the temperature readings was only accurate to within 0.5 °C owing to the electrical noises. The mathematical model is anticipated to contribute considerable errors. The simpli-fications and assumptions made in the formulation limit the predictive capability of the model. As a consequence, errors arise from the difference in the resulting temperatures between the complex experimental system and the simplified model. Moreover, the nu-merical technique (finite element method) itself is an approximation, and hence there is an additional numerical error. The analysis of the experimental results by the model probably constitutes the most errors. The temperature measurements were corrected so that they could be utilized by the mathematical model. This procedure introduced errors in the derived thermal property. The calibrations of the overall heat-transfer coefficient employed an assumed function of the parameter whose constants were obtained by fitting the results from the reference specimens (Pyrex and BN) with known thermal conductivities. The calibrated expression was assumed to apply to the results from the unknown specimen (the P B N crucible) which might have a different extent of thermal insulation. Errors in this process could come from the assumptions involved in the overall heat-transfer coefficient and the inaccuracy of the thermal conductivities of the reference specimens. The most apparent errors, however, He in the process of matching the model-calculated results with the measurements in the sense that it was subjective in choosing the 'best fit' criterion. A quantitative assessment of the percent error in the derived thermal diffusivity is difficult to do given the above sources of errors. A proper analysis of the errors is beyond the scope of this work. 198 5.5 T h e r m a l Conduc t i v i t y of the P B N Cruc ib le To convert OLPBN to the thermal conductivity, kpBN, a knowledge of the density and specific heat capacity of the P B N crucible must be available. The density was taken to be constant having a value of 2185 kg/m3 for P B N [2] since this information was not available as a function of temperature. The specific heat capacity was obtained by fitting the data points available for P B N [2] as given in Figure 5.88. The corresponding fitted equation for the specific heat capacity is given by CPPBN = 7 6 7 - 7 7 + 2.99 T - 3 x 1(T 3 T 2 + 1.22 x 1(T 6 T3 (5.59) where C P P B N is in J/kg.°C. Utilizing the derived ctpBN, kpsN was calculated using kpBN = CtPBN C P P B N PPBN (5.60) The calculated kpBN is presented in Figure 5.89 as a function of temperature. As a result of describing OLPBN in two equations, kpsN is given by two equations in the two temperature ranges. The corresponding equations for the line shown in the figure are kpBN = 4.132 x 10 _ 1 + 3.1 x 10~3 T + 1.81 x 10" 6 T 2 + 6 x 10~8 T 3 (5.61) 110 < T < 677 °C and kpBN = -19.35 + 1.116 x 10 _ 1 T - 9.28 x 10" 5 T 2 + 2.76 x 1 0 - 8 T3 (5.62) 677 <T < 820 °C where T and kpsN are in °C and W/m.°C, respectively. The line has been extended to cover from 0 °C to 1000 °C. The thermal-conductivity curve shown in Figure 5.89 appears to have a maximum above 1000 °C. At high temperatures, phonon scattering is known to limit the thermal 199 0 • data point fitted line 200 400 600 Temperature 800 1000 Figure 5.88: Specific heat capacity of P B N [2] (equation 5.59) 200 0 200 400 600 800 1000 Temperature Figure 5.89: Thermal conductivity of the P B N crucible in the thickness direction calcu-lated using ocpBNt PPBN and C P P B N (equations 5.61 and 5.62) 201 conductivity of electrically insulating solids. This might give rise to the characteristic maximum of the material as a function of temperature [25]. 5.6 Compar i son of the T h e r m a l Conduct iv i t ies The comparison between the thermal-conductivity values of the P B N crucible in the thickness direction (predominantly "c" direction) obtained in this work and the published values for P B N in the "a" direction and "c" direction is shown in Table 5.19. The values listed in the table show that the thermal conductivity of the crucible material is appreciably lower than the published values for P B N in both the "a" and "c" directions at low temperatures. At high temperatures, the crucible values are between the published values for the "a" and "c" directions. The thermal conductivities of Pyrex, B N , P B N in the "a" and "c" directions, and the P B N crucible in the thickness direction are shown in Figure 5.90 for comparison. The thermal conductivities of Pyrex, P B N in the "a" and "c" directions are calculated from the following equations: kpyrex = 9.36 x 10 _ 1 + 1.005 x 1 0 - 3 T + 6.2 x 10~8 T2 (5.63) kpBN..a„ = 101.1 - 4.264 x 10~2 T (5.64) kpBN..c„ = 1-654 + 8.57 x 10~4 T (5.65) The equations were obtained by curve-fitting the data available in the literature [88,2]. The thermal conductivities of Pyrex and P B N in the "c" direction are the lowest, followed by the P B N crucible in the thickness direction, B N , and P B N in the "a" direction. Since the P B N crucible through thickness is predominantly parallel to the "c" direc-tion, it is expected to have a thermal conductivity similar to and slightly higher than that of P B N in the "c" direction and lower than that in the "a" direction. However, the shape 202 Temperature P B N T H E R M A L C O N D U C T I V I T Y ( W/m.°C ) (°C) "a" Direction "c" Direction the Crucible in the Thickness Direction ( Predominantly "c" Direction ) Reference [2] Present Results 25 500 1000 104.5 71.6 62.7 1.672 2.09 2.51 0.49 9.92 27.05 Reference [25] Present Results 20 800 63 63 1.5 2.9 0.48 24.67 Table 5.19: Comparison of the thermal conductivity values of P B N in the "a" and "c" directions with the values obtained in this investigation for the P B N crucible in the thickness direction 203 Figure 5.90: Comparison of the thermal conductivities of Pyrex [88], B N [88], P B N in the "a" direction [2], P B N in the "c" direction [2], and the P B N crucible in the thickness direction (this work) (equations 5.63, 5.51, 5.64, 5.65, and 5.61 and 5.62, respectively) 204 of the line is different from those of the "c" and "a" directions. At low temperatures, the P B N crucible curve coincides with the "c" direction curve but increases much more rapidly than the "c" direction curve with increasing temperature. As anticipated, the line never reaches that of the "a" direction. The rapid increase in the thermal conductivity of the crucible material could be attributed to several factors. First, since the material consists of layers and has a positive thermal expansion coefficient, the gaps between the layers might close up due to thermal expansion as the temperature is increased which gives a better contact between the layers. This should result in an increase in the thermal conductivity of the material. The second factor could be associated with the crucible being used a number of times, the thermal conductivity changing with use. A third factor could be related to the method used to derive the thermal conductivity, in which a 'best fit' criterion was arbitrarily selected to fit the data. The values of the thermal conductivity obtained are also dependent on the assumptions and simplifications of the model, and the uncertainties of the values of the physical parameters used. Nevertheless, as mentioned earlier, The method used in the measurements is not standard. Even with standard methods, different values are often obtained for the same material by different investigators especially at high temperatures, indicating that the standard methods are also in need of further improvement to increase the reliability and accuracy of results. 5.7 Relative Thermal Resistance Analysis of the Mathematical Model A n analysis of the thermal resistances in the mathematical model was carried out to see the relative extent of heat transfer. Figure 5.91 shows the thermal resistance circuit considered in the model domain. There are two convective/radiative resistances, one along the side surface (Rni) and the other along the end surface of the liquid gallium 205 convection/ radiation to domain T C I conduction through liquid gallium TC2 conduction through specimen 1 convection/ —\f\f\/ radiation to domain Figure 5.91: Thermal resistance circuit in the model domain used for the determination of the thermal conductivity of the PBN crucible 206 (Rh2)- There is a conductive resistance across the liquid gallium to the location of T C I ( R k G a i ) , from the location of T C I to that of TC2 (RkGa2 )> a n d fr°m t n e location of TC2 to the liquid gallium/specimen interface (RkGai)- The last is the conductive resistance from the interface to the surface of the specimen (Rk,). The resistances are given as Rin Rh2 RkGai RkGa2 Rk. Ai and A2 are the areas along the side and end boundary, respectively. A is the area across the width of the domain. L\ is the distance between the end boundary of the liquid gallium and location 1, which is the same as the distance between location 2 and the interface. L2 is the distance between location 1 and location 2. L is the thickness of the specimen. The overall heat-transfer coefficient, h, and the thermal conductivities of the liquid gallium and specimen, kca and ks, are taken to be at 500 °C. The values are h = 7.2x 10~4 W/m2.°C kGa = 350 W/m.°C kpyrex - 1.5 W/m.°C kBN = 29.5 W/m.°C kpBN = 9.9 W/m.°C 1 hAt (5.66) 1 hA2 (5.67) U kaa A (5.68) L2 kGa A (5.69) L kg A (5.70) 207 Using these values, the resistances are calculated below: ^ O . O O O T ^ x 0.005= 2 7 7 7 7 8 ' ° / W * • = IMOn ' x 0.029 = 4 7 8 9 3 ° ° ' W Rka = ^00? = 0.001143 °C/W G a i 350 x 0.005 ' 0 024 = 35070005 = 0 0 1 3 7 1 4 ° ° I W OT0305 = 0 1 3 3 3 3 3 ^ / W for Pyrex ^ = 29.5°x0001005 = 0 0 0 6 7 8 0 ° C < W f O T B N Rk, = g ^ Q Q Q g = 0.020202 °CjW for the P B N crucible The thermal circuit as well as the resistance values is shown in Figure 5.92. The total resistances are heat gain from surrounding: Rh = i ? / > 1 ^ = 40850 °C/W Rh2 + Rhr heat conducting through liquid gallium: RkGa = RkGai + RkGa2 + ^fc G o i = 0.016 °CjW heat conducting through specimen: Rk, = 0.133333 °C/W for Pyrex or Rkt = 0.006780 °C/W for B N or Rk, = 0.020202 °C/W for the P B N crucible It turns out that the convective/radiative heat gain from the surrounding has the highest resistance which may imply that the experimental system was well insulated. The resistances across the liquid gallium and specimen are several orders of magnitude lower than that of the convective/radiative mode. Among the materials used, Pyrex is the most resistant to conduction while B N is the least. Depending on the specimen material, the liquid gallium may have a smaller or a larger thermal resistance compared to the specimen. 208 convection/ radiation to domain conduction through liquid gallium conduction through specimen 277778 °C/W 47893 °C/W 2 E 0.001143 °CJW :0.013714 °C/W aEO.001143 °C/W ,0.133333 °CJW for Pyrex 10.006780 °C/W for BN 0.020202 °C/W for the PBN crucible Figure 5.92: Thermal circuit with relative resistances of the mathematical model used for the determination of the thermal conductivity of the PBN crucible Chapte r 6 S u m m a r y and Conclusions The objective of this research was to characterize the materials relevant to the Czochral-ski crystal growth system, in particular the thermal conductivities of the melt and the crucible. Experimental measurements and mathematical modelling were used to examine the effect of an applied magnetic field of 0.099 T on the heat transfer in a germanium melt and to determine the thermal conductivity of a pyrolytic boron nitride (PBN) crucible in the thickness direction. The summary and conclusions of this work are listed below: 1. Applying a vertical magnetic field to a germanium melt in the C Z growth config-uration increases radial temperature gradients in the melt, due to suppression of radial fluid flow. 2. The measured radial temperature gradients in the melt were fitted to a steady-state conduction-dominated mathematical model of heat transfer, by adjusting the effec-tive thermal conductivity of the melt. Both thermal insulating and non-insulating boundary conditions were considered at the bottom surface of the melt. 3. It was found that the effective thermal conductivity of the melt, in the presence of a vertical magnetic field of 0.099 T, was lower than the effective melt thermal conductivity without a field by a factor of seven. This was the case for both the insulating and non-insulating boundary conditions at the bottom melt surface. This result shows that fluid flow in the melt increases the radial heat transfer in the melt, by at least a factor of seven. 209 210 4. The thermal diffusivity through the wall of a pyrolytic boron nitride crucible, used in growing GaAs crystals, was obtained in the temperature range of 0 to 1000 °C by matching the results calculated by a transient heat-transfer mathematical model with the measured thermal responses. The thermal conductivity through the wall of the crucible, calculated from the thermal diffusivity, was found to be strongly temperature-dependent. The thermal conductivity was 0.48 Wjm.°C at 20 °C, increasing at an accelerating rate to approximately 23 W/m.°C at 700 °C, then slowly rising nearly linearly to 27.05 Wjm.°C at 1000 °C. The values at high tem-peratures (above 300 °C) differ markedly from the P B N values in the "a" direction and to a less extent in the "c" direction. However, these values are between those of P B N in the "a" direction and "c" direction, consistent with the P B N crucible in the thickness direction being predominantly parallel to the "c" direction. 5. The published values of the thermal conductivity of P B N in the "c" direction are higher than the present measured values at 20 °C. Wi th increasing temperatures the published values drop below the measured values, becoming appreciably lower above 500 °C. The measured low values at the low temperatures are probably attributed to delamination of the P B N crucible material with use. It thus appears likely that the thermal conductivity across a given crucible wall decreases with the number of crystals grown using that crucible. Chapter 7 Recommendations for Future Work The temperature measurements in the CZ germanium melt contained a large scatter band, and they were obtained for only a few melt depths. Hence, more measurements should be made at more locations in the melt to better characterize the temperature distributions and gradients. Once the temperature fields in the melt are clearly estab-lished, a conduction-dominated heat-transfer mathematical model which includes more complex boundary conditions and temperature-dependent heat-transport properties can be developed to analyze the influence of magnetic field on the effective thermal conduc-tivity of the melt. Different magnetic field strengths should be applied to determine the dependence of the reduction factor on the magnetic field strength. Similar investigation can also be carried out using different semiconductor melts. The method used in the determination of the thermal conductivity of the pyrolytic boron nitride crucible should be considered as a first attempt in generating the ther-mal conductivity of the material. Improvement of the experimental procedure and the mathematical model to obtain more reliable results can be achieved by conducting more experiments on standard specimens. 211 Appendix A Derivations of the Finite Element Equations A . l Cylindrical Coordinates (r, z) The content of this section is a modification of the derivations done in reference [15]. Similar derivations can also be found in many standard books on the finite element method. The Galerkin method is used to derive the elemental matrices and vectors for equation In order to solve equation A.71 using the Galerkin method, it is necessary that the dependent variable T within an element be piecewise approximated using trial or shape functions Nj by T « « £ Tj Nj (r,z) (A.72) where Tj is the temperature at node j, and n is the number of elemental nodes. Substi-tuting equation A.72 into equation A.71 will only give an approximation to the equality of equation A.71, resulting in an error, known as residual R, given by The weighted averages of the residual i i ! is then set to zero using the shape functions as the weighting functions. Mathematically, R is made orthogonal to the weighting functions: NiRrdrdz = 0 for t = l , 2 , . . . , n (A.74) 212 213 Substituting equation A.73 into equation A.74 with a slight rearrangement of the r term and the deletion of the summation sign gives 1.1 * {* (4 - 4 (t)) f} * - - <-) where i = 1,2,... , n indicating the row number, and j = 1,2,... , n indicating the column number in the matrix. Equation A.75 is integrated by parts to reduce the second-order derivatives to first-order derivatives. In general the integration by parts in two dimensions is given by where T is the boundary surrounding the domain of interest in (£, rj) coordinates, and n„ is the rj component of the unit normal vector to the boundary n. Equation A.75 after integration by parts becomes r dN- r dN-- / r Ni k Tj nrdT — / r Ni k T , n z oT = 0 (A.76) Jr dr Jr dz ' The boundary terms in equation A.76 ~ k T j ~&fnr ~ kTj^dfnz = ~ k T j i t give the Fourier law of heat conduction and represent the natural boundary condition which can be either specified heat flow at the boundary, or convective heat exchange with ambient surrounding. Only convective heat transfer is considered here. Applying the convective boundary condition dN--kTJ^- = h(TjNj-Tco) 214 with a constant heat-transfer coefficient h and ambient temperature results in I.Lr{lTi( dNi dNj dNi 6NA _ A r A 7 dTA , . dr dr dz dz j dt + J^rhNi (Tj Nj — T^) dT = 0 (A.77) Radiative boundary condition can be included by linearizing the fourth power of tem-perature and introducing a constant radiative heat-transfer coefficient hT given by hr = Fve (T2 +Tl>)(T + T00) where F is the radiation view factor, a is the Stefan-Boltzman constant, and e is the emissivity of the boundary surface. hr can be directly added to h so that h becomes the sum of the convective heat-transfer coefficient hc and the radiative heat-transfer coefficient hr. Equation A.77 is transformed into the form of matrix as [l<M] {T} + [cM] {T} = {/<e>} (A.78) where |iir(e)j is the thermal stiffness matrix, [CM] is the capacitance matrix, { /^} is the thermal (force) load vector, {T} is the nodal temperature vector, and | r | is the vector containing time derivatives of nodal temperatures, of the element. [l<M] and {/<«)} can be written as [KM] = [K^] + [4e)] + [KM] and {/(e)} = {/ie)} + {/r(e)} where [KM] is the conduction matrix, [ifj^] is the convection matrix, [.K^] is the linearized radiation matrix, {/i^ } is the convective load (force) vector, {/,le }^ is the linearized radiative load (force) vector, of the element. 215 Equation A.78 provides n simultaneous equations for the elemental nodal tempera-tures. The elemental matrices and vectors are ^ ^ / J r ^ S ^ + f f f l drdz r " r 3 (A.79) C\f =JzJrrpCpNiNjdrdz fl? =IrrhcNiT00ldT ft? =JrrhrNiT002 aT The next step is to choose the type of element and shape functions. Numerical integration must be employed to evaluate the integrals. For transient heat transfer, time-stepping scheme is applied to approximate the transient response. A.1.1 Type of Element and Shape Functions In order to ascertain finite element convergence as the element size is decreased, the type of shape function used must satisfy two criteria: continuity and completeness. To solve a differential equation that has the highest order of derivative equal to rn, a shape function continuous up to the (m — 1) derivative is required; in the finite element jargon, the shape function must satisfy C m _ 1 continuity. Completeness requires that the polynomial used as the shape function be a complete polynomial to the order m. Equation A.77 has a first-order derivative as the highest order of derivative, i.e. m — 1. Thus the shape function used must have at minimum C° continuity and be a complete polynomial up to the order 1. Therefore a linear shape function could be used in solving the differential equation. The type of element used here is an 8-noded isoparametric quadratic element, as shown in Figure A.93a. Correspondingly, 8 isoparametric quadratic shape functions are 216 Figure A.93: Element type: (a) 8-noded quadratic isoparametric element; (b) Quadratic isoparametric boundary element 217 used to approximate the temperature at any point in the element by 8 L i=l The shape functions are a function of £ and r) coordinates given by ^^TiNi^rj) (A.80) (A.81) J V i =-(l-0(l-'?)(l + e + ' / ) /4 N2 = -(l+t)(l-v)(l-t + r,)/4 J V 3 = -(l+0(l+'/)(l-e-'?)/4 N4 = -(l-Q(l+r,)(l+t-r,)/4: N5 = (l-e)(l-v)/2 N6 = (l-v2)(l+0/2 N7 = (l-e)(l+v)/2 J V 8 = ( l- i ? 2 ) ( l -0 /2 where —1 < £, T/ < 1. Boundaries of the domain are accounted for using a 3-noded one-dimensional quadratic element with its associated shape functions given by « i = - i ) ( l - i ) ) /2 * - , ( ! + , ) /2 ( A 8 2 ) J V 3 = 1 - »72 The one-dimensional element is shown in Figure A.93b. A.1.2 Mapping Functions and Numerical Integration The domain is in terms of r and z coordinates and the element shape functions are in terms of £ and n coordinates. This requires that the element in the (r, z) coordinates be mapped into the (£, vj) coordinates for numerical integration. The mapping function 218 must be quadratic since the curved boundaries of the element in the r-z plane need three points for their unique locations. The mapping function also must equal unity at its node and zero at all other nodes. Since the shape function satisfies all these characteristics, it is used as the mapping function, and such an element is termed isoparametric. Figure A.94 illustrates the coordinate and geometrical transformation of the element where r = Y.rtN,(i,V) * = X > i V i ( & i ? ) 1=1 t=i The numerical integration method used is the Gaussian quadrature which involves evaluating an integrand at a number of sampling points. Each evaluation is multiplied by a weighting value and the summation over all the sampling points is equal to the integral. The Gaussian-quadrature formula in one dimension is / f(0 dC = Wi + W 2 / (&) + . . . + Wmf(U) = £ Wiftfi) or in two-dimensional form, where / is a polynomial, and m\ and m2 are the number of sampling points in the £ and •q direction respectively. The locations of the sampling points and the corresponding weighting coefficients are given in Table A.20. Higher accuracy is obtained by employing more sampling points. If a polynomial expression is being numerically integrated in one dimension, it can be shown that for n sampling points In unknowns (/,• and £,) can be solved, and thus a polynomial of degree (2ra —1) could be constructed and exactly integrated. The corresponding error is of the order O (S2") where 8 is the point spacing [57]. Thus for a quadratic shape function the degree of the polynomial is 2. However during the formulation of the problem the shape function is multiplied by itself (Ni Nj), and therefore the function to be integrated 219 Figure A.94: Mapping of the quadratic isoparametric element from the (r, z) coordinates into the (£, r/) coordinates: (a) domain element; (b) boundary element 220 Sampling Location Number of Sampling Points Weighting Coefficients - 0.5773502691 + 0.5773502691 m = 2 1.0000000000 1.0000000000 - 0.7745966692 0.0000000000 + 0.7745966692 m = 3 0.5555555555 0.8888888888 0.5555555555 Table A.20: Sampling locations and weighting coefficients of the Gaussian-quadrature formula is an incomplete 4 t / l-order polynomial. The optimum number of sampling points is 3 (i.e. 2n — 1 = 4, .-. n = = 2.5 < 3). For two dimensions the number of integration points is 3 in the £ direction and 3 in the 77 direction, and thus a total of 9 sampling points is used (32 = 9). The larger the number of sampling points is, the longer the solving times will be. A reduced number of sampling points can be used and still achieve the accuracy of the larger sampling point integration. Using fewer than the optimum number of sampling points is known as reduced integration [58]. Reduced integration using 4 sampling points gives comparable results to nine-point integration if stability criterion is satisfied. Stability criterion for reduced integration requires that more than one element be present in the domain [57]. Nine sampling points are preferred when there are temperature-dependent thermo-physical properties and heat-transfer coefficients. More sampling points allow the variation in thermal conductivity, density, specific heat capacity, etc. with temperature to be better modelled over the domain. For a domain with constant thermo-physical properties, the reduced integration is often used to achieve shorter solving times. 221 A.1.2.1 Conduction Matrix The elemental conduction matrix dNi dNj dNi dN, + dr dz dz dz j must be numerically integrated. Equation A.83 can be written in matrix form as (A.83) = J J^rk [B]T [B] drdz where [B] a/Vi dN-dr dr dN dN-dz dz dN„ dr dN„ dz (A.84) (A.85) and [B]T is the transpose of [B]. Since the shape functions are a function of £ and 77, it is necessary to express and dr dz in terms of £ and 77. The limits of the integration is changed to -1 and +1 in the £-77 space accordingly. Using the chain rule of differentiation, the derivatives of the shape functions with respect to £ and 77 are given as dNi__dNi_dr_ dNj dz d£ dr d£+ dz d£ dNj__dNi_dr_ dNj dz dr) dr dr] dz dr) or in matrix form dNj dr dz f -v. dNj dNj < * > = dZ ae < dr - = [ J] < dN; dr dz dNj dNj . dn . L dn dr, J dz . I dz ) (A.86) where [J] is the Jacobian matrix defined as 222 The Jacobian matrix must be inverted in order to find the first derivatives of the shape functions with respect to r and z. After can be expressed in terms of £ and r) as ft r inverting the Jacobian matrix, the and dNj f ~\ dNj "ill h2 ' r \ dNj dt y < * dNj \ dz J dNj dr] . hi I22 dNj (A.88) where i.y is the content of the inverse of the Jacobian matrix. 7,-j's are given by i n = i l 2 = " hl = -I22 — dNj det[J] jr( ~' dr] E * dNi det[J] ^ d£ dNi (A.89) det [J] 1 E n n 5 > drj dNi det[J] £ •* d£ where det [J] is the determinant of the Jacobian matrix defined as (A.90) \i=l / \l=l / \l=l The expressions for the and are then substituted into the [B] matrix giving it in terms of £ and rj: [B] Jn-ae + Ji2-a^-r , 7- dNi hi -o£- + h2-af J11 -gj + JI2-QJ- ••• hi-px -r h T 9N2 1 T 9N7 hi-df + h2-af 11 sc -»21 - § £ + -»22 —" 12 dn Correspondingly, dr dz must be converted to the £-77 coordinates using dr dz = det [J] d£ dt] (A.91) Transforming the elemental conduction matrix from the (r, z) to (£, rj) coordinates changes equation A.84 to [tfW] = J1 rk [B]T [B] det[J) dtdr) (A.92) 223 Equation A.92 is numerically integrated using the Gaussian quadrature resulting in Wi and Wj = weighting constants £ and T) — coordinates at sampling points 3 sampling points are used in the £ direction, and 3 in the r] direction (m = 3). Thus, a total of 9 sampling points is employed in the element for this case. A.1.2.2 Convective/Radiative Matrices and Vectors Numerical integration of the following integrals m m rp [Kie)} = E E Wi Wi r^) k [ B ( M ] [%,„)] MA*,*)) (A-93) i=l j=lwhere m number of sampling points in one coordinate direction is necessary in order to evaluate the convective/radiative matrices and vectors of the element. Writing them in matrix form gives {fie)}= I rhc {N} dT where {N}T is the transpose of {N}. 224 In this case, numerical integration requires that dT be put in terms of r\. Before doing so, dT is defined in terms of r and z using the Pythagoras' theorem as dT - Vdr2 + dz2 Differentiating dT with respect to r\ yields dT ( dr dr) ^\drj) {dr}. Mapping into the rj space gives dT ' 3 dNiY / » dNi drj (A.94) The term on the right-hand side of equation A.94 can be taken as the boundary Jacobian, A, rewritten as A = 2MB dr\ Therefore dT = A drj The integrals to be carried out become [Kne)] =Jrrhc {N} {N}T Adrj {fie)} = jrrhc {N} T^Adrj [l<ie)] = f r h r {N} {N}T Adr] {/r(e)} = Jr rhr {N} T^Adr, or in terms of the Gaussian quadrature for numerical integration, = E W rM hc {#(,)} {NM} A(ni) (A.95) 225 m {fie)} = E Wi r ( „ ) /zc A(,.) (A.96) t=l [ffW] = E W i K {NM} {NM} A M ( A . 9 7 ) m {/^ } = E Wi TM hr {NM} T o o , A ( r , ) ( A . 9 8 ) •=1 A.1.2 .3 Capaci tance M a t r i x The elemental capacitance matrix integral to be numerically integrated is C\f = j J rpCpNiNjdrdz In matrix form, [C<e)] = J J rpCp {N} {N}T drdz The shape functions used for the elemental capacitance matrix are the same as those used for the elemental conduction matrix. The Jacobian matrix, [J], and its determinant are given by equation A.87 and A.90, respectively. Mapping into the £-r? space gives [ C ^ ] = J1 J1 rpCp {N} {N}T det[J] d£dV In the Gaussian quadrature form for numerical integration, m m rj-i [C ( e ) ] = E E ^ W3 rUiitlj)PCp del [j(e,ils)] (A.99) A .1 .3 T ime-S tepp ing Techniques For steady-state heat transfer, the equation to be solved [K] m = {/} (A.100) 226 is algebraic in nature, and hence can be evaluated directly using matrix algebra. However, for transient heat-transfer problems, the following equation must be solved: [K] {T} + [C] {T} = {/} (A.101) Equation A . 101 involves the time derivatives of temperature and is generally difficult to solve analytically. To overcome this problem, a recurrence technique is employed in which trial functions are used to discretize the time domain and the solution at each time step is obtained using the temperatures from the previous time step(s). The number of trial functions used affects the accuracy and stability of the solution. In general the greater is the number of shape (trial) functions used the better will the accuracy and stability be [12]. The recurrence method which uses two shape functions is known as a two-point (linear) recurrence scheme, while a three-point (quadratic) scheme uses three shape functions. The two schemes are often used in the following way. The two-point scheme which only requires the temperatures from the previous time step is used to obtain the solution at the very first time step. The three-point scheme which requires the temperatures from the two previous time steps is then used for all subsequent calculations. In the time-stepping method, temperature is approximated by {T} « £ Nt {Tt} (A.102) i where {T}} is the temperature vector at time /. The shape functions, Ni, are a function of time only. A.1.3.1 Two-Point Recurrence Scheme The shape functions are linear as shown in Figure A.95. The dimensionless time variable 227 Figure A.95: Linear shape functions (two-point recurrence formulae) for the time domain [57] 228 £ and the shape functions are o < C < i C = £ N, = l-C; Ni = -£i (A-103) N,+i = C; Nl+1 = ± Thus, {J} « {Ti} N, + {Tl+1} Nl+1 (A.104) where the subscripts I and / + 1 denote current and next time step, respectively. Substituting equation A.104 into equation A.101 and applying the weighted residual method yields Jl W {[K] ({2}} Ni + {Tl+1} Nl+l) + [C] ({T,} Nt + {Tl+1} 7Y J + 1 ) - {/}} d{ = 0 (A.105) for the time domain where W is the weighting function. Collecting the Tn and Tn+i terms in equation A.105 gives [K] £ WNtdC + [C] f W N, dc] {Ti} + [K] £ WNI+1 d( + [C] £ WN,+i dc] {Tl+i} - £ W {/} d( = 0 (A.106) Substituting the shape functions and their time derivatives from equation A.103 into equation A.106 results in [ m / V ( . - 0 J ( + [ q / > ( - i ) « ] ( r , ) + [K] J*WCdC + [C] <*c] {T,+i} - j T W {f} <*C = 0 (A.107) A l l the terms in equation A.107 are then divided by JQ W d£, and by introducing a new variable 6 where Jo W(d( e = Jo Wd( 229 equation A . 107 becomes '[K] (l-0) + [C] ( - ^ ) ] {T,} + [K] 0 + [C] (—)] {Tl+1} - {/} = 0 (A.108) where T r r / p W{f} d( Jo1 Wd( The same interpolation used for temperature can be applied to the vector {/} to give {f} = {fi+i} 0 + {fi} (1 - 8) Since the convective/radiative boundary conditions do not change significantly during time stepping, {//} and {/z+i} are approximately equal and thus {/} f» {/}. A value of 6 = | , known as the Crank-Nicholson method, generally gives the best solution [39]. Inserting this value of 0 in equation A.108, the final form of the matrix equation to be solved is M { 2 - % ^ } + . [ C ] { ^ ^ i } « { / } (A.109) A.1 .3 .2 Three-Poin t Recurrence Scheme In the three-point scheme, three interpolation functions are used and hence the shape functions are quadratic as shown in Figure A.96. The dimensionless time variable and the shape functions are - i < C < i c = & i V / _ 1 = - | c ( i - C ) ; NI.1 = (-\ + c) & Nt = ( l - 0 ( i + 0; Nt = - | £ (A.110) A i + 1 = | C ( i + 0; ^ = (1 + 0 h 230 Figure A.96: Quadratic shape functions (three-point recurrence formulae) for the time domain [57] 231 Temperature is thus approximated by {T} = N,-r + {T,} N, + {T,+1} Nl+1 (A.111) where the subscripts / — 1, / and / + 1 denote previous, current and next time step, respectively. Substituting Equation A . I l l into Equation A.101 and applying the weighted residual method gives £ W { [K] ( { 2 U } J V , - ! + {Ti} N, + {Tl+1} Nl+1) + [C] ({r,_!} TY/.! + {T,} N, + {7} + 1} Nl+1) - {/}} dC = 0 (A.112) Going through the same derivation as done in the two-point recurrence scheme, Equa-tion A.112 can be written as [K] At2 (I + B - 7) + [C] At (1 - 7)] {T,.!} + [K] At2 (I - 2/3 + 7 ) + [C] At (1 - 2 7 ) ] {2}} + [[K] At2 8 + [C] A i 7 ] {Tl+1} - {/} A * 2 = 0 (A.113) where 7 0 _ Hi w(c + $) d( Hi Wd( ji, w ( | c ( i + Q) d( I—i wac Hi w {/} d( {/} Hi Wd( The same interpolation used for temperature is applied to the vector {/} to give 1 7 } = { / < - i } ( } + / ? - 7 ) + { / / } Q - 2 / ? + 7 ) + U + i } / ? 232 Since the convective/radiative boundary conditions do not change significantly during time stepping, {//} and are approximately equal and thus {/} {/}. The values of 7 and /3 used are for the Dupont form with 7 = 1 and /? = | . The method is known as the Dupont three-level technique and has been successfully used in the past [12]. The final form of equation A.113 using the Dupont values is [ K ] | 3 n » + n - , i + [ c ] „ { / } (A.iH) A.2 Cartesian Coordinates (r, z) The same procedure carried out in the cylindrical system is applied to equation d2T d2T dT The derived equations in the cartesian coordinates are only slightly different from those in the cylindrical coordinates. The r terms in the elemental matrices and vectors derived in the cylindrical system are removed from the integrals and the (r, z) coordinates are replaced by the (x, y) coordinates. The resulting elemental matrices and vectors in the cartesian system are K $ = S y L k ( ^ d ^ + ^^)dxdy Kj$ = J r he Ni Nj dT KM = Jr hr Ni Nj dT C\f=jy!xPCpNiNjdxdy fhi — Ir hc Ni T O C l dT fi^frKNiT^dT or in numerical integration form m m j, K}] = E E Wi W3 k [B^] [B^] det[JUi^} (A.117) i=l j = l (A.116) 233 [4 e ) ] = E W i h {NM} {NM}T AM (A.118) 1=1 [Kie)] = jtwihr {NM} {NM}T A(„.) (A.119) »=I [ c « ] = E E ^ W , - , C P { % ^ i ) } T d e t [ J « " - « ) ] ( A - 1 2 0 ) m {/ie)} = E fcc { * < „ > } ^ A(„.) (A.121) i = l m {/r(e)} = E W i K {NM} Too2 A ( , ) (A.122) 1=1 j = i Appendix B Analytical Solutions The analytical solutions derived below are used to check the finite element calculations. They are listed according to the order presented in chapter 3. B. l Steady-State Conduction B . l . l One-Dimensional Solution The solution is for a hollow tube with fixed temperatures at its inner and outer surface. The equation to be solved is Integrating equation B.123 gives dT r — = const! (B.124) dr where const} is a constant. Integrating equation B.124 results in T = cons^ ln(r) + const2 (B.125) where const2 is another constant. consti and const2 are determined by applying boundary conditions. The boundary conditions are T(r t ) = T( 234 235 T(ra) = T0 where r, and ra are the inner and outer radius, and Ti and TQ are the inner and outer temperature, respectively. After substituting the boundary conditions into equation B.125, the solution is [27] T(r) =Ti + (T0 - Ti) — ^ (B.126) B.1 .2 Two-Dimens iona l Solutions The solution is for a cylinder 0 < r < b and 0 < z < c, when the boundary at r = b is at temperature f(z) and the boundaries at z = 0 and z = c are at zero temperature. The mathematical formulation is given as + 0 = 0 m 0 < r < 6 , 0<z<c T = Hz) at r = b (B.127) T = 0 at z = 0 and z = c Using the separation of variables and applying the boundary conditions, the resulting solution is [54] T(r,z) = - f ) ^ 4 MVmz) [C sm{r,mz!)f(z')dz' (B.128) where m7r ?/m = C and IQQ is the modified Bessel function of order 0 of the first kind. 236 B.1.2.1 B o u n d a r y Temperature Profi le 1 Given that / M = T - z (B.129) f(z) = Tmid c - Tmid z \c< z <c equation B.128 becomes T ( r , 2 ) = _ £ _ . | _ s.n ( _ j - _ cos ( — j j C m = l 0 < Z < -c ~ ~ 2 £ ^ e ) + - £ _ „ , . ( ! = £ ) } ( B . 1 3 0 ) — c < z < c 2 B . l . 2 . 2 B o u n d a r y Temperature Profi le 2 Given that f(z) - Tmid sin (~^) equation B.128 becomes B . 2 Steady-State Convect ion — Two-Dimens iona l Solu t ion The solution is for a cylinder 0 < r < b and 0 < z < c, when the boundary at z — 0 is at temperature f(r), boundary at z = c is kept at zero temperature, and that at r = b 237 dissipates heat by convection into a medium at zero temperature. The mathematical formulation is given as 0 + r - f + 0 = O i n 0 < r < 6 , 0 < z < c % + = 0 at r = i * (B.132) T = f(r) at 2 = 0 where T — 0 at z — c Given that f(r) = T^o = constant, using the separation of variables and applying the boundary conditions, the resulting solution is [54] H r ' Z ) ~ b + sinh(/?mc) J 0(/3m6) ( B ' 1 3 3 ) where /?m's are the positive roots of A„ Ji(0ra&) = HJ0(l3mb) and Jo() and Ji() are the Bessel functions of order 0 and 1 of the first kind, respectively. B.3 Transient Conduction — One-Dimensional Solutions B.3.1 Case 1 The equation to be solved and the boundary conditions are * ^ = I « £ a i n O < * < L , t > 0 fX = 0 at x = 0 and x = L, t> 0 (B.134) T = F(x) for * = 0, in 0 < x < L 238 The corresponding solution using the separation of variables after applying the boundary conditions is [54] T(x,t) = ±r / F(x')dx'++±r £ <Talimt.cos{Pmx) / F(x') cos(Pmx')dx' (B.135) L Jo L Jx>=o where TH 7T Pm = m = 1,2,3, Given that F(x) is linear, i.e. F(x) = T0x the solution becomes cos(/?mL) 1 T(x,t) = ^ + ^ £ e - ^ . c o s ( / ? m * ) j i - L sin(/?mL) + Pm Pr, (B.136) B.3 .2 Case 2 The equation to be solved and the boundary conditions are £ g a = I « £ f l . in 0 < * < o o , r > 0 T = 0 at i = 0 , i > 0 (B.137) T = F(x) for t = 0, in 0 < x < oo Given that .F(a;) = T 0 = constant, the corresponding solution using the separation of variables after applying the boundary conditions is [54] T(x,t)=T0erf^=j (B.138) where er/() is the error function. 239 B.4 Transient Convection — One-Dimensional Solutions B.4.1 Case 1: Cartesian System The equation to be solved and the boundary conditions are ^ - = Z ^ « 0 < x < L, t > 0 f£ = 0 at x = 0, i > 0 g + HT = 0 at i = I , i > 0 (B.139) T = F(x) for t = 0, in 0 < x < L Using the separation of variables and applying the boundary conditions, the solution is [54] °° •> Q2 4- H2 rL T(x,t) = 2 £ e~al3mt L{/32+H2) + H co<^x) £=/(*') cos((3mx')dx' (B.140) where /?m's are the positive roots of /9m tan(/3 r aL) = H For the case where F(x) — To = constant, the solution becomes B.4.2 Case 2: Cylindrical System For a body with negligible internal resistance where it is assumed that no thermal gradient exists within the body, the transient response is governed by [27] -VpCp^ = hA(T-T00) where A and V are the surface area and volume of the body respectively, and T is a function of time only. When Too = 0, the equation can be directly integrated and the 240 resulting solution is where To is the initial temperature of the body. Appendix C Calculation of the Rate of Temperature Increase in the Liquid Gallium inside the Middle Compartment Assuming No Heat Loss Total rate of heat input into middle compartment, dqtotai/dt = P.R, where I = electrical current w 5 amperes, and R = resistance of Pt wire « 0.009064 ohm. Thus, dqtotai/dt « 5 2 x 0.009064 « 0.2266 watt. Assuming heat stays in liquid gallium inside middle compartment, Rate of heat input into middle compartment, dqtotai/dt = m.Cp.dT/dt, where m = mass of liquid gallium in middle compartment w 0.00287 kg, Cp = specific heat capacity of liquid gallium w 372 J/kg.°C, and dT/dt — rate of temperature rise in middle compartment. Therefore, dT/dt w 0.2266/(0.00287 x 372) » 12.75 °C/minute . 241 Appendix D Calculation of the Rate of Temperature Increase in the Liquid Gallium inside the Side Compartment Assuming Uni-Directional Heat Flow across the Specimen Rate of heat input into middle compartment, dqi/dt = rrii .Cp.dTi / dt, where mi = mass of liquid gallium in middle compartment « 0.00287 kg, Cp — specific heat capacity of liquid gallium w 372 J/kg.°C, and dT\jdt = rate of temperature rise in middle compartment w 1.5 °C/minute . Thus, dqx/dt « 0.00287 x 372 x 1.5/60 » 0.0267 watt. Total rate of heat input into middle compartment, dqtotai/dt = P.R, where I = electrical current « 5 amperes, and R = resistance of Pt wire ~ 0.009064 ohm. Thus, dqtotaildt « 5 2 x 0.009064 « 0.2266 watt. Assuming heat only flows through specimen plates, Rate of heat input into side compartment, dq-ijdt — 1/2 x (dqtotai/dt - dqi/dt). Thus, dq2/dt w 1/2 x (0.2266 - 0.0267) » 0.1'watt. But dq2/dt = m2.Cp.dT2/dt, where m2 — mass of liquid gallium in side compartment « 0.035 kg, and dT2/dt = rate of temperature rise in side compartment. Therefore, dT2/dt » 0.1/ (0.035 x 372) » 0.461 0 C/minute. 242 Bibliography [1] Standard Method of Test for Thermal Conductivity of Materials by Means of the Guarded Hot Plate (C177-71), Annual Book of ASTM Standards. Volume 18, A S T M , American Society for Testing and Materials, Philadelphia, 1974. [2] Summary of Pyrolytic Boron Nitride Properties. Union Carbide Corporation. [3] A . S. Jordan. A n Evaluation of the Thermal and Elastic Constants Affecting GaAs Crystal Growth. Journal of Crystal Growth, 49:631-642, 1980. [4] A . S. Jordan. Estimated Thermal Diffusivity, Prandtl Number and Grashof Number of Molten GaAs, InP, and GaSb. Journal of Crystal Growth, 71:551-558, 1985. [5] A . S. Jordan. Some Thermal and Mechanical Properties of InP Essential to Crystal Growth Modeling. Journal of Crystal Growth, 71:559-565, 1985. [6] A . S. Jordan, A . R. Von Neida and R. Caruso. The Theoretical and Experimental Fundamentals of Decreasing Dislocations in Melt Grown GaAs and InP. Journal of Crystal Growth, 76:243-262, 1986. [7] A . S. Jordan, A . R. Von Neida and R. Caruso. The Theory and Practice of Dislo-cation Reduction in GaAs and InP. Journal of Crystal Growth, 70:555-573, 1984. [8] A . S. Jordan, R. Caruso, A . R. Von Neida, and J . W. Nielsen. A Comparative Study of Thermal Stress Induced Dislocation Generation in Pulled GaAs, InP, and Si Crystals. Journal of Applied Physics, 52(5):3331-3336, May 1981. [9] A . S. Jordan, R. Caruso, and A . R. Von Neida. A Thermoelastic Analysis of Dis-location Generation in Pulled GaAs Crystals. The Bell System Technical Journal, 59(4):593-637, Apri l 1980. [10] A . W. Moore. Characterization of Pyrolytic Boron Nitride for Semiconductor Ma-terials Processing. Journal of Crystal Growth, 106:6-15, 1990. [11] Alan J . Chapman. Heat Transfer. Macmillan PubHshing Company, 866 Third Avenue, New York, New York 10022, Fourth edition, 1984. [12] B . G . Thomas, I. V . Samarasekera, and J . K . Brimacombe. Comparison of Numerical Modeling Techniques for Complex, Two-Dimensional, Transient Heat-Conduction Problems. Metallurgical Transactions B, 15B:307-318, June 1984. 243 244 [13] Brian R. Pamplin, editor. Crystal Growth. Volume 16 of International Series in the Science of the Solid State, Pergamon Press, Oxford, Second edition, 1980. [14] C. E . Harrison. M. A. Sc. Thesis. The University of British Columbia, Apri l 1984. [15] C. L . Parfeniuk. M. A. Sc. Thesis. The University of British Columbia, November 1990. [16] C. Schvezov, I. V . Samarasekera and F . Weinberg. Calculation of the Shear Stress Distribution in L E C Gallium Arsenide for Different Growth Conditions. Journal of Crystal Growth, 92:479-488, 1988. [17] C. Schvezov, I. V . Samarasekera and F . Weinberg. Calculation of the Thermal and Stress Fields in L E C Gallium Arsenide after Separation from the Melt. Journal of Crystal Growth, 92:489-497, 1988. [18] C. Schvezov, I. V . Samarasekera and F . Weinberg. Mathematical Modelling of the Liquid Encapsulated Czochralski Growth of Gallium Arsenide I: Heat Flow Model. Journal of Crystal Growth, 84:212-218, 1987. [19] C. Schvezov, I. V . Samarasekera and F. Weinberg. Mathematical Modelling of the Liquid Encapsulated Czochralski Growth of Gallium Arsenide II: Stress Model. Journal of Crystal Growth, 84:219-230, 1987. [20] C. Schvezov, I. V . Samarasekera and F . Weinberg. The Effect of Crystal Radius Fluctuations on the Stress Field in L E C Gallium Arsenide. Journal of Crystal Growth, 85:142-147, 1987. [21] D. E . Bornside, T. A . Kinney, R. A . Brown and George K i m . Finite Element/ Newton Method for the Analysis of Czochralski Crystal Growth with Diffuse-Grey Radiative Heat Transfer. International Journal for Numerical Methods in Engineer-ing, 30:133-154, 1990. [22] D. T. J . Hurle. Analytical Representation of the Meniscus in Czochralski Growth. Journal of Crystal Growth, 63:13-17, 1983. [23] D. W. Pepper and A . J . Baker. Finite Differences versus Finite Elements. Handbook of Numerical Heat Transfer, John Wiley and Sons, Inc., 1988. [24] E . Bill ig. Growth of Monocrystals of Germanium from an Undercooled Melt. Proc. R. Soc. Lond. A, 229:346-363, 1955. [25] E . K . Sichel and R. E . Miller. Thermal Conductivity of Highly Oriented Pyrolytic Boron Nitride. Proc. l^th International Conference on Thermal Conductivity, 11-17, 1975. 245 [26] F . Dupret, Y . Ryckmans, P. Wouters and M . J . Crochet. Numerical Calculation of the Global Heat Transfer in a Czochralski Furnace. Journal of Crystal Growth, 79:84-91, 1986. [27] F . Kreith and W. Z. Black. Basic Heat Transfer. Harper and Row, New York, 1980. [28] G. Strang and G. J . Fix. An Analysis of the Finite Element Method. Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1973. [29] G . Williams and W. E . Reusser. Heat Transfer in Silicon Czochralski Crystal Growth. Journal of Crystal Growth, 64:448-460, 1983. [30] H . Kopetsch. A Numerical Method for the Time-dependent Stefan Problem in Czochralski Crystal Growth. Journal of Crystal Growth, 88:71-86, 1988. [31] H . S. Carslaw and J . C. Jeager. Conduction of Heat in Solids. Oxford University Press, Glasgow, 1959. [32] Haruo Emori, Takao Matsumura, Tosliio Kikuta and Tsuguo Fukuda. Effect of Ambient Gas on Undoped L E C GaAs Crystal. Japanese Journal of Applied Physics, 22(11):1652-1655, November 1983. [33] Hiroki Kohda, Kohji Yamada, Hideo Nakanishi, Takashi Kobayashi, Jiro Osaka and Keigo Hoshikawa. Crystal Growth of Completely Dislocation-Free and Striation-Free GaAs. Journal of Crystal Growth, 71:813-816, 1985. [34] Huebner and Thornton. The Finite Element Method for Engineers. John Wiley and Sons, Inc., New York, Second edition, 1982. [35] I. Hatta. Thermal Diffusivity Measurements of Thin Films and Multilayered Com-posites. International Journal of Thermophysics, ll(2):293-303, 1990. [36] J . J . Derby and R. A . Brown. Thermal-Capillary Analysis of Czochralski and Liquid Encapsulated Czochralski Crystal Growth I: Simulation. Journal of Crystal Growth, 74:605-624, 1986. [37] J . J . Derby, L . J . Atherton, P. D. Thomas, and R. A . Brown. Finite-Element Methods for Analysis of the Dynamics and Control of Czochralski Crystal Growth. Journal of Scientific Computing, 4(2):297-343, 1987. [38] J . J . Derby, R. A . Brown, F . T. Geyling, A . S. Jordan, and G. A . Nikolakopoulou. Finite Element Analysis of a Thermal-Capillary Model for Liquid Encapsulated Czochralski Growth. J. Electrochem. Soc, 32(2):470-482, 1985. 246 [39] J . M . Bettencourt, 0 . C. Zienkiewicz and G. Cantin. Consistent Use of Finite Elements in Time and the Performance of Various Recurrence Schemes for the Heat Diffusion Equation. International Journal for Numerical Methods in Engineering, 17:931-938, 1981. [40] J . Osaka, H . Kohda, T. Kobayashi and K . Hoshikawa. Japanese Journal of Applied Physics, 23:L195, 1984. [41] John C. Lambropoulos and Christopher N . Delametter. The Effect of Interface Shape on Thermal Stress during Czochralski Crystal Growth. Journal of Crystal Growth, 92:390-396, 1988. [42] K . Hoshikawa, H . Kohda and H . Hirata. Japanese Journal of Applied Physics, 23:L37, 1984. [43] K . Kakimoto and T. Hibiya. Temperature Dependence of Viscosity of Molten GaAs by an Oscillating Cup Method. Appl. Phys. Lett., 50(18):1249-50, 1987. [44] K . Matsumoto, H . Morishita, M . Sasaki, S. Nishine, M . Yokogawa, M . Sekinobu, K . Tada and S. Akai. L E C Growth of 3-Inch Diameter Undoped SI GaAs with Good Reproducibility Using an Improved Hot Zone. Semi-Insulating III- V Materials Third Conference, Kah-Nee-Ta, Oregon, 175-177, 1984. [45] K . Terashima and T. Fukuda. A New Magnetic-Field Applied Pulling Apparatus for L E C GaAs Single Crystal Growth. Journal of Crystal Growth, 63:423-425, 1983. [46] Kazutaka Terashima, Johji Nishio, Shoichi Washizuka and Masayuki Watanabe. Magnetic Field Effect on Residual Impurity Concentrations for L E C GaAs Crystal Growth. Journal of Crystal Growth, 84:247-252, 1987. [47] Kenji Hiraga, Takeo Oku, Makoto Hirabayashi, Toshitsugu Matsuda and Toshio Hi -rai. Fivefold Multiply-Twinned Precipitates in Chemically Vapour-Deposited Boron Nitride Studied by Transmission Electron Microscopy. Journal of Materials Science Letters, 8:130-134, 1989. [48] L . J . Atherton, J . J . Derby and R. A . Brown. Radiative Heat Exchange in Czochral-ski Crystal Growth. Journal of Crystal Growth, 84:57-78, 1987. [49] L . N . Hjellming and J . S. Walker. Modelling of Melt Motion in a Czochralski Crystal Puller with an Axial Magnetic Field. Journal of Crystal Growth, 79:96-98, 1986. [50] M . J . Crochet, P. J . Wouters, F . T. Geyling and A . S. Jordan. Finite-Element Simulation of Czochralski Bulk Flow. Journal of Crystal Growth, 65:153-165, 1983. 247 [51] M . Mihelcic and K . Wingerath. Numerical Simulations of the Czochralski Bulk Flow in an Axial Magnetic Field: Effects on the Flow and Temperature Oscillations in the Melt. Journal of Crystal Growth, 71:163-168, 1985. [52] M . Mihelcic and K . Wingerath. Three-Dimensional Simulations of the Czochralski Bulk Flow in a Stationary Transverse Field and in a Vertical Magnetic Field: Effects on the Asymmetry of the Flow and Temperature Distribution in the Si Melt. Journal of Crystal Growth, 82:318-326, 1987. [53] M . Mihelcic, K . Wingerath and Chr. Pirron. Three-Dimensional Simulations of the Czochralski Bulk Flow. Journal of Crystal Growth, 69:473-488, 1984. [54] M . Necati Ozisik. Heat Conduction. John Wiley and Sons, Inc., New York, 1980. [55] N . Naito and C. H . Hsueh. Residual Stress and Strain in Pyrolytic Boron Nitride Resulting from Thermal Anisotropy. Journal of Materials Science, 23:1901-1905, 1988. [56] Nobuyuki Kobayashi. Heat Transfer in Czochralski Crystal Growth. Volume 6 of Preparation and Properties of Solid State Materials, Marcel Dekker, Inc., 1981. [57] O. C. Zienkiewicz. The Finite Element Method. McGraw-Hill Book Company, New York, 1977. [58] O. C. Zienkiewicz and Morgan. Finite Elements and Approximation. John Wiley and Sons, Inc., New York, 1983. [59] P. A . Ramachandran and M . P. Dudukovic. Simulation of Temperature Distribution in Crystals Grown by Czochralski Method. Journal of Crystal Growth, 71:399-408, 1985. [60] P. A . Sackinger, R. A . Brown, and J . J. Derby. A Finite Element Method for Analysis of Fluid Flow, Heat Transfer and Free Interfaces in Czochralski Crystal Growth. International Journal for Numerical Methods in Fluids, 9:453-492, 1989. [61] P. D. Thomas, J . J . Derby, L . J . Atherton and R. A . Brown. Dynamics of Liquid-Encapsulated Czochralsski Growth of Gallium Arsenide: Comparing Model and Experiment. Journal of Crystal Growth, 96:135-152, 1989. [62] P. Nicodeme, F . Dupret and M . J . Crochet. Effect of Geometrical Parameters upon the L E C Growth of GaAs Crystals. Semi-Insulating III- V Materials Fifth Confer-ence, Malmo, Sweden, 465-470, 1988. 248 63] P. Nicodeme, F . Dupret, M . J . Crochet, J . P. Farges, and G. Nagel. Numerical Simulation of Heat Transfer in L E C Growth of Gallium Arsenide. Journal of Crystal Growth, 97:162-172, 1989. 64] P. Sabhapathy and M . E. Salcudean. Numerical Analysis of Heat Transfer in L E C Growth of GaAs. Journal of Crystal Growth, 97:125-135, 1989. 65] P. Sabhapathy and M . E . Salcudean. Numerical Study of Flow and Heat Transfer in L E C Growth of GaAs with an Axial Magnetic Field. Journal of Crystal Growth, 104:371-388, 1990. [66] R. A . Brown, P. A . Sackinger, P. D. Thomas, J . J . Derby and L . J . Atherton. Inter-Disciplinary Issues in Materials Processing and Manufacturing. ASME, 331, 1987. [67] R. Balasubramaniam and S. Ostrach. Transport Phenomena near the Interface of a Czochralski-Grown Crystal. Journal of Crystal Growth, 88:263-281, 1988. [68] R. D. Cook, D. S. Malkus and M . E . Plesha. Concepts and Applications of Finite Element Analysis. John Wiley and Sons, Inc., New York, Third edition, 1989. [69] R. K . Srivastava, P. A . Ramachandran and M . P. Dudukovic. Interface Shape in Czochralski Grown Crystals: Effect of Conduction and Radiation. Journal of Crystal Growth, 73:487-504, 1985. [70] R. K . Srivastava, P. A . Ramachandran and M . P. Dudukovic. Simulation of Jet Cooling Effects on Czochralski Crystal Growth. Journal of Crystal Growth, 76:395-407, 1986. [71] R. Lamprecht, D. Schwabe, A . Scharmann and E . Schultheiss. Experiments on Buoyant, Thermocapillary, and Forced Convection in Czochralski Configuration. Journal of Crystal Growth, 65:143-152, 1983. [72] Robert A . Brown. Perspectives on Modelling Bulk Crystal Growth of Semiconduc-tors from the Melt. A paper distributed at the F. Weinberg International Symposium on Solidification Processing, June 1990. [73] Robert A . Brown. Theory of Transport Processes in Single Crystal Growth from the Melt. AIChE Journal, 34(6):881-911, June 1988. [74] Sybil P. Parker, editor. McGraw-Hill Dictionary of Scientific and Technical Terms. McGraw-Hill , Inc., New York, Fourth edition, 1989. 75] T. R. AuCoin, R. L. Ross, M . J . Wade and R. O. Savage. Solid State Technology, 59:, 1979. 249 [76] Tadashi Kimura, Tooru Katsumata, Masato Nakajima and Tsuguo Fukuda. The Ef-fect of Strong Magnetic Field on Homogeneity in L E C GaAs Single Crystal. Journal of Crystal Growth, 79:264-270, 1986. [77] Takao Tsukada, Nobuyuki Imaishi and Mitsunori Hozawa. Theoretical Study of the Flow and Temperature Fields in CZ Single Crystal Growth. Journal of Chemical Engineering of Japan, 21(2):184-191, July 1987. [78] Takashi Fujii, Minoru Eguchi, Tomoki Inada and Tsuguo Fukuda. Effects of Crucible Material on Undoped L E C GaAs Crystals. Journal of Crystal Growth, 84:237-240, 1987. [79] Toshitsugu Matsuda, Naoki Uno, Hiroyuki Nakae and Toshio Hirai. Synthesis and Structure of Chemically Vapour-Deposited Boron Nitride. Journal of Materials Sci-ence, 21:649-658, 1986. [80] W. E. Langlois. A Parameter Sensitivity Study for Czochralski Bulk Flow of Silicon. Journal of Crystal Growth, 56:15-19, 1982. [81] W. E . Langlois. Computer Simulation of Czochralski Melt Convection in a Magnetic Field. Journal of Crystal Growth, 70:73-77, 1984. [82] W. E. Langlois and J. S. Walker. Czochralski Crystal Growth in an Axial Magnetic Field. Proc. Second BAIL Conference, Boole, Dublin, 299-304, 1982. [83] W. E. Langlois and Ki-Jun Lee. Czochralski Crystal Growth in an Axial Magnetic Field: Effects of Joule Heating. Journal of Crystal Growth, 62:481-486, 1983. [84] W. E . Langlois, L . N . Hjellming and J . S. Walker. Effects of the Finite Electrical Conductivity of the Crystal on Hydromagnetic Czochralski Flow. Journal of Crystal Growth, 83:51-61, 1987. [85] W. R. Wilcox and R. L . Duty. Macroscopic Interface Shape during Solidification. Journal of Heat Transfer, 88c:45-51, 1966. [86] W. Zulehner. Czochralski Growth of Silicon. Journal of Crystal Growth, 65:189-213, 1983. [87] Willie Fajber and James Fajber. Preliminary Temperature Measurements of Liquid Germanium under a Magnetic Field. Department of Metals and Materials Engineer-ing, The University of British Columbia, 1990. [88] Y . S. Touloukian and C. Y . Ho, editor. Thermophysical Properties of Matter: Theory of Thermal Conductivity of Nonmetallic Solids. Volume 2 of The TPRC Data Series, IFI/Plenum Data Corporation, New York, 1970.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Heat flow in the melt and crucible for crystal growth
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Heat flow in the melt and crucible for crystal growth Kho, Rowin Wisadi 1991
pdf
Page Metadata
Item Metadata
Title | Heat flow in the melt and crucible for crystal growth |
Creator |
Kho, Rowin Wisadi |
Publisher | University of British Columbia |
Date Issued | 1991 |
Description | A magnetic field is imposed on the melt during Czochralski crystal growth to control fluid flow. The influence of an applied magnetic field on heat transfer in a germanium melt has been investigated in this work. In heat-transfer modelling of the Czochralski process, values of thermal conductivity of the crucible material, pyrolytic boron nitride, are required as a function of temperature. The thermal conductivity of a pyrolytic boron nitride crucible material has also been determined in this work. Both the investigation of the influence of the applied magnetic field and the determination of the thermal conductivity have involved experimental temperature measurements and mathematical modelling of heat transfer. The finite element equations for the two-dimensional heat-conduction equation have been derived using the Galerkin method of weighted residuals. They have been validated by comparing the finite element solutions with the corresponding analytical solutions. The finite element results are in excellent agreement with the analytical solutions. A steady-state conduction-dominated mathematical model has been developed to analyze the temperature measurements obtained within the germanium melt in a Czochralski crystal growth configuration with and without an applied magnetic field of 0.099 tesla. The effect of the applied magnetic field on the heat transfer in the melt has been determined by fitting the model-calculated results to the measurements. The effective thermal conductivity of the melt has been found to decrease by a factor of seven due to the magnetic field. The thermal conductivity across a pyrolytic boron nitride crucible plate (through crucible thickness) has been determined from measurements of temperature responses in liquid gallium positioned on both sides of the plate, in conjunction with a transient mathematical model which simulates the thermal responses. By matching the simulated thermal responses with the measurements, the thermal conductivity of the pyrolytic boron nitride crucible plate has been obtained as a function of temperature. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-11-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0078450 |
URI | http://hdl.handle.net/2429/29991 |
Degree |
Master of Applied Science - MASc |
Program |
Materials Engineering |
Affiliation |
Applied Science, Faculty of Materials Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1991_A7 K56.pdf [ 11.79MB ]
- Metadata
- JSON: 831-1.0078450.json
- JSON-LD: 831-1.0078450-ld.json
- RDF/XML (Pretty): 831-1.0078450-rdf.xml
- RDF/JSON: 831-1.0078450-rdf.json
- Turtle: 831-1.0078450-turtle.txt
- N-Triples: 831-1.0078450-rdf-ntriples.txt
- Original Record: 831-1.0078450-source.json
- Full Text
- 831-1.0078450-fulltext.txt
- Citation
- 831-1.0078450.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0078450/manifest