Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Mathematical modelling of the vertical Bridgman growth of cadmium telluride Parfeniuk, Christopher Luke 1990

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1991_A7 P37.pdf [ 15.49MB ]
Metadata
JSON: 831-1.0078458.json
JSON-LD: 831-1.0078458-ld.json
RDF/XML (Pretty): 831-1.0078458-rdf.xml
RDF/JSON: 831-1.0078458-rdf.json
Turtle: 831-1.0078458-turtle.txt
N-Triples: 831-1.0078458-rdf-ntriples.txt
Original Record: 831-1.0078458-source.json
Full Text
831-1.0078458-fulltext.txt
Citation
831-1.0078458.ris

Full Text

MATHEMATICAL MODELLING OF T H E VERTICAL BRIDGMAN GROWTH OF CADMIUM TELLURIDE By Christopher Luke Parfeniuk B. A. Sc. , The University of British Columbia, 1988 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF APPLIED SCIENCE in T H E FACULTY OF GRADUATE STUDIES METALS AND MATERIALS ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA November 1990 © Christopher Luke Parfeniuk , 1990 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 scholarly 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. Metals and Materials Engineering The Universit}' of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: Abstract The temperature and stress fields for the Vertical Bridgman growth of CdTe were calcu-lated using the finite element method. The critical resolved shear stress (CRSS) of CdTe was measured at elevated temperatures. The measured CRSS and the calculated shear stress were compared to determine the onset of dislocation formation in the crystal. The effect of furnace gradient, crucible velocity, hot furnace temperature and cold furnace temperature on the solid/liquid interface shape and calculated resolved shear stress were examined using the model. Etch pit density measurments were carried out in order to verify the stress model. 11 Table of Contents Abstract ii List of Figures viii List of Tables xvii List of Symbols xix Acknowledgement xxiii Dedication xxiv 1 Introduction 1 1.1 Crystal Growth in Crucibles 1 1.2 Vertical Bridgman Growth 2 2 Literature Review . 9 2.1 Thermal Models 9 2.2 Stress Models and Dislocation Studies 10 2.3 CRSS of CdTe 12 3 Objectives 15 4 Critical Resolved Shear Stress Experimental Procedure 17 4.1 Sample Preparation 17 4.1.1 Polishing 17 iii 4.1.2 Crystal Orientation Measurements 18 4.1.3 Shaping the Tensile Samples 18 4.2 Tensile Testing Apparatus 19 4.3 Tensile Testing Procedure 20 5 Etch Pit Experimental Procedure 22 5.1 Etching 22 5.2 Image Analysis 24 6 Critical Resolved Shear Stress Experimental Results 26 6.1 Crystal Orientations 26 6.2 Tensile Tests 26 6.3 Critical Resolved Shear Stress ' 31 7 Etch Pit Density and Distribution 33 8 Interface Position and Shape 56 9 Mathematical Models 65 9.1 Scope of Models and Assumptions 65 9.2 Idealized Domain 66 9.3 Temperature Model 68 9.3.1 Boundary Conditions 72 9.3.2 Phase Transformation Solution 74 9.3.3 Solidification Gap 75 9.4 Stress Model 77 9.4.1 Critical Resolved Shear Stress 79 9.4.2 Boundary Conditions 79 iv 9.5 Thermo-Physical Properties 80 10 Modelling Results 83 10.1 Comparisons to Analytical Solutions 83 10.1.1 Conduction 83 10.1.2 Convection 85 10.1.3 Transient 90 10.1.4 Phase Change 94 10.1.5 Stress/Strain 99 10.2 Sensitivity Analysis 100 10.2.1 Mesh Size 100 10.2.2 Time Step Size 103 10.2.3 Latent Heat 103 10.2.4 Conductivity 109 10.3 Variables Examined and Operating Ranges 116 10.4 Model Predictions Assuming No Fluid Flow 118 10.4.1 Effects of Crucible Velocity and Furnace Gradient 119 10.4.2 Effects of Cold and Hot Zone Temperature 122 10.4.3 Excess Shear Stress Results 134 10.5 Model Predictions using a Bulk Heat Transfer Coefficient Seven Times the Static Value 138 10.5.1 Effects of Crucible Velocity and Furnace Gradient 140 10.5.2 Effects of Cold and Hot Zone Temperature 142 10.5.3 Excess Shear Stress Results 151 11 Discussion 156 11.1 Discussion of Experimental Results 156 v 11.2 Discussion of Model Predictions Assuming no Fluid Flow 157 11.3 Discussion of Model Predictions using a Bulk Heat Transfer Coefficient Seven Times the Value with no Fluid Flow 158 11.4 Comparison Between the Two Sets of Model Results 159 12 Summary and Conclusions 160 13 Recommendations for Future Work 163 Bibliography 164 A Formulation of Thermal F E M Model 167 A. l Galerkins Method 167 A.2 Finite Element Method 172 A.3 Shape (Trial) Function 172 A.4 Numerical Integration and Mapping Functions 175 A.5 Conduction 178 A.6 Convection 180 A.7 Time Stepping Techniques 181 A.7.1 Capacitance Matrix 182 A.7.2 Time Stepping Using Recurrence Technique 182 A.7.3 Two Point Recurrence Scheme 182 A.7.4 Three Point Recurrence Scheme 185 A.8 Phase Transformation Solution 188 A.9 Solidification Gap 192 B Formulation of F E M for Stress Fields 195 vi C Grashof Number Calculation 203 C l Grashof Number for Tin 203 C.2 Grashof Number for CdTe 204 C.2.1 Theoretical CdTe viscosity calculation 205 C.2.2 Grashof number 208 D Relative Thermal Resitance 209 vn List of Figures 1.1 The growth of spurious crystals nucleated at the intersection of the solid-liquid interface with the crucible wall, (a) convex interface, the spurious crystals are grown toward the crucible surface, (b) concave interface, spu-rious crystals grow inward and dominate crystal growth 3 1.2 Vertical Bridgman Growth 4 1.3 Phase Diagram of CdTe [1] 5 1.4 Representation of Solid CdTe Phase Region [1] 6 1.5 Different Ampule Shapes for Vertical Bridgman Crystal Growth (a) Con-ical Bottom (b) Capillary tip attached on bottom 7 2.6 Critical resolved shear stress of CdTe [15] 13 2.7 Stress/strain curves of CdTe single crystals of identical orientations de-formed at the temperatures indicated. Orientation of compression axis of sample is shown in stereographic triangle [16] 14 4.8 Tensile sample(A) and tensile test apparatus(B) 21 5.9 Crystal structure of CdTe (Inoue et al [19] ) 23 6.10 Wulf net for 111 plane in Cdo.96Z1io.04Te with force direction for each tensile test 27 6.11 Stress-Displacement curve for a tensile test of Cdo.96Zno.04Te at 932K. . . 29 6.12 Stress Displacement Curves for tensile tests of Cdo.96Zno.04Te 30 6.13 Dependence of critical resolved shear stress on temperature for Cdo.96Zno.01Te. 32 viii 7.14 Location of slices in boules 521 and 555 37 7.15 Etch pit density map for sample 521s01, Mag x 3 38 7.16 Areas of sample shown in Figures 7.17a, to 7.17c 39 7.17aEtch surface of sample 521s01 at area A, Mag x 80 40 7.17bEtch surface of sample 521s01 at area B, Mag x 80 40 7.17cEtch surface of sample 521s01 at area C, Mag x 80 . . . 41 7.18 Etch pit density map for sample 521sll, Mag x 3 42 7.19 Areas of sample shown in Figures 7.20a to 7.20f 43 7.20aEtch surface of sample 521sll at area A, Mag x 80 44 7.20bEtch surface of sample 521sll at area B, Mag x 80 44 7.20cEtch surface of sample 521sll at area C, Mag x 80 45 7.20dEtch surface of sample 521sll at area D, Mag x 80 45 7.20eEtch surface of sample 521sll at area E, Mag x 80 46 7.20f Etch surface of sample 521sll at area F, Mag x 80 46 7.21 Etch pit density map for sample 555sll, Mag x 3 47 7.22 Areas of sample shown in Figures 7.23a to 7.23f 48 7.23aEtch surface of sample 555sll at area A, Mag x 80 49 7.23bEtch surface of sample 555sll at area B, Mag x 80 49 7.23cEtch surface of sample 555sll at area C, Mag x 80 50 7.23dEtch surface of sample 555sll at area D, Mag x 80 . . . . ' 50 7.23eEtch surface of sample 555sll at area E, Mag x 80 51 7.23f Etch surface of sample 555sll at area F, Mag x 80 51 7.24 Etch pit density map for sample 555s28, Mag x 3 52 7.25 Areas of sample shown in Figures 7.26a to 7.26d 53 7.26aEtch surface of sample 555s28 at area A, Mag x 80 54 7.26bEtcli surface of sample 555s28 at area B, Mag x 80 54 ix 7.26cEtch surface of sample 555s28 at area C, Mag x 80 55 7.26dEtch surface of sample 555s28 at area D Mag x 80 55 8.27 Cdo.96Zno.04Te boule grown to examine interface profile 58 8.28 Top of boule were multicrystalline growth occured 59 8.29 Second Cdo.9eZno.04Te boule grown to examine interface profile 60 8.30 SIMS analysis for Zn, Cd and Te at Area 1 62 8.31 SIMS analysis for Zn, Cd and Te at Area 2 63 8.32 SIMS analysis for Zn, Cd and Te at Area 3 64 9.33 Model representation of domain 69 9.34 Eight noded quadratic isoparametric element 71 9.35 Furnace Temperature as seen by side of crucible (A) and as seen by the bottom and top of the crucible (B) 73 9.36 Boundary elements that are a soHdification gap 76 9.37 Strain and Stress involved in the analysis of axi-symmetric solids 78 10.38Domain used for testing the conduction matrix 84 10.39Temperature profile for a small temperature difference between inside and outside radius (T(r=o.oim) = 0.0JO T( r =i.6im) — l-OA') 86 10.40Temperature profile for a large temperature difference between inside and outside radius (T^o.oim) = 0.0IO T( r = 1.6im) = 60.0/0 8 7 10.41Temperature profile for a small temperature difference between inside and outside radius (T(r=o.oim) = 0.0/0 0^0 = 1-0K) 90 10.42Temperature profile for a large temperature difference between inside and outside radius (r ( r = 0.oim) = 0.0JO ^ = lOO.OiT) 91 10.43Transient response at r = 1.4 m for F.E.M. formulation 95 x 10.44Domain used for testing the latent heat solution 97 10.45Comparison to analytical solution for a phase change 98 10.46Comparison between finite element solution and analytical solution for thermal stress in a cylinder. Temperature variation: T = —r2K, Con-straint Fz = 0 99 10.47The three mesh densities used to examine the models sensitivity. (a)Fine mesh density, size approximately 1.25mm by 1.5mm. (b) Normal mesh density, size approximately 2.5mm by 3.0mm. (c) Coarse mesh density, size approximately 5.0mm by 6.0mm 101 10.48Axial interface location during crystal growth for different mesh densities (Interface as measured at center of crystal) 102 10.49Axial interface location during crystal growth for different time step sizes (Interface as measured at center of crystal) 104 10.50Axial interface location during crystal growth for different time step sizes (Interface as measured at side of crystal) 105 10.51Axial interface location for crystal growth including latent heat and ex-cluding latent heat (Interface as measured at center of costal) 106 10.52Axial interface location for crystal growth including latent heat and ex-cluding latent heat (Interface as measured at side of crystal) 107 10.53Interface shapes at 70 hours, 100 hours and 105 hours from the start of crystal growth, (a) including latent heat, (b) excluding latent heat . . . 108 10.54Interface shapes at 70, 100 and 105 hours for constant conductivity, (a) ksolid > hiq (b) ksoud = kuq, (c) ksoud < klig 110 10.55Interface shapes at 70 hours for CdTe conductivity, varying the value of the liquid conductivity intercept 112 xi 10.56Interface shapes at 70 hours varying the shape of the solid conductivity curve, (a) normal, (b) linear, (c) constant solid conductivity 113 10.57Interface shapes at 70 hours varying the shape of the liquid conductivity line, (a) normal, (b) constant, (c) increasing 114 10.58Interface shapes at 70 hours with the liquid conductivity increased by 7 times to simulate convective fluid flow, (a) static conductivity, (b) 7 times the static conductivity 116 10.59Crystal grown at 2.78 x 10~4 mm/sec with a furnace gradient of 1 de-gree/mm. (a) Isotherm contours after 2.88 x 105 sec of growth, (b) Solid/liquid interface shape and position at the lines indicated at inter-vals of 3.6 x 104 sec of growth, (c) MRSS contours after 2.88 x 105 sec growth, in MPa 120 10.60Dependence of interface curvature on furnace gradient and crucible veloc-ity (V) 123 10.61 Variation on interface curvature with distance from the start of solidifica-tion for different furnace temperature gradients (G) 124 10.62Isotherms for crystals grown at two different gradients with furnace loca-tions at indicated positions, V = 2.78 x 10 - 4 mm/sec (a) 7.02 x 105 sec, G = 0.1 degree/mm (b) 1.98 x 105 sec, G = 1 degree/mm 125 10.63Variation of interface velocity with distance from the start of solidification for the furnace gradients (G) shown. Imposed crucible velocity indicated by CV 126 10.64MRSS* as a function of furnace gradient for the growth velocities (V) shown. 127 10.65MRSSt as a function of furnace gradient for the growth velocities (V) shown. 128 xi i 10.66Crystal grown at 2.78 x 10 - 4 mm/sec with various furnace gradients. The furnace location relative the the crystals are indicated in the Figure, (a) MRSS contours after 2.16 x 105 sec of growth with a 4 dergree/mm gradi-ent, (b) MRSS contours after 2.52 x 105 sec of growth with a 1 dergree/mm gradient, (c) MRSS contours after 3.96 x 105 sec of growth with a 0.5 der-gree/mm gradient 129 10.67Dependence of interface curvature as a function of cold furnace temper-ature. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V = 2.78xl0-4mm/sec 131 10.68MRSS* and MRSSt as a function of cold furnace temperature. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V = 2.78xl0~4mm/sec . . 132 10.69Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm, hot zone temperature of 1380 K, (a) cold zone temperature of 1360 K, (b) cold zone temperature of 1345 K, (c) cold zone temperature of 1335 K 133 10.70Dependence of interface curvature as a function of hot furnace temper-ature. Cold zone temperature is 1335 K, G = 1.0 degrees/mm, V = 2.78 xl0~ 4 mm/sec 135 10.71MRSS* and MRSSt as a function of hot furnace temperature. Cold zone temperature is 1335 K, G = 1.0 degrees/mm, V = 2.78 x l O - 4 mm/sec . 136 10.72Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm, cold zone temperature of 1335 K, (a) hot zone temperature of 1390 K, (b) hot zone temperature of 1380 K, (c) hot zone temperature of 1370 K 137 xm 10.73(a) MRSS contours in MPa (b) Excess stress above the CRSS (MPa). G = 0.5 degree/mm, V = 2.78 x 10~4 mm/sec, Hot Zone = 1370 K, Cold Zone = 1360 K, 3.6 x 105 sec from start of crystal growth. (crcrss = 0.0232 exp (^M)MPa) 139 10.74Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm using 7 x the normal liquid conductivity, (a) Isotherm contours after 2.88 X 105 sec of growth, (b) Solid/liquid interface shape and position at the lines indicated at intervals of 3.6 x 104 sec of growth, (c) MRSS contours after 2.88 x 105 sec growth, in MPa 141 10.75Dependence of interface curvature on furnace gradient and crucible veloc-ity (V) using 7 x the normal liquid conductivity 143 10.76Variation on interface curvature with distance from the start of solidifica-tion for different furnace temperature gradients (G) using 7 x the normal liquid conductivity 144 10.77Variation of interface velocity with distance from the start of solidification for the furnace gradients (G) shown 7 X the normal liquid conductivity used. Imposed crucible velocity indicated by CV 145 10.78MRSS* as a function of furnace gradient for the growth velocities (V) shown. 7 X the normal liquid conductivity 146 10.79Dependence of interface curvature as a function of cold furnace tempera-ture using 7 x the normal liquid conductivity. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V = 2.78x 10~4mm/sec 148 10.80MRSS* as a function of cold furnace temperature using 7 x the normal liquid conductivity. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V = 2.78xl0"4 mm/sec 149 xiv 10.81Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm using 7 x the normal liquid conductivity, hot zone temperature of 1380 K, (a) cold zone temperature of 1360 K, (b) cold zone temperature of 1345 K, (c) cold zone temperature of 1335 K 150 10.82Dependence of interface curvature as a function of hot furnace temperature using 7 x the normal liquid conductivity. Cold zone temperature is 1335 K, G = 1.0 degrees/mm, V = 2.78 xl0~ 4 mm/sec 151 10.83MRSS* as a function of hot furnace temperature using 7 x the normal liq-uid conductivity. Cold zone temperature is 1335 K, G = 1.0 degrees/mm, V = 2.78xlO-4 mm/sec 152 10.84Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm using 7 x the normal liquid conductivity, cold zone temperature of 1335 K, (a) hot zone temperature of 1390 K, (b) hot zone temperature of 1380 K, (c) hot zone temperature of 1370 K 153 10.85(a) MRSS contours in MPa (b) Excess stress above the CRSS (MPa) using 7 x the normal liquid conductivity. G = 0.5 degree/mm, V = 2.78 x 10~4 mm/sec, Hot Zone = 1370 K, Cold Zone = 1360 K, 3.6 x 105 sec from start of crystal growth. ( a c r „ = 0.0232 exp MPa) 155 A.86 Domain for a three-dimensional field problem 168 A.87 Element Shape a) Quadratic isoparametric element b)Quadratic isopara-metric boundary element 174 A.88 Mapping of quadratic element from 7' and z space into £ and r\ space . . 176 A.89 Linear shape function for time 183 A.90 Quadratic shape function for time 186 A.91 Variation of heat capacity and enthalpy through a phase change . . . . : 189 xv A.92 Variation of fraction liquid / with temperature contours (interface shape) 191 A. 93 Boundary elements that are a solidification gap 194 B. 94 Strain and Stress involved in the analysis of axi-symmetric solids 196 B. 95 Quadratic isoparametric element for stress analysis 197 C. 96 The reduced viscosities of liquid metals and their dependence on reduced temperature 207 D. 97 Thermal circuit of Vertical Bridgman boundary conditions 210 D.98 Thermal circuit of Vertical Bridgman boundary conditions with resistances.212 xvi List of Tables 5.1 Etchent to determine the A and B side of the 111 face 23 5.2 Etching Reagents 24 5.3 Procedure for image analysis 25 6.4 Schmidt Factors 28 6.5 Critical Resolved Shear Stress of Cdo.96Zno.04Te 28 9.6 Assumptions used in temperature model 67 9.7 Assumptions used in stress model 67 9.8 Resolved Sheat Stress Components in (001) Crystals [13] 80 9.9 CdTe Properties 81 9.10 Quartz Properties 82 lO.llBoundary conditions and thermal conductivity used for comparing con-duction model to analytical solution 83 10.12F.E.M. and analytical temperatures for a small temperature difference between inside and outside radius (T(r=0.oim) == 0.0/0 Z(r=i.6im) = 1.0K) 88 10.13F.E.M. and analytical temperature profile for a large temperature differ-ence between inside and outside radius (T( r = 0 . 0 1 m) = 0.0/O ZV=i.6im) = 60.0/0 8 9 10.14Boundary conditions and thermo-physical properties used for comparing convection model to analytical solution, f HTC = heat transfer equation 89 xvii 10.15F.E.M. and analytical temperature profile for a small temperature differ-ence between inside and outside radius (T(r=o.oim) — 0.0K, = 1.0K) . 92 10.16F.E.M. and analytical Temperature profile for a large temperature differ-ence between inside and outside radius (T(r=o.oim) = 1.0K, = 100.0A') 93 10.17Boundaiy conditions and thermo-physical properties used for comparing transient model to analytical solution 94 10.18Boundary conditions and thermo-physical properties used for comparing latent heat inclusion to analytical solution 96 10.19CdTe Model Operating Parameters used in the sensitivity analysis . . . . 100 10.20CdTe Model Operating Parameters for examining the effect of furnace gradient and crucible velocity in the interface curvature, interface velocity and MRSS 117 10.21CdTe Model Operating Parameters for examining the effect of hot and cold zone temperatures on the interface curvature and MRSS 118 A.22 Sampling locations and weight coefficients of the Gaussian Quadrature Formula 178 xvm List of Symbols Vz T MRSS CRSS A t [ ] { } [ r 1 r k p cP V T Ni, Nj R Interface velocity in the z direction Shear stress Maximum resolved shear stress Critical resolved shear stress Excess shear stress, MRSS — CRSS Rectangular or square matrix Column, row matrix Matrix transpose Matrix inverse Three dimension isotropic material Surface area Conductivity Density Specific heat Gradient operator Radial direction Axial direction Angular direction Temperature Trial (shape) function Residual xix Wi Weighting function £,?7 Isoparametric coordinates n r , nz normal component to surface in r and z directions hc convective heat transfer coefficient Ambient temperature F Radiation view factor a Stefan-Boltzman constant a Absorptivity of surface hr Radiation heat transfer coefficient [K] — [Kc] + [Kn] Global stiffness matrix [iQ Global conduction matrix Global convection matrix {T} Global temperature matrix {T} Global temperature transient matrix ( [C] Global capacitance matrix {/} = {A} + {/'"} Global force matrix {h} Global convection force matrix {Tr} Global radiation force matrix [Kf Local (elemental) stiffness matrix {TV Local (elemental) temperature matrix {Ty Local (elemental) temperature matrix [cr Local (elemental) capacitance matrix {fV Local (elemental) force matrix Entry i,j of local stiffness matrix X X [J] Jacobian [B] • B matrix A Operator for mapping surface element H enthaply / (T) Fraction liquid L Latent heat of solidification {a} Global nodal displacment vector Ui, Vj displacments of node i {e} Strain vector {a} Stress vector [D] Materials properties matrix E Youngs Modulas v Poissons ratio Lip Potential energy of a structure {e0} , {&o} initial strains and inital stresses {F} Body forces {(f)} surface tractions {P} Loads applied to nodes by external agencies MRSS* Largest MRSS at side of the crystal by the interface MRSSt Largest MRSS at side of the crystal adjacent to were the gradient furnace and cold furnace join Gr Grasliof number H Fluid viscosity j3 Fluid coefficient of thermal expansion g acceleration of gravity xxi D characteristic dimension of body Th Homologous temperature M Molecular weight 8 Atomic radius TM Melting temperature n Particle number density R Gas constant | Energy parameter T* Reduced temperature V* Reduced volume fx* Reduced viscosity X X l l Acknowledgement I would like to thank Dr Fred Weinberg, Dr Indira Samerasakera Dr Carlos Schvezov and Lee Li for their interest, suggestions and help over the course of this project. The assistance of the profesional and technical staff at UBC, in particular, W. Fajber, M. Mager and L. Fredericks is greatly appreciated. I especially thank my parents, Walter and Ruth Parfeniuk, for their love, encourag-ment and understanding that they have given me over the years. I would like to thank my girlfriend, Robin Larson, for making my life complete. Thanks are also extended to my fellow graduate students, especially Derek Riehm, Dave Tripp, Bernardo Hernandez-Morales, Ed Chong and Steve Cockcroft for their friendship and help during this project. I would like to express my gratitude to the British Columbia Science Council for financial support and Johnson Matthey Electronics for material support during my studies. xxin Dedication I dedicate this to my good friends Brian and Doris Williamson. Their courage and strength during his ongoing cancer problems have been an inspiration to my work on this thesis. xxiv Chapter 1 Introduction Cadmium telluride (CdTe) and Cdi-^Zn^Te are II-VI compounds used in 7-ray and X-ray spectrometers and as a substrate for Hg 1_ ; rCd xTe infrared detectors [1,2]. The large band gap of CdTe allows 7-ray and X-ray spectrometers to be run at room temperature as opposed to the silicon or germanium equivalent which must be cooled with liquid nitrogen to avoid excessive thermal currents. The high atomic number of CdTe (z=50) allows a high interaction of photons with the CdTe spectrometer unlike silicon (z=32) or germanium (z=16). The crystal growth methods for the different materials are as different as the materials themselves. Czochralski is the prefered crystal growth method for silicon and germanium while CdTe is grown using a simple liquid/crucible technique. The zinc addition to CdTe is used to strengthen the lattice and is usually added in amounts of less than four atomic percent. The presence of zinc has a negligible effect 011 the thermal diffusivity [9]. 1.1 Crystal Growth in Crucibles Growth of a single crystal from a melt contained in a container is a simple method when compared to Czochralski crystal growth. This technique is also more economical since the growth apparatus is simple and the growth process requires little or no operator supervision. Disadvantages with crystal growth in a crucible arise from the possible nucleation of new grains at the liquid-crucible boundary. The growth of the spurious grains depends on the solid-liquid interface shape (Figure 1.1). If the solid liquid interface 1 Chapter 1. Introduction 2 is convex to the melt the new crystals will grow towards the crucible and will stop after they reach it. If the solid liquid interface is concave to the melt, any new crystals will grow towards the inside of the melt and will dominate crystal growth. The actual control of the interface shape depends on the ratio of radial to axial heat flow. Zero radial heat flow will give a flat solid-liquid interface. If the radial heat loss is greater than the axial heat loss a concave interface will form. For the case were the axial heat loss is greater than the radial heat loss a convex interface will form. 1.2 Vertical Bridgman Growth The Vertical Bridgman method of crystal growth (Figure 1.2) is conducted in a furnace that has four distinct zones: Reservoir zone, Hot zone, Gradient zone and Cold Zone. The crystal is grown in a quartz crucible with the inside surface coated with graphite. The graphite prevents liquid CdTe from reacting and sticking to the quartz. Extra cadmium is present at the top of the quartz crucible. The extra cadmium in the reservoir zone is used to control the over-pressure of cadmium vapor present during crystal growth by adjusting the reservoir zone temperature. The cadmium over-pressure in turn controls the stoichiometry of the CdTe melt [3]. The phase diagram of CdTe is shown in Figure 1.3. The solid phase CdTe region on an enlarged scale is shown in Figure 1.4. For each point on the liquidus there is a corresponding pressure of the volatile component (Cd or Te). To grow at a specific com-position the exact over-pressure of Cd must be maintained during the crystal growth. As solidification takes place segregation will make the melt more concentrated in either Cd or Te depending on the composition of the material being grown. The reaction between the liquid and the gas will push the melt toward the proper stoichometric ratio (equa-tion 1.1). Since the vapor-liquid equilibrium is achieved much faster than the segregation Chapter 1. Introduction 3 1 (a) (b) Figure 1.1: The growth of spurious crystals nucleated at the intersection of the solid-liquid interface with the crucible wall, (a) convex interface, the spurious crystals are grown toward the crucible surface, (b) concave interface, spurious crystals grow in an dominate crystal growth. Chapter 1. Introduction 4 FURNACE T E M P E R A T U R E - LOW TEMPERATURE FURNACE HIGH TEMPERATURE FURNACE CdTe o o MODIFIED BRIDGMAN 'RESERVOIR T E M P E R A T U R E SOLIDIF ICATION T E M P E R A T U R E Figure 1.2: Vertical Bridgman Growth. due to solidification, proper stoichiometry in the melt is maintained. CdTe{liq) ^ Cd(g) + ^Te2(g) (1.1) In general to grow a single crystal the bottom of the crucible is formed into a conical shape or a capillary tip is added to a flat bottom crucible as shown in Figure 1.5. The crucible with the conical bottom leads to supercooling of a small volume of liquid at the tip, as the crucible is lowered through the furnace, at the start of solidification. Only a few crystals nucleate, and it is hoped that one orientation will predominate, resulting in a single crystal having the orientation of the grain which was dominant. The capillary Chapter 1. Introduction Figure 1.3: Phase Diagram of CdTe [1] Chapter 1. Introduction 6 liquid I0>8 I0 1 7 lO 1 6 lO 1 6 lO 1 7 Excess Cd, (atoms/cm3) Excess Te, (atoms/cm3) Figure 1.4: Representation of Solid CdTe Phase Region [1] (a) (b) Figure 1.5: Different Ampule Shapes for Vertical Bridgman Crystal Growth (a) Conical Bottom (b) Capillary tip attached on bottom Chapter 1. Introduction 8 tip serves the same purpose as the tip of the conical bottom. In this case the rapid increase in diameter of the crystal as it leaves the capillary tip could lead to nucleation and growth of other orientations on the flat bottom surface of the crucible. As growth proceeds, in both crucible geometries, nucleation may occur at the crucible walls during growth leading to multicrystal growth. The process of growing a CdTe crystal consists of the following steps. Prior to crystal growth weighed amounts of cadmium and tellurium are placed in a quartz ampoule and vacuum sealed. The crucible containing the metals is then heated to approximately 800°C and held for 24 hours to allow the cadmium and tellurium to react with one another. Following this the crucible is then heated to 1150°C to melt the reacted cadmium telluride. Finally the crucible is lowered through the furnace at a rate between 1.39 x 10~4 and 5.56 x 10 - 4 mm/sec during which time the CdTe solidifies progressively from the bottom of the container. Chapter 2 Literature Review 2.1 Thermal Models All reported mathematical models dealing with Vertical Bridgman growth of CdTe [5,7,9, 11], CdHgTe [8] and GaAs [10] employ the assumption of axi-symmetry about the vertical axis. Without this assumption a three dimensional model would be required which would need more computer memory and longer solution times to get the same solution as the axi-symmetric case. The finite element method was used to solve the mathematical model for most situations [5,7,9,10,11] whilst in one case the finite difference method has been used [8]. The assumption of quasi-steady state has been employed in most of the reported studies. Quasi-steady state ignores the intital transient [7,9] or it is assumed that the velocity of the interface is equal to the furnace velocity [5,10]. The transient can be approximated in the quasi-steady state case by assuming that the interface velocity is equal to the furnace velocity. In this case the transient changes from ^ to Vz^, where Vz is the velocity of the interface in the vertical (z) direction. Solving the heat equation including the transient has been done for the Horizontal Bridgman [6] and the Vertical Bridgman method [11]. The assumption of quasi-steady state with the interface velocity equal to the furnace velocity is not satisfactory since the interface velocity will change according to its location to the end of the crucible. When the interface is at the top of the crucible there will be no thermal mass ahead of it. This 9 Chapter 2. Literature Review 10 will cause the interface velocity to increase as it approches the end of the crucible. The boundary conditions at the crystal wall are convective and radiative heat transfer with the radiative heat transfer dominating. The methods for incorporating surface heat transfer are heat flux [5,8] and radiation using a heat transfer coefficient approximation [7,9,10,11]. The full derivation of the radiation heat transfer approximation is given in Appendix A.5. It has been shown that view factors should be utilized to compute radiation heat flux [8,11] to accurately account for the boundary heat transfer. The radiation that reaches the quartz crucible surface from the furnace is partially absorbed by the quartz, the remainder being transmitted directly to the CdTe [8,10,11]. Natural convection in the liquid CdTe during growth was found to have little impact upon the temperature field [11]. This is consistent with a stable melt which is hot at the top and cold at the bottom. The solid liquid interface was calculated to be concave upward when the conductivity of the solid is less than that of the liquid [5,7,8,10,11]. The interface was found to be convex upward for the case when the conductivity of the solid is less than that of the liquid and the conductivity in the solid and liquid are constant [9]. The interface could be changed from concave to convex by using a lower temperature in the cold zone [8]. However, using a lower temperature in the cold zone moves the isotherms closer together in the solid. The crucible material and crucible wall thickness have a marginal affect on the thermal field [10]. 2.2 Stress Models and Dislocation Studies Stress analysis of the crystal during Vertical Bridgman Growth was carried out assuming that the stress is linearly related to the strain over the range of defomation considered [10, 12]. The stress analysis was then used to determine whether dislocations are generated in the crystal and to estimate the dislocation density. This was accomplished in previous Chapter 2. Literature Review 11 studies using three different methods. In one case plastic deformation and dislocation density was assumed to be directly proportional to the local shear stress T [12]. In a second case the local stress was resolved along the slip plane in the slip directions, and the maximum shear stress (MRSS) was selected. The critical resolved shear stress (CRSS) for the local crystal temperature was subtracted from the MRSS to obtain the shear stress A T associated with the plastic deformation. The dislocation density could then be estimated by assuming the density increases linearly with c ^ s [10]. The third case assumes the dislocation density is proportional to the local plastic deformation thus A T [13]. The three procedures for estimating the dislocation density are listed below: local shear stress, r MRSS - CRSS A T CRSS °r CRSS A T The results of the calculations showed that the maximum stress was near the crys-tal/quartz interface and close to the solid/liquid interface [10,12]. After growth disloca-tion densities in CdTe were observed on etched {111} surfaces of transverse sections and the densities were found to be a factor of 10 to 20 times higher at the sides of the crystal compared to the center [9]. This is consistent with the predictions made by the mathe-matical models. The stress level on the crystal was reduced by using a low temperature Chapter 2. Literature Review 12 gradient during growth combined with a low growth rate [10]. An etch pit investigation of CdTe crystals grown at different gradients [14] indicated that at a high temperature gradient of 25°C/cm the etch pit density is greater than 106 /cm 2 with the etch-pits exhibiting a subgrain structure. When the temperature gradient is less than 10°C/cm the etch pit density dropped to less than 105/cm2 with no subgrain structure. 2.3 CRSS of CdTe There are very few deformation studies of CdTe at high temperatures. Gutmanas et al. [15] measured the CRSS of CdTe in the range of -173°C to 227°C. The CdTe sample size used was 4x3x8 mm3 which were cut and mechanically polished to size. The stress strain curves were obtained for compression, using a cross head speed of 50 /xm/min. The results, giving the CRSS as a function of temperature in the range of —173°C to 227°C are shown in Figure 2.6. Hall and Vander Sande [16] examined the plastic deformation of CdTe up to a tem-perature of 500°C. The sample size was 2.5 x 2.5 x 5 mm3. The samples where formed using a wire string saw, then mechanically and finally chemically polished. As with Gut-manas the samples were tested in compression giving the results shown in Figure 2.7. The stress/strain curves show a strong decrease in work hardening with increasing tem-perature with little work hardening at 500°C, where CdTe becomes ductile. In both sets of observations the maximum temperature is low, at 500°C. The tem-perature range of interest in crystal growth is 800°C to 1080°C which is well above the range examined by previous investigations. As a result the CRSS above 800°C can only be estimated by extrapolation of the data in Figure 2.6. Chapter 2. Literature Review 13 Figure 2.6: Critical resolved shear stress of CdTe [15] Chapter 2. Literature Review 14 Figure 2.7: Stress/strain curves of CdTe single crystals of identical orientations deformed at the temperatures indicated. Orientation of compression axis of sample is shown in stereographic triangle [16] Chapter 3 Objectives The objective of this research programme was to investigate dislocation formation during crystal growth of Cdo.96Zno.04Te in a Vertical Bridgman furnace. In order to calculate the operating conditions that create dislocations the critical resolve shear stress of CdTe must be known. Any mathematical model used must be verified using experimental analysis. Taking the preceding into consideration results in the following methodology: 1. A tensile test apparatus was built to measure the critical resolved shear stress (CRSS) of Cdo.96Z1io.04Te as a function of temperature in the range of 460°C to 870°C to determine the stress required for dislocation formation and multiplication. 2. Interface profile experiments were carried out to verify the heat transfer model. 3. Etch pit densities studies were carried out on Cdo.96Zno.04Te slices to verifj' the thermal stress model. 4. A finite element temperature model was written to determine the thermal history of the system during crystal growth. 5. A finite element model was written to calculate the thermal stresses. The maximum resolved shear stresses (MRSS) were calculated from the stress components and compared to the CRSS. The dislocation density was then determined from the excess stress A T . 15 Chapter 3. Objectives 16 6. The influence of the furnace gradient and crucible velocity on the interface shape, interface movement and MRSS were examined. The influence of hot and cold zone temperature on the interface shape and MRSS were examined. Chapter 4 Critical Resolved Shear Stress Experimental Procedure Small plates of Cdo.96Zno.04Te of various orientations were obtained from Johnson Matthey Electronic Materials Ltd. The zinc content of the samples was 4 atomic percent. The orientation of each plate was determined by the X-Ray back reflection Laue technique. The plates were then shaped into small tensile specimens and deformed in tension at temperatures between 460°C and 870°C. 4.1 Sample Preparation 4.1.1 Polishing The samples received from Johnson Matthey were 22.2 x 7.94 x 1.9mm3 rectangular single crystals wafers with unpolished surfaces. One surface of the sample was mechanical-chemically polished using a 90% ethanol - 10% bromine solution for chemical polishing and lintless cloth for mechanical polishing. The process consisted of wrapping a lintless cloth around a piece of glass. The glass wrapped cloth was placed into a glass dish to contain the chemical polishing solution. The surface was polished by applying the polish solution to the cloth and moving the sample in circular motion across the cloth. During polishing very little pressure was applied to avoid breaking the sample. Extra polishing solution was added as the cloth dried out. After polishing the surface of the sample was very smooth with an occasional ripple. 17 Chapter 4. Critical Resolved Shear Stress Experimental Procedure 18 4.1.2 Crystal Orientation Measurements Back reflection Laue X-Ray pictures of the polished surface were obtained using a Phillips X-Ray diffraction machine with Fe radiation [40]. The Cdo.96Zno.04Te were placed on cardboard using two sided tape, the polished surface was face up. The cardboard was taped to the sample holder on the X-ray diffraction machine. The sample surface was positioned 3 cm from the film holder and irradiated for 20 minutes at 30 kV and 10 mA. After irradiation, the X-Ray film was developed and the orientation determined using a Greninger chart and Wulff net as described in reference [40]. 4.1.3 Shaping the Tensile Samples The Cdo.96Zno.04Te samples were shaped into tensile specimens having the configuration shown in Figure 4.8 A. The samples fitted into the test sample holder for tensile de-formation shown in Figure 4.8 B. This configuration of samples and holder allowed the assembly components to expand when heated to test temperature without straining the sample. The samples were shaped by grinding using the procedure described below. The process requires care since Cdo.96Zno.04Te is very brittle at room temperature and breaks readily. Two glass slides were cut to the same dimensions of a sample and heated on a hot plate. Lakeside resin was melted onto one surface of each slide. A sandwich assembly was then made with the Cdo.96Zno.04Te samples in the centre and the glass slides on the outside joined by resin. The assembly was then cooled to room temperature. The assembly was subsequently mounted vertically on a grinding machine, and the upper edge of the assembly progressively ground using a 51 mm diameter fast speed diamond impregnated grinding wheel. The assembly was moved horizontally as the edge was ground in each pass with vertical steps between passes of about 0.6 mm. This process Chapter 4. Critical Resolved Shear Stress Experimental Procedure 19 was continued until the final circular contour was reached as shown in Figure 4.8 A. The assembly was then turned over and the opposite side ground in the same way. After grinding the assembly was heated to melt the resin, and then separated. Wax from the samples was removed with an ethyl alcohol solution. 4.2 Tensile Testing Apparatus The tensile testing apparatus shown in Figure 4.8 was made from a nickel super-alloy in order to withstand the high temperatures that are required for the tests. A piece of nickel super-alloy was used on each grip (B) to prevent the samples from falling but of the grips prior to the test. In the tensile test the load axis is parallel to the longitudinal axis of the sample but is not coincident with it. The seperation between the axes is 2 mm. This will give rise to eccentric loading which could affect the measured stress. The gripping system has considerable lateral movement which can accomodate any eccentric loading. Accordingly the small eccentric loading should not significantly affect the stress values measured. A chromel-alumel thermocouple was welded to each grip to measure the furnace temperature during the tensile test. Molybdenum rods were used to connect the tensile grips to the load cell and cross head on the Instron. The instron tester used a type E or CT load cell with a cross head speed of 0.85 x 10 - 3 cm/s for all tests. A resistance tube furnace was used to heat the sample. The furnace temperature was recorded on a chart recorder and a microprocessor furnace controller was used for temperature control. One of the thermocouples that was attached to the tensile grip was also used as the furnace control thermocouple. Helium gas was used in the furnace to prevent the Cdo .96Z1io.04Te from disassociating at the test temperatures. Chapter 4. Critical Resolved Shear Stress Experimental Procedure 20 4.3 Tensile Testing Procedure The chart recorder was calibrated using a voltage generator before each test. Calibration of the Instron was done using a dead weight. The sample was placed between the tensile grips and the furnace was put over the assembly. The gas flow rate was 30 cc/min for all tests. The furnace was turned on and set to the temperature set point. During heat up the crosshead was moved apart to compensate for thermal expansion of the molybdenum rods. After the furnace had reached steady state the sample was deformed with a cross head speed of 0.85 x 10~3 cm/s. As the sample was loaded, the crystal self-aligned in the grips, then deformed elastically and plastically. Chapter 4. Critical Resolved Shear Stress Experimental Procedure 21 Figure 4.8: Tensile sample(A) and tensile test apparatus(B). Chapter 5 Etch Pit Experimental Procedure Slices of Cdi_ xZn xTe that had been cut along the {111} plane were obtained from John-son Matthey Electronic Materials Ltd. The samples were polished and etched at Johnson Matthey. The polishing and etching procedure used is standard industrial practice to determine the etch pit density on the (111) plane [18]. 5.1 Etching The structure is oriented such that the {111} planes are horizontal. From Figure 5.9 it is evident that the top (111) surface, represented by the open circles is different from the bottom surface represented by the closed circles. We define the top (111) surface (A) as made of Te atoms and the bottom surface (B) ( i l l ) as Cd atoms. For delineating dislocations by etch pits the B surface is used, the etch preferentially attacking Cd atoms. For a given crystal wafer, the B side is delineated by using an orientation etch. The etch consists of the etching solutions and etching times listed in Table 5.1. The sample does not have to be polished prior to etching. The orientation etch produces a bright shiny finish on the A surface and a dull finish on the B surface. Having determined the B surface of the wafer, the sample is polished on both sides with a 5% bromine-rnethanol solution The polished samples were etched in the solution developed by Nakagawa [18] given in Table 5.2. The etching procedure consisted of immersing the sample in the Nakagawa etch for 60 seconds followed by immersion in 1:1 mixture of reagent E and water for 1 22 Chapter 5. Etch Pit Experimental Procedure 23 (III) SURFACE till") SURFACE Figure 5.9: Crystal structure of CdTe (Inoue et al [19]) Step Components Amount Time 1 H 2 0 2 parts H 2 0 2 1 part 2 min etch HC1 1 part 2 HC1 1 part H 2 0 2 1 part 30 sec quench 3 H 2 0 30 sec wash Table 5.1: Etchent to determine the A and B side of the 111 face Chapter 5. Etch Pit Experimental Procedure 24 Name Components Amount Nakagawa etch H 2 0 2 parts H 2 0 2 2 parts HF 3 parts Reagent E HN0 3 10 ml H 2 0 20 ml K 2 C r 2 0 7 2 grams Table 5.2: Etching Reagents second and washing in water for 30 seconds. The composition of reagent E is given in Table 5.2. 5.2 Image Analysis A Leitz image analyzer was used to map the dislocation etch pit density and distribution over the entire surface of a CdTe sample. The programs used for CdTe were those which were developed for mapping GaAs samples. The image analyser consists of an optical microscope with a videocamera detecting and recording the optical image. An area of the etched sample is observed under the microscope using an 8x objective lens, and the image is divided into a range of grey levels using a microprocessor attached to the videocamera. The area in the image (162 fxm xl62 jim) having a grey level above a set level is then measured and recorded in a file. The sample is then automatical!}' translated to an adjacent area and the process repeated. The camera automatical!)' focusses on the surface features when the sample is translated. To convert the total dark area measured to etch pit density, the average dark area of an etch pit must be determined. This is done by selecting six areas on the sample at random, measuring the grey area above the fixed level with the image analyser, then manually counting the etch pits in the area to determine the density. This establishes the effective gray area of the average etch pit on Chapter 5. Etch Pit Experimental Procedure 25 Step Task 1 Set lighting and detection level so that only areas that are dislocations are detected. 2 Image analyzer collects data. 3 Manually detects area associated with dislocations and count dislocations to get an area per dislocation ration. Repeat this at 6 random locations. Table 5.3: Procedure for image analysis the sample. The average gray area of an etch pit is then used to determined the etch pit density for each area scanned and the results plotted to produce an etch pit density map over the entire surface. The steps in mapping the etch pit density are summarized in Table 5.3. In measuring the etch pit density, a number of factors can markedly change the observed density. A stain or dark film in an area of the sample surface can markedly change the gray levels detected in the area. This will indicate an erroneous sharp increase in the etch pit density in the area. If the entire surface is not a {111} plane, the stray crystal surface will not etch and give zero etch pit density in this area. At the edges of the sample, the surface is rounded, which results in smaller pits as the surface deviates from a {111} plane giving lower etch pit densities at the edges. Surface stains can generally be detected visually and removed or allowed for in the etch pit density image. Stray crystal surfaces without etch pits are readily identified by their size and shape in the image. If the etch pit density at the edges falls off due to rounding, the density can be determined by manualy counting the pits in the region. Chapter 6 Critical Resolved Shear Stress Experimental Results 6.1 Crystal Orientations The orientation stress axis relative to the samples was plotted on a Wulff net in Fig-ure 6.10. All test samples were oriented such that the large face was parallel to {111} and the [T01] direction was inclined to the stress axis. Table 6.5 gives the angle from the [T01] direction to the force axis. Using the sample orientations the Schmidt factors for the samples were calculated to determine which slip plane and slip direction acted during crystal yield (Table 6.4). 6.2 Tensile Tests The stress-displacement curve for sample 4, tested at 932 K is shown in Figure 6.11. The stress was determined from the measured load and the minimum cross-sectional area of the test sample (Table 6.5). The displacement refers to the movement of the crosshead during the test. The initial portion of the stress-displacement curve is observed to be linear to near 5 MPa stress, which corresponds to the elastic deformation. The slope of the curve then decreases significantly above 5 MPa indicating work hardening. The yield is determined by drawing the best lines through the elastic and initial plastic portion of the curve, as indicated by the dashed lines and taking the point of intersection as the 3'ield stress. The CRSS is then the resolved yield stress on the operating slip plane and direction. 26 Chapter 6. Critical Resolved Shear Stress Experimental Results 27 Figure 6.10: Wulf net for 111 plane in Cdo.9eZno.04Te with force direction for each tensile test Chapter 6. Critical Resolved Shear Stress Experimental Results Orientation Sample Plane Direction 1 2 3 4 5 6 111 TlO 0 0 0 0 0 0 OTl 0 0 0 0 0 0 101 0 0 0 0 0 0 T i l TOT 0.3448 -0.0627 0 -0.0734 0 0 OTl -0.1914 0.4624 0.5104 0.1710 0.4045 0.4045 110 -0.1598 0.2342 -0.4019 -0.0920 -0.4045 -0.4045 T i l Oil -0.3287 -0.3834 -0.4176 -0.0923 -0.4096 -0.4096 TOT -0.2202 0.0762 0.3731 0 0 0 1T0 0.5590 0.2991 0.4292 -0.2875 0.4096 0.4096 1T1 Oil -0.3519 0.0610 -0.0277 0.0739 0 0 TOT -0.2357 -0.0121 -0.2988 0 0 0 TTO -0.1093 0.0758 -0.0254 0.3747 0 0 Table 6.4: Schmidt Factors Sample Temperature Orientation Dimensions (cm) CRSS (K) Planef Angle to p l ] | length width (MPa) 1 738 (111) 42 2.92 1.91 4.78 2 764 (111) 10 2.95 1.30 3.88 3 817 (111) 3 2.84 1.27 5.22 4 932 (111) 49 2.82 1.27 1.72 5 1024 (111) 0 2.90 1.32 0.98 6 1136.5 (111) 0 2.95 1.30 0.86 Table 6.5: Critical Resolved Shear Stress of Cdo.9eZno.04Te. f Plane parallel to the surface of the sample JAngle between the stress axis and the [TOl] direction Chapter 6. Critical Resolved Shear Stress Experimental Results 29 0.35 Cross Head Displacement (cm) Figure 6.11: Stress-Displacement curve for a tensile test of Cdo.96Zno.04Te at 932K. The stress-displacement curves for tests 1 to 6*are shown in Figure 6.12. The curves for 738, 764 and 817K are very similar showing steep elastic curves and extensive work hardening in the plastic region. The yield points are indicated by the diamonds. At 932 K and above, the slope of both the elastic and plastic portions of the curve are shallower than in the lower temperature curves. The oscillations present in the curve at 1136 K are associated with the temperature controller, and thus not considered significant. Chapter 6. Critical Resolved Shear Stress Experimental Results 16 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Cross Head Displacement (cm) Figure 6.12: Stress Displacement Curves for tensile tests of Cdo.9eZno.04Te Chapter 6. Critical Resolved Shear Stress Experimental Results 31 6.3 Critical Resolved Shear Stress The values of the CRSS for samples 1 to 6, determined from the yield stress are listed in Table 6.5. The measured CRSS values are plotted on a logarithmic plot against the reciprocal of the test temperature, in Figure 6.13. The best fit line for the six data points is shown, and corresponds to the following equation: Reported values of the CRSS [15] have been extrapolated above 500 K giving the dashed line shown in Figure 6.13. The line is appreciably lower than the present measured values below about 1000 K and higher above 1000 K. The CRSS decreases with increasing temperature for all measured values except sample 3 (817 K). The high value is attributed to the conducting of the test before the sample had reached temperature equilibrium. The value of CRSS extrapolated to the melting point is 0.43 MPa which is much lower than the CRSS of GaAs at its melting point (approximately 0.9 MPa). o-, crss = 0.0232 (in MPa) Chapter 6. Critical Resolved Shear Stress Experimental Results 32 Temperature (K) 1414 1314 1214 1114 1014 10 914 814 714 0 0 1 CO Pi u Melting Point Best Fit - - -J Previous Study 0.1 10 11 10,000/T (1/K) 12 13 14 Figure 6.13: Dependence of critical resolved shear stress on temperature for Cdo.96Zno.04Te. Chapter 7 Etch Pit Density and Distribution Etch, pit density maps were done on four samples; slices 1 and 11 from boule 521 (521s01, 521sll) and slices 11 and 28 from boule 555 (555sll , 555s28). Each slice was cut at a 15 degree tilt from the axial direction. Figure 7.14 shows the relative location of the slices in the boule. The large faces of the slices are (111) planes in the crystal. The EPD map for slice 521s01 is shown in Figure 7.15. The pit density is indicated by the colours, the darker colours corresponding to higher densities, as given in the figure. Regions 1 and 2 are white, which would correspond to an etch pit density less than 105 c m - 2 in these regions. However this is not the case; area 1 is a stray crystal in which the etched face is not (111) and therefore etch pits do not form. The white strip, area 2, is an artifact resulting from the image analyser focusing on the bottom of the etch pits instead of the top surface. The flat surface, 3, of the sample is the bottom of the crystal which was adjacent to the bottom surface of the crucible. The remaining curved edges, 4, are intersections of the cut surface with the vertical outer walls of the boule. The EPD on the (111) surface of the sample ranges between 1 x 105 to 2.2 x 105 cm - 2 , the lowest density occuring in the center of the sample, and the highest density at the outer curved edge. The appearance of the etched surface at three positions in the sample shown in Figure 7.16 are given in Figure 7.17a, b, c. In the three areas (A - C) in Figure 7.17a, b, c, the etch pits are the same size and form a cellular structure, in which many etch pits are concentrated along cell walls with the remainder being randomly distributed within the cells. The etch pit density at the centre (area A, Figure 7.17a) is appreciable 33 Chapter 7. Etch Pit Density and Distribution 34 lower than that at the outside edge (area C, Figure 7.17c) which is consistent with the image analysis results of 1.2 x 105 c m - 2 at the center and 2.2 x 105 c m - 2 at the edge respectively. The pits aligned along straight lines in area B (Figure 7.17b) are considered to be artifacts resulting from residual polishing lines on the polished surface. The EPD of sample 521sll is shown in Figure 7.18 in which edges 1 and 2 are intersections of the sectioning plane with the center surface of the boule, edge 3 is a vertical cut through the centre of the slice, edge 4 is at the bottom of the boule. The slice has a low etch pit density at the centre and higher density at edge 4. Moving from the centre to sides 1 and 2 the density increases and at edge 1 decreases. The etched surface was photographed at six areas (A-F) in the sample at the positions shown in Figure 7.19 The photographs were taken to determine the distribution of etch pits in a given area, and to qualitatively compare the observed pit density to the image analysis results. The results for the six areas are shown in Figure 7.20a-f. As in sample 521s01 all six areas show a strong cellular pattern in which the pits are concentrated along the cell walls. There is appreciable variation in the etch pit density, cell wall thickness and density of pits in the cell walls, cell diameter, and pit density inside the cells in the areas examined. The lowest dislocation density was observed in area A (Figure 7.20a) at the centre of the boule with a density of 1.4 x 105 cm - 2 . The cells are large, approximately 0.2 mm in diameter, and most of the pits are in the cell walls. The etch pit distribution at the edge of the sample, areas B , C and D are shown in Figure 7.20b, c, and d respectively. The etch pit density is clearly larger than that of region A, estimated as twice the density, with smaller etch pits, and multiple etch pit lines along cell walls. The light-region along the edge of the sample in the image analysis results is an artifact due to rounding of the edges during polishing. The appearance of the etch pit density in area E is shown in Figure 7.20e. The density is clearly higher in this area than the others examined, with much less pronounced cell structure and wide Chapter 7. Etch Pit Density and Distribution 35 cell boundaries. The line of higher etch pit density, area F, is shown in Figure 7.20f. This is similar to area E, but with larger pits and lower density than E consistent with the image analysis results. The image analysis results of the etch pit density of sample 555sll is shown in Fig-ure 7.21. Side 1 is the top of the crystal, sides 2 and 3 are the intersection of the (111) surface with the vertical walls of the crystal, and side 4 is a cut edge through the cen-tre of the crystal. High etch pit densities were measured adjacent to sides 1 and 3 and halfway up side 2. A thin band of higher density is also observed along side 4. Six areas A-F (Figures 7.23a to f) in the sample were photograhed, their locations being shown in Figure 7.22. As with all previous observations the six areas show a cellular pattern with the etch pits concentrated at the cell walls. There is a large variation in the etch pit density, and thickness of the cell walls. Areas A, B and C (Figures 7.23a - c) show a low etch pit density of 1.5 x 105 cm - 2 . The cell size is approximatley 0.2 mm for area A, B and C. There are some scratches at area C (Figure 7.23c) due to polishing. The scratches were not dark enough to be detected as etch pits. The etch pit distribution at the edge of the sample, areas D, E and F are shown in Figure 7.23d to f. The etch pit density present in areas D, E and F are twice that seen in areas A, B and C. The EPD map for slice 555s28 is shown in Figure 7.24. The three large white regions 1, 2 and 3 are stray crystals and therefore do not etch. The flat surface of the sample, 4, is the top of the crystal. The remaining curved edges, 5 and 6, are intersections of the cut surface with the vertical outer walls of the boule. The white streaks adjacent to area 3 and edge 6 are due to the image analysis focusing on the bottom of the pits instead of the top surface. For areas A-D (Figure 7.26a to 7.26d) were examine as shown in Figure 7.25. Area A has the lowest dislocation density in the sample at 1.5 x 105 cm - 2 . The etch pits are in a cellular structure with the cells having a diameter of 0.2 mm. The region of the crystal adjacent to crucible are shown in areas B and C (Figures 7.26b and Chapter 7. Etch Pit Density and Distribution 36 c). The etch pit density is approximatley double that found in area A (2.5 x 105 cm - 2). The diameter of the cells are less than 0.05 mm. Area D (Figure 7.26d) is located at a region close to a stray crystal. The etch pit density appears to be higher than that found in area B and C. Only a small portion of area D has a cellular structure. The remainder of the surface consists of a high density of etch pits relative uniformly distributed on the surface. Chapter 7. Etch Pit Density and Distribution slice # 28 Figure 7.14: Location of slices in boules 521 and 555 Figure 7.15: Etch pit density map for sample 521s01, Mag x 3 Chapter 7. Etch Pit Density and Distribution Figure 7.16: Areas of sample shown in Figures 7.17a, to 7.17c Chapter 7. Etch Pit Density and Distribution •10 Figure 7.17a: Etch surface of sample 521s01 at area A, Mag x 80 Figure 7.17b: Etch surface of sample 521s01 at area B, Mag x 80 Chapter 7. Etch Pit Density and Distribution 41 Figure 7.17c: Etch surface of sample 521s01 at area C , Mag x 80 Figure 7.18: Etch pit density map for sample 521sll, Mag x 3 Figure 7.19: Areas of sample shown in Figures 7.20a to 7.20f Chapter 7. Etch Pit Density and Distribution •1-1 Figure 7.20a: Etch surface of sample 521sll at area A, Mag x 80 Figure 7.20b: Etch surface of sample 521sll at area B, Mag x 80 Chapter 7. Etch Pit Density and Distribution 45 Figure 7.20d: Etch surface of sample 521sll at area D, Mag x 80 Figure 7.20f: Etch surface of sample 521sll at area F, Mag x 80 Chapter 7. Etch Pit Density and Distribution 47 Figure 7.21: Etch pit density map for sample 555sll, Mag x 3 Figure 7.22: Areas of sample shown in Figures 7.23a to 7.23f Chapter 7. Etch Pit Density and Distribution 49 Figure 7.23b: Etch surface of sample 555sll at area B, Mag x 80 Chapter 7. Etch Pit Density and Distribution 50 Figure 7.23c: Etch surface of sample 555sll at area C, Mag x 80 Figure 7.23d: Etch surface of sample 555sll at area D, Mag x 80 Chapter 7. Etch Pit Density and Distribution Figure 7.23f: Etch surface of sample 555sll at area F, Mag x 80 Chapter 7. Etch Pit Density and Distribution 52 Figure 7.24: Etch pit density map for sample 555s28, Mag x 3 Chapter 7. Etch Pit Density and Distribution slice # 28 t • i H i f ' Figure 7.25: Areas of sample shown in Figures 7.26a to 7.26d Chapter 7. Etch Pit Density and Distribution 54 Figure 7.26a: Etch surface of sample 555s28 at area A, Mag x 80 Figure 7.26b: Etch surface of sample 555s28 at area B, Mag x 80 Figure 7.26d: Etch surface of sample 555s28 at area D Mag x 80 Chapter 8 Interface Position and Shape The heat transfer mathematical model predicts the position and shape of the solid/liquid interface during solidification. Accordingly, if these interface characteristics were mea-sured, they could be used to verify the model predictions. A number of techniques to determine interface position and shape were considered. The experiments were carried out by Johnson Mattheys and some limitations applied. The interface position might be determined from thermocouple measurments of CdTe during solidification. This was not done, as the thermocouples could not be brought out of the quartz container without increasing the possibility of the quartz exploding during growth due to pressure from the Cd vapour. Radioactive tracers were considered, to mark the interface position but this was not adopted due to contamination problems if the quartz failed during growth. Attempts were made to mark the interface shape and position using the Zn that is added to the CdTe, and interrupting growth at given times during solidification. Two crystals of CdZnTe were grown, containing 4 atomic percent Zn, and the growth rate interrupted by a sharp drop in the ciystal position of 2 mm, or a sharp change in growth velocity from 1.0 mm/hr to 3 mm/hr. It was hoped that these interruptions would trap a con-centration or depletion laj'er of Zn at the interface which would then be delineated by either etching, SEM/WDS, or SIMS. It was not known if Zn was segregated with k > 1 or k < 1. since a pseudo-binary phase diagram of the CdTe-Zn system could not be found. Two crystals were grown in order to determine if the profile measurements were repeatable. Both crystal were grown using the same experimental conditions. 56 Chapter 8. Interface Position and Shape 57 After growth the CdZnTe crystal was sectioned on a plane containing the vertical crystal axis, and then was etched using nitric acid. The nitric acid preferentialy removes Te atoms. The crystal was not a single crystal, as shown in Figure 8.27(a). A schematic diagram of the surface is shown in Figure 8.27(b) delineating the grain boundaries and the Areas 1, 2 and 3 examined in detail. The etched surface exhibited a line, readily visible with the naked eye, marked A, B, C, D in Figure 8.27(b) which could correspond to the instantaneous increase in growth velocity. The line is flat in the center of the crystal and steeply inclined to the horizontal near the side walls. There was also a band across the sample near the top surface marked E, F, G. in Figure 8.27(b), as shown in Figure 8.28. This is associated with the 2 mm drop near the end of crystal growth. The shape of the interface is difficult to see since the the side that is between E and F appears to be concave while the side between F and G appears to be convex. The second crystal that was grown in this set of experiments had an accidental furnace shut down near the end of crystal growth. The resulting crystal is shown in Figure 8.29(a). The result of shutting the furnace off was a clear marking of the interface near the end of solidification. Figure 8.29(b) is a closeup of the interface region when the furnace shut off. The interface is clearly convex at the top of the crystal. No line was present that was similar to A, B, C, D in the crystal shown in Figure 8.27 If the line ABCD in figure 8.27(b) delineates the interface as a result of a change in growth rate, then there should be a large change in Zn concentration, to produce the preferential etching along the line. This was investigated in Areas 1, 2 and 3 in Figure 8.27(b). Initial attempts were made by scanning across the line using the WDS analysis system of an SEM. The Zn signal was found to be too low to obtain significant composition measurements. SIMS observations were then made across the line, since SIMS has appreciably higher sensitivity than WDS. In general the results from the SIMS observations were inconclusive. For Area 1 five Figure 8.27: Cdo.96Zno.04Te boule grown to examine interface profile Chapter 8. Interface Position and Shape 59 Figure 8.28: Top of boule were multicrystalline growth occured Chapter 8. Interface Position and Shape 61 scans across the line were made (Figure 8.30). In all cases the signal from the side toward the crystal centre was lower than the outer side indicating lower electrical conductivity. The signal was enhanced by flooding the surface with oxygen. In three of the scans there was a significant drop in signal of 6 4 Z n + in a region about 20 fxm in width at the demarcation line. The scans for 1 1 4 C d + and 1 3 0 T e + were constant across the line. This suggests a drop in Zn concentration at the line. In the other two scans (Figure 8.31 and Figure 8.32) all three elements showed a sharp drop at the line. In the other scans all three elements showed a sharp drop at the line followed by a slow rise, suggesting the changes are due to surface anomalies. In one scan for Area 2, a depression at the demarcation line for 6 4 Z n + was observed, but when the scan was repeated adjacent to the initial scan with an e~ flood on the surface the depression disappeared. In Area 3 two scans were made across the demarcation line with an e~ flood of the surface. In both cases there were no significant increases of decreases of the 6 4 Z n + signal at the demarcation line. In summary, the results indicate that there may be a 6 4 Z n + depletion at the demarca-tion line but this remains only a possibility, since the depletion is observed in one crystal and not the other. Thus these results will not be compared to model prediction since the comparison would not be conclusive. Chapter 8. Interface Position and Shape 62 Scan Zn Cd Te Inside ! Outside Zn Cd Te Zn Zn Cd Te Zn Figure 8.30: SIMS analysis for Zn, Cd and Te at Area 1 Chapter 8. Interface Position and Shape 63 Figure 8.31: SIMS analysis for Zn, Cd and Te at Area 2 Figure 8.32: SIMS analysis for Zn, Cd and Te at Area 3 Chapter 9 Mathematical Models 9.1 Scope of Models and Assumptions Models were developed to predict the temperature and stress that occur during crystal growth. The temperature distribution in the solid and liquid CdTe during solidifica-tion was calculated with the model. During solidification the CdTe passes from a high temperature furnace, through a temperature gradient region where solidification occurs, and finally to a cold zone furnace where the solid CdTe cools. Cooling from the cold zone furnace to room temperature was not considered. A model to determine the stress distribution in the solid CdTe after solidification was developed, based on the calculated temperature distribution in the crystal. A number of assumptions were used to simplify the model. These are listed in Table 9.6 for the thermal model and Table 9.7 for the stress model. The tables also include comments on the assumptions and their validity. A major assumption made is that the furnace, crucible and Cdi_ xZn xTe are axisymmetric. In general Vertical Bridgman growth systems are axisymmetric. Derivations due to non symmetric heating elements or crucibles which are not centered make the production of single crystals more difficult and therefore are corrected. Radiative heat transfer between the furnace and the crucible, is an important bound-ary condition in the model. Radiation heat transfer also occurs between the crucible and the Cdi_ xZn xTe when the metal seperates from the crucible due to shrinkage during cooling. The extent of the heat transfer in both cases depends on the orientation of 65 Chapter 9. Mathematical Models 66 the thermal source and the absorbing surface. If the source and absorbing surfaces are parallel, a maximum amount of heat is transferred between the surfaces. If the surfaces are not parallel radiation view factors are calculated to estimate the heat transfer in this case. For any surface that is parallel to a radiative source, it is assumed that the temperature directly adjacent to the source will dominate, and radiation from portions of the emitting surface not adjacent to the source will be lower due to the low view factor. The Cdi_ xZn xTe was assumed to deform elastically in order to convert the calculated thermal strains to local stresses in the growing crystal. The mechanical properties of CdTe/Cdi_ xZn xTe are non isotropic. Youngs modulus and Poissons ratio are dependent on the crystal orientation of the crystal relative to the deformation directions [33,34]. We note that in many cases large grained polycrystals are grown instead of single crystals. This will make the calculations less rigorous. The stress model assumes an axisymmetric stress distribution. The alternative of a three dimensional solution would require ex-cessive computing time and memory. The axisymmetric assumption will give the same general location of the stresses and the variation across the crystal as the three dimen-sional model. The inside of the quartz crucible containing the Cdi_ xZn xTe is coated with graphite. This prevents the metal from reacting and sticking to the quartz surface. As the crystal cools it separates from the quartz container due to greater volume contration. The graphite ensures that separation occurs easily without local deformation of the crystal. 9.2 Idealized Domain The crystal growth domain adopted for the mathematical model is shown in Figure 9.33. It is noted that the model assumptions make the domain a simplified version of actual growth conditions. For example, the furnace temperatures do not change in the model as Chapter 9. Mathematical Models 67 Assumption Validity Comment Temperature Axisymmetric Excellent Normal practice Transient included Excellent Best approximation for actual heat transfer View factors for surfaces perpendicular to furnace Excellent Angle of absorbing surface to radiation source has large effect on heat transfer No view factors for surfaces parallel to furnace Good Furnace temperature directly adjacent to surface will be rate controlling Table 9.6: Assumptions used in temperature model Assumption Validity Comment Linear Elastic Stresses Fair Gives a good estimate Isotropic Mechanical Properties Unknown Mechanical properties are dependent on crystal orientation. The assumption gives a reasonable estimate Single crystal Poor Often large grained polycrystals grown No CdTe sticking to walls Excellent Industry standard Table 9.7: Assumptions used in stress model Chapter 9. Mathematical Models 68 the mass of the crucible and metal move through the thermal field. The top end of the quartz crucible is assumed to be an insulated boundary, which is not entirely the case. 9.3 Temperature Model The full derivations of the heat transfer differential equations are given in Appendix A. The heat transfer equation was solved in cylinder coordinates using Galerkins method. The simplified heat transfer equation is: where r and z are the radial coordinates; T is the temperature; k is the thermal conduc-tivity; p is the density; Cp is the specific heat; and t is time. The mesh is moving with the domain thus there is no bulk flow term in the heat transfer equation. Solving gives the matrix equation: where [K] is the thermal stiffness matrix, [C] the capacitance matrix, {/} the thermal force matrix, {T} the vector containing nodal temperatures, and {T} is the transient temperature vector. The thermal stiffness matrix consists of a conduction term ([/Cj]), a convection term ([/<"/,]), and a radiation term ([/C])-The thermal force matrix consists of a convective term ({fh}) and radiative term ({/,}). (9.2) [Al {T} + [C] {T} = {/} (9.3) [K] = [Kc] + [Kh] + [Kr] {/} = {fh} + {fr} The individual elements of the matrices and the vectors are: Chapter 9. Mathematical Models 69 f u r n a c e t e m p e r a t u r e Figure 9.33: Model representation of domain Chapter 9. Mathematical Models 70 Khi. = Jrrhc NiNjdT KTij = I rYT NiNjdT dj = 2TT / / r p C p NiN3dr dz flu = J rhNiT^dT fn = J NiT^dT were Ni,Nj are the shape (Trial) functions; h is the convective heat transfer coefficient; Too the ambient temperature; hr the radiation heat transfer coefficient; T a surface area; and dT the surface integral. Appendix A contains the full derivation of the heat transfer model. The type of shape function (element type) chosen has important implication on the accuracy of the solution. A quadratic isoparametric element was employed to give in-creased accuracy for the stress analysis (Appendix A.3). The quadratic isoparametric element consists of eight nodes with four sides. The sides of the element can be curved due to the quadratic nature of the shape functions. A quadratic isoparametric element is shown in Figure 9.34. The surface integrals are applied over one or more sides of the element. Numerical integration is required since the integration in equation 9.3 cannot be done exactly. The numerical integration method used is Gaussian Quadrature. Gaussian Quadrature involves evaluating the formula using a number of sampling points. Nine sam-pling points were used for the temperature model since the thermophysical properties are temperature dependant. The more sampling points allow the variation in conductivity, density, etc to be better modelled over the domain. Solving equation 9.3 including the transient is generally difficult analytically and cannot be done for non-linear situations. Thus the solution technique that is used is a Chapter 9. Mathematical Models 71 Figure 9.34: Eight noded quadratic isoparametric element Chapter 9. Mathematical Models 72 trial function (finite element) discretization in the time domain. This method is known as a recurrence technique since each time step is solved using the previous time step temperatures as initial conditions. The number of trial functions used, will affect the accuracy and stability of the solution. In general the greater the number of shape (trial) functions used the better the accuracy and stability [28]. A two point (linear) scheme was used for the initial solution with a three point (quadratic) scheme used for all subsequent calculations. 9.3.1 Boundary Conditions Two of the boundary conditions in the thermal model deal with radiation and forced convection. The dominant mode of heat transfer during crystal growth between the furnace and the crucible is radiation. For non parallel surfaces a view factor is required to determine the heat transfer, as is the case for the horizontal top and bottom of surfaces of the crucible. For these surfaces the view factor for an annular ring radiating on a cylindrical band is used [31]. With a radiation view factor an assumption is required as to how the top and bottom surfaces see the furnace temperature profile. The top and bottom surfaces see the furnace temperature profile changing as a step function, as shown in Figure 9.35(B). instead of linearly. Thus for any boundary element on the top or bottom of the crucible the radiation view factor for each temperature step in the gradient furnace has to be calculated. For the sides of the crucible it was assumed that the temperature changes linearly (Figure 9.35(A)). A side element will have a radiation view factor of one, and only see the furnace temperature that is directly parallel to it. During radiative heat transfer some radiation will be absorbed and some transmitted through the body. Quartz has a low absorptivity, thus most of the radiation is transmitted to the CdTe. Thus there are two radiative boundaries, the surface of the quartz where a fraction of the radiation is absorbed and the surface of the CdTe where the remainder of Chapter 9. Mathematical Models 73 Hot Furnace Temperature Figure 9.35: Furnace Temperature as seen by side of crucible (A) and as seen by the bottom and top of the crucible (B) the radiation is absorbed. The outside surface of the CdTe is assumed to be a black bod}r thus absorbing all the incoming radiation. All radiative boundary surfaces are assumed to have zero reflectivity. It is assumed that a gas is forced through the furnace during crystal growth. The temperature is assumed to be 10°C below the furnace temperature. A forced convective heat transfer coefficient is used to approximate the forced convection [32]. The dominant mechanism for heat transfer at the high temperatures present during crystal growth is the radiation that is transmitted through the quartz. A calculation Chapter 9. Mathematical Models 74 of the thermal resistances for heat transfer at the boundary under typical conditions is given in Appendix D. 9.3.2 Phase Transformation Solution In a heat transfer model of solidification a means must be found to accurately introduce the latent heat of solidification. If the latent heat is not properly applied during any time step, the resulting temperature fields will be incorrect. A simple variation of the enthalpy method [21,23,24] is used to incorporate the latent heat (Appendix A.8). The enthalpy of a material can be expressed as the sum of the specific and latent heats. The latent heat evolved is a function of the local fraction of liquid [25]. For this method the enthalpy is written as H= !J(Cp + f(T)*L)dT . (9.4) JTT where / (T) =fraction liquid present between solidus and liquidus L =latent heat of solidification Equation 9.4 can be rewritten to separate out the specific heat term. H = [TCpdT+ fT f (T) * LdT JTT JTT The heat transfer equation can now be rewritten as: k d ( 8T\ d2T dT d rT Applying Galerkins method with the latent heat included, gives an extra term that is included into the final matrix form of the equation: [Kl]=2TJJ rp^LN.Njdrdz Chapter 9. Mathematical Models 75 CdTe has a univariant solidification temperature. In order to use this method of including latent heat a solidus/liquidus gap of 0.5 degrees is assumed. The fraction liquid between solidus and liquidus (/) is assumed to be linear with temperature. 9.3.3 Solidification Gap Since the volume shrinkage of Cdi_ xZn xTe on cooling is larger than quartz, by one order of magnitude between 600 K and 1364 K, the gap will form between the metal and quartz on cooling. This will change the heat flux between the quartz and the metal. Before solidification there is good contact between the melt and the crucible, with conduction as the mode of heat transport. After solidification and gap formation the mode of heat exchange changes from conduction to radiation. The solidification gap is modelled as if it were a radiative boundary, with the ambient temperature given by the temperature of the adjacent boundary (Figure 9.36). The estimated shrinkage gap is small (/xm), as a result nodes at each side of the gap occupy the same r-z space. The convective matrix and vector for the right side boundary are: The value for T/ e / t is the temperature of the left side boundary. It is evaluated at the left side elements gaussian integration points that are'directly opposite to the right side element gaussian integration points. The left side of the boundary is evaluated as the right side is: Before solidification the value of the heat transfer coefficient is inflated to approximate conduction. A value of 8000.0 J / m 2 s K for the heat transfer coefficient across the gap Chapter 9. Mathematical Models 76 CdTe Visualization Figure 9.36: Boundary elements that are a solidification gap before solidification was found to approximate conduction [30]. Between the liquidus and solidus temperature the inflated heat transfer coefficient is ramped down to the normal value. The normal value of the heat transfer coefficient for radiation heat transfer is temperature dependent. The formulais equation A.19 in Appendix A. Using this method the quartz and CdTe have two separate temperature models which are linked to each other by the solidifaction gap boundary. Chapter 9. Mathematical Models 77 9.4 Stress Model The stresses and strains that occur in cylindrical coordinates are shown in Figure 9.37. The stress/strain calculations are based on the displacement that occurs in the domain that is being modelled. Using the same quadratic isoparametric element as used in the temperature analysis there are eight nodes, each node describe the displacment in two principle directions r and z. The general matrix equation as formulated in Appendix B is: [K] {«} - {/} = 0 where [K] is the stiffness matrix: [K] = 2*JJr [B]T [D] [B] dr dz and {/} is the force vector: {f} = 2irJJj[B]T [D]{e0}drdz The [B] matrix consists of the differential operators and [D] is the material properties matrix (see Appendix B). The nodal displacment vector is {a} and {eo} is the thermal strain vector which is described as: {*} = f 1 < > = < coo 0 V. J where a is the thermal expansion coefficient and A T is the temperature change that is causing the thermal stress. The temperature changes are calculated at each Gaussian Chapter 9. Mathematical Models Figure 9.37: Strain and Stress involved in the analysis of axi-symmetric solids Chapter 9. Mathematical Models 79 integration point. The relationship between displacment and strain is: e2 {'} = { = [B] {a} (9.5) The relationship between stress and strain is given by: {*} = { were {eo} is the thermal (initial) strain, {<70} is the initial stress. 9.4.1 Critical Resolved Shear Stress The local stresses in the crystal were determined from the thermal strains. The stresses were then resolved on the 12 slip systems in CdTe {111} <110>. This was done using the procedure developed for GaAs [13], CdTe and GaAs having the same structure. The maximum resolved shear stress (MRSS) was selected from the operating slip system and compared to the CRSS for CdTe at the temperature being considered. The excess stress A T = MRSS — CRSS was then assumed to generate and multiply dislocations locally, the number of dislocations roughly increasing with A T . The formulas to resolve the stresses are given in Table 9.8. 9.4.2 Boundary Conditions The boundaries of the crystal are considered to be free surfaces for the stress model. The inside and outside of the quartz crucible are also assumed to be free surfaces. This Chapter 9. Mathematical Models 80 Plane Direction Resolved Shear Stress Mode (in) TlO OTl i o f ^ {- (crr - ae) cos (29) + y/2 rTZ cos (9 + f) } & {- (ar - a6) y/2 sin9 sin + f) + (az - ae) + rrz cos9} {(o-r ~ °~e) V2 sin9 sin [B + f ) + (az - ae) - T t z cos9) la Ha Ilia (Til) TlO OTl 10T ^ {(ar - ag) y/2 cos9 cos (0 + f) - (ar - az) - rrz sin9J & {- (crr - ae) \/2 sin9 cos (& + \ ) + (ar - az) - rrz cos9} ^ {- (ar - ae) cos (29) + y/2 rrz cos (9 + f) } IVa Vb la (Til) TlO OTl 10T ^ {- (ar - ae) \/2 sin9 sin [9 + f) + (az - a0) - T t z cos9J & {(ar - a6) y/2 sin9 sin (o + \ ) + (az - ae) + rrz cos9) 4{-(ar- a0) cos (29) + y/2 rTZ cos (6 + *)} lib Ilia la (111) TlO OTl i o T ^ {- (°"r - <T$) y/2 sin9 cos (9 + f) + (ar - az) + rrz cos9J &{(ar-a0)y/2 cos9 cos ((9 + \ ) - (ar - az) + rrz sin9} *§{-(o-r~ ae) cos (29) - y/2 rrz cos [9 + f ) } Va IVb lb Table 9.8: Resolved Sheat Stress Components in (001) Crystals [13] assumption that the inside surface is free arises from the liquid being unable to apply a pressure to the quartz and the solid CdTe shrinking faster than the quartz. The surface of the CdTe is also assumed to be a free surface. In order to avoid any free body movement [20] a node at the top of the quartz and CdTe are pinned. The value of Youngs Modulus is lowered if the CdTe is liquid. This is to approximate the sudden drop in stress at the solid/liquid interface. If the material is liquid Youngs Modulus is set to an arbitrary value of 1 MPa. Between the liquidus and solidus Youngs Modulus is linearly ramped up to the normal value. 9.5 Thermo-Physical Properties An extensive literature search was done to find the thermo-physical properties of CdTe and quartz. Table 9.9 gives the CdTe thermo-physical properties selected and used in Chapter 9. Mathematical Models 81 Property Units Temperature Polynomial Conductivity W / K m T < 500K 10.701-1.3 x i o - 2 r 500K < T < 1200K 5.337 - 2.263 x 10~3T 1200K < T < 1364K -38.31 + 6.78 x 10~'2r -2.811 x 10~ 5T 2 T > 1364K 9.987 - 5.676 x 10~3T Heat of Fusion J / g T = 1364 209.2 Specific Heat J / g K 298K < T < 1600K 0.1778+ 1.663 x 10~4T -6.607 x 10- 8 T 2 +1.703 x 1 0 - n T 3 density g / m 3 - 6.2 x 106 HTC of Gap f J / m 2 s K T > 1364K 8000.0 Forced Conv. HTC J / m 2 s K - 50.0 Poissons Ratio - - 0.459 Young Modulus MPa - 3.98 x 105 Expansion Coefficient T < 525K 1.031 x l O " 6 + 2.054 x 10- 8T -19.779 x 10- 1 2 T 2 T > 525K 6.3629 x 10~6 Table 9.9: CdTe Properties. jHTC = heat transfer coefficient the model and Table 9.10 the quartz thermo-physical properties. The low temperature thermal conductivity (T < 300K) of CdTe is well established [1]; the high temperature values (T > 1173A") have been measured experimentally [9]. The values of the thermal conductivity between 300K and 1173K were approximated by a straight line connecting the known values. The heat capacity of CdTe is well known between 300K and 1364K [36]. Above the melting temperature the heat capacity was assumed to belong to the same curve as the lower temperature. The thermal expansion coefficient of CdTe has been measured up to 400K above which it is constant [37]. Poissons ratio within the {111} plane is used for all stress calculations [33,34]. Chapter 9. Mathematical Models 82 Property Units Temperature Polynomial Conductivity W/K m 298K < T < 1600K 1.25 + 1.05 X 10~3T -3.95 x 10~7T2 +7.33 x 10- n T 3 Specific Heat J / g K 298K < T < 1600K 0.760 + 6.34 x 10~4T -2.39 x 10~7T2 +4.44 x 10- n T 3 density g / m3 - 2.2 x 106 Quartz Absorptivity - - 0.3 Poissons Ratio - - 0.16 Young Modulus MPa - 7.17 x 104 Expansion Coefficient T < 1250K 4.028 x 10-7 + 5.466 x 10-1 0T -4.623 x 10- 3T 2 T > 1250K 3.637 x 10-'' Table 9.10: Quartz Properties. fHTC = heat transfer coefficient Chapter 10 Modelling Results 10.1 Comparisons to Analytical Solutions 10.1.1 Conduction Considering only steady state heat conduction, the governing equation can be solved according to fixed temperature or zero flux boundary conditions giving: [IQ {T} = 0 (10.6) The finite element solution is compared to the one dimensional analytical solution for fixed temperatures at an inside radius of 0.01m and outside radius of 1.61 m. Two cases were examined, one for a low temperature difference between the inside and outside radius and the second for a high temperature difference between the inside and outside radius. The domain used for the calculations is shown in Figure 10.38. The boundary conditions are given in Table 10.11. The calculations in the two cases were made to ensure that the percentage error in Units Small A T Large A T Tr=o.01m K 0.0 0.0 T r = 1.61m K 1.0 60.0 k W/ m K 1.0 1.0 Table 10.11: Boundary conditions and thermal conductivity used for comparing conduc-tion model to analytical solution 83 Chapter 10. Modelling Results z A 0.0250 --0.0125 0.0 dT dT dz = 0 0.01 0.21 0.41 0.61 0.81 1.01 1.21 1.41 1.61 Figure 10.38: Domain used for testing the conduction matrix Chapter 10. Modelling Results 85 the F.E.M. solution did not change with a sharp increase in temperature at the bound-ary. Calculated temperature profiles for the low temperature difference is shown in Figure 10.39 and the large temperature difference in Figure 10.40. The associated errors of the temperature profiles are given in Tables 10.12 and 10.13. The percentage error for the two cases are the same for the first four significant figures. The error in the F.E.M. formulation increases towards the inside of the cylinder up to a maximum of 15%. The errors are relatively low (< 5%) between 7- = 0.31 m and r = 1.61 m. Close to the z axis (r = 0.01 m to r = 0.31 m) the errors increase rapidly. The high error close to the z axis is due to the rapidly changing temperature gradient, which the quadratic F.E.M. solution is unable to fit accurately. For all cases the error that is associated with the region close to the z axis can be reduced by decreasing the mesh size locally if the temperature in the radial direction has a large gradient. 10.1.2 Convection As in the conduction matrix the temperatures from the current time step are calculated at the Gauss points and are used for the calculation to solve the following equation: [ [IQ + [Kh] ] {T} = {fh} (10.7) Equation 10.7 describes combined conductive and convective heat transfer. The finite element solution is compared to the one dimensional analytical solution for fixed tem-peratures at the inside boundary and a fixed ambient temperature of the outside radius. The domain used for the test is shown in Figure 10.38. The boundary conditions and physical properties are given in Table 10.14. The resulting calculations for low A T is shown in Figure 10.41 and a large A T is shown in Figure 10.42, the associated errors are given on Tables 10.15 and 10.16. The temperatures that are predicted along the convective boundary are within 1.5% of the Chapter 10. Modelling Results 0 . 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Radial position (m) Figure 10.39: Temperature profile for a small temperature difference between inside outside radius (T( r = 0 .oim) = 0-OK, T^=imm) = 1.0K) Chapter 10. Modelling Results 87 Figure 10.40: Temperature profile for a large temperature difference between inside and outside radius (T( r =o.oim) = 0-0K, T(r-i.6im) — 60.OA") Chapter 10. Modelling Results 88 radial position Temperature (K) % Error (m) F.E.M. Analytical 0.01 0 -2.0E-17 0 0.11 0.401507 0.471896 14.91624 0.21 0.552072 0.599149 7.857434 0.31 0.637445 0.675794 5.674785 0.41 0.699102 0.730816 4.339544 0.51 0.747082 0.773767 3.448776 0.61 0.786494 0.809003 2.782372 0.71 0.819882 0.838878 2.264488 0.81 0.848876 0.864809 1.842482 0.91 0.874482 0.887719 1.491131 1.01 0.897422 0.908237 1.190791 1.11 0.91819 0.926816 0.930783 1.21 0.937167 0.943792 0.701991 1.31 0.954636 0.959419 0.498560 1.41 0.970819 0.973896 0.315957 1.51 0.985893 0.987380 0.150656 1.61 1 " 1 0 Table 10.12: F.E.M. and analytical temperatures for a small temperature difference be-tween inside and outside radius (T( r = 0 .oim) = 0.0K, T ( r = 1 6 1 m ) = 1.0K) Chapter 10. Modelling Results 89 radial position Temperature (K) % Error (m) F.E.M. Analytical 0.01 0 3.3E-15 0 0.11 24.09041 28.31377 14.91625 0.21 33.12432 35.94898 7.857429 0.31 38.24667 40.54769 5.674851 0.41 41.94614 43.84896 4.339478 0.51 44.82492 46.42605 3.448761 0.61 47.18963 48.54020 2.782376 0.71 49.19289 50.33269 2.264536 0.81 50.93256 51.88859 1.842466 0.91 52.46894 53.26314 1.491077 1.01 53.84529 54.49423 1.190844 1.11 55.09137 55.60899 0.930828 1.21 56.23003 56.62754 0.701960 1.31 57.27813 57.56515 0.498612 1.41 58.24915 58.43376 0.315928 1.51 59.15356 59.24283 0.150679 1.61 60 60 0 Table 10.13: F.E.M. and analytical temperature profile for a large temperature difference between inside and outside radius (T( r = 0 .oim) = O.Oiff, T ( r = 1 6 1 m ) = 60.OA') Units Low A T Large A T Tr=0.01m K 0.0 0.0 T K 1.0 100.0 Conductivity W/ K m 1.0 1.0 Convective HTCf J / m 2 s K 1.0 1.0 Table 10.14: Boundary conditions and thermo-physical properties used for comparing convection model to analytical solution. fHTC = heat transfer equation Chapter 10. Modelling Results 90 Figure 10.41: Temperature profile for a small temperature difference between inside and outside radius (T( r = 0 .oim) = 0.0A", T^ = 1.0A") analytical solution for both cases. The errors within the domain are of the same order of magnitude as in the conduction formulation. 10.1.3 Transient A routine was written to find the transient response to the heat conduction equation (equation A.35). The transient F.E.M. formulation was tested using the domain shown in Figure 10.38 with the boundary conditions and material properties given in Table 10.17. The results of the finite element solution and the corresponding analytical solution [29] Chapter 10. Modelling Results 91 Figure 10.42: Temperature profile for a large temperature difference between inside and outside radius (T( r=0.oim) — 0.0A", — 100.OA') Chapter 10. Modelling Results 92 radial position Temperature (K) % Error (m) F.E.M. Analytical 0.01 0 0 0 0.11 0.35324 0.420497 -15.9947 0.21 0.485705 0.533890 -9.02535 0.31 0.560814 0.602187 -6.87052 0.41 0.61506 0.651215 -5.55203 0.51 0.657272 0.689488 -4.67257 0.61 0.691946 0.720886 -4.01462 0.71 0.72132 0.747507 -3.50334 0.81 0.746829 0.770614 -3.08662 0.91 0.769357 0.791028 -2.73969 1.01 0.789538 0.809312 -2.44332 1.11 0.80781 0.825867 -2.18653 1.21 0.824506 0.840994 -1.96060 1.31 0.839874 0.854919 -1.75986 1.41 0.854112 0.867819 -1.57952 1.51 0.867374 0.879835 -1.41630 1.61 0.879785 0.891080 -1.26757 Table 10.15: F.E.M. and analytical temperature profile for a small temperature difference between inside and outside radius (T(r=0.oim) = 0-0K, Too — 1.0K) Chapter 10. Modelling Results 93 radial position Temperature (K) % Error (m) F.E.M. Analytical 0.01 1 1 0 0.11 35.97074 42.62923 -15.6195 0.21 49.08477 53.85515 -8.85781 0.31 56.52060 60.61655 -6.75714 0.41 61.89093 65.47035 -5.46723 0.51 66.0699 69.25939 -4.60514 0.61 69.50262 72.36780 -3.95919 0.71 72.41064 75.00327 -3.45668 0.81 74.93602 77.29088 -3.04674 0.91 77.16631 79.31185 -2.70519 1.01 79.16427 81.12190 -2.41319 1.11 80.97314 82.76092 -2.16017 1.21 82.62607 84.25846 -1.93735 1.31 84.14753 85.63702 -1.73930 1.41 85.55712 86.91412 -1.56131 1.51 86.87000 88.10368 -1.40025 1.61 88.09872 89.21693 -1.25335 Table 10.16: F.E.M. and analytical Temperature profile for a large temperature difference between inside and outside radius (T( r = 0.oiT O) = 1-0K, = 100.OA') Chapter 10. Modelling Results 94 Units Value Initial Temperature K 0.0 Tr=1.6m K 100.0 Conductivity W / m K 1.0 P g / m3 1.0 cP J / g K 100.0 At sec 2.0 Table 10.17: Boundary conditions and thermo-physical properties used for comparing transient model to analytical solution are shown in Figure 10.43 for the transient response at r = 1.4m. There is an initial oscillation in the FEM solution between 0 and 8 seconds. At times greater than 8 seconds the time stepping increment used gives an approximation that is very close to the analytical solution. 10.1.4 Phase Change The fraction liquid will be assumed to be linear with decreasing temperature between the solidus and liquidus. Thus the fraction liquid (/) varies according to: T Tgoiidug f = (Tliqt idus <T <T. solidus Tliquidus Tsolidus f — 1 (T > Tuquidus) f — 0 (T < TsoHjus) The thermo-physical properties used in the comparision are given in Table 10.18. The domain used is given in Figure 10.44. A sudden and large temperature change was used at the surface to examine how the model would respond to such a change. The initial temperature of the body was 1800 K. The temperature at z = 0 m was set to 1700 K at time > 0 seconds. The sides of the cylinder were considered to be insulated. The thermal history was examined at three locations: z = 0.5 m, z = 1.0 m Chapter 10. Modelling Results 0 2 4 6 8 , 10 12 14 16 18 20 Time (sec) Figure 10.43: Transient response at r = 1.4 m for F.E.M. formulation Chapter 10. Modelling Results 96 Units Value Initial Temperature K 1800.0 Liquidus Temperature K 1800.0 Solidus Temperature K 1799.0 Tr=1.6m K 1700.0 Conductivity W/ m K 1.0 P g/m3 1.0 cP J/gK 100.0 Latent heat J/g 100.0 A* sec 2.0 Table 10.18: Boundary conditions and thermo-physical properties used for comparing latent heat inclusion to analytical solution and z = 2.0 m. The analytical solution [29] and finite element solution are shown in Figure 10.45. Close to the quenching interface (z = 0.5 m) the calculated temperature drops faster than the analytical solution. As time increases the calculated temperature approaches the analytical solution. The deviation from the analytical solution is due to the rapid quench that is applied at the boundary. The calculated finite element solution is unable to include the latent heat in the elements close to the interface since the element changes from being higher than the solidification temperature to less than the solidification temperature in one time step. At z = 1.0m the FEM solution initially rises before following the analytical solution. The temperature rise (oscillation) occurs from the sudden decrease in temperature just in front its location. The quadratic temperature shape function can only fit a certain range of temperature profiles. The nodes at the boundary are instantly put to the quenching temperature causing the quadratic shape function to fit what is essentially a step function in the temperature domain. At z = 2.0m the analytical and FEM solution have the same thermal history. Chapter 10. Modelling Results 97 1 2 5 . 4 1 2 5 . 2 1.2 1.0 0 . 8 0.6 0 . 4 0 . 2 0 . 0 I 1 1 0 . 0 0 . 0 2 5 0 0 . 0 1 2 5 Figure 10.44: Domain used for testing the latent heat solution Chapter 10. Modelling Results Figure 10.45: Comparison to analytical solution for a phase change Chapter 10. Modelling Results 99 00 0.00015 0.0001 0.00005 -0.00005 -0.0001 --0.00015 — Stress r. Analytical — Stress theta, Analytical " • Stress z, Analytical - O Stress r, F.E.M. • Stress theta, F.E.M. - A Stress z, F.E.M. JOT A ^ r ^ ^ ^ .A , A ' A - 1 1 l l i l l 0.001 0.002 0.003 0.004 0.005 0.006 Radial Distance ( m ) 0.007 0.008 Figure 10.46: Comparison between finite element solution and analytical solution for thermal stress in a cylinder. Temperature variation: T = —r2K, Constraint Fz = 0 10.1.5 Stress/Strain The thermal stress formulation was compared to an analytical solution [35] for an imposed radial tempurature gradient of —r2K with zero axial force (Fz = 0). Figure 10.46 gives the comparison between the analytical and finite element formulation. The mesh used in the analysis is given in Figure 10.38. The finite element solution is very close to the analytical solution except for the stress in the z direction close to the centre of the cylinder. Chapter 10. Modelling Results 100 Parameter Units Values Ampoule diameter mm 50.0 Charge height mm 90.0 Quartz wall thickness mm 2.75 Hot zone temperature Kelvin 1380.0 Furnace gradient degrees/mm 1.0 Cold furnace temperature Kelvin 1335.0 Ampoule velocity mm/sec 2.78 x 10~4 Table 10.19: CdTe Model Operating Parameters used in the sensitivity analysis 10.2 Sensitivity Analysis A sensitivity analysis was conducted to determine the optimum operating parameters. The parameters examined are mesh size, time step size, latent heat inclusion, solid/liquid thermal conductvity, and specific heat. The model operating conditions used in the sensitivity analysis are given in Table 10.19 10.2.1 Mesh Size Three different mesh sizes were used to examine the sensitivity of the model to mesh size (Figure 10.47). Figure 10.48 show the axial interface positions at the center of the crystal for the three mesh densities. The interface positions coincide for most of the crystal for the three mesh sizes, at the top of the crystal the interface position moves faster for finer mesh densities. At the start of solidification the interface position exhibits oscillations for all three mesh densities. The mesh size of 2.5mm by 3.0mm was used (Figure 10.47(b)) since it gave a solution close to the higher mesh density but with a reduced computational time. Chapter 10. Modelling Results 101 1 x x x x X \ \ \ \ \ \ k \ . A i l l A A \ \ \ X X \ A.„...\ (a) (b) (c) Figure 10.47: The three mesh densities used to examine the models sensitivity. (a)Fine mesh density, size approximatley 1.25mm by 1.5mm. (b) Normal mesh density, size approximatley 2.5mm by 3.0mm. (c) Coarse mesh density, size approximatley 5.0mm by 6.0mm. Chapter 10. Modelling Results 102 0.1 20 30 40 50 60 70 80 90 100 110 Time (hours) Figure 10.48: Axial interface location during crystal growth for different mesh densities (Interface as measured at center of crystal) Chapter 10. Modelling Results 103 10.2.2 Time Step Size Five different time step sizes (150 sec, 300 sec, 600 sec, 1200 sec, 2400 sec) were used to determine the optimum time step size for the calculations. The results are shown in Figure 10.49 in which it is observed that the step size has little effect on the inter-face position except for 2400 sec which is slightly lower than the others. The interface movement at the side of the crystal was identical for all time step sizes except for the v. largest. At the side of the crystal, Figure 10.50, the difference between 2400 sec and the others is increased. Some oscillations of the interface position is observed in Figure 10.49. Oscillations occur at most time steps except for 1200 second and 150 second time step size. The oscillations are not a function of time step since they do not lessen with smaller time step size. A time step size of 600 seconds was used to reduce computational time while giving the same solutions as the smaller time steps. 10.2.3 Latent Heat The effect of latent heat was examined. The oscillation which occurs at the start of crystal growth are associated with the latent heat component in the model. If the latent heat is not included the oscillations disappear or are markedly reduced. This is shown in Figure 10.51 in which the interface position is shown with and without the latent heat component in the model at the crystal centre. The corresponding interface position at the side wall of the crystal is Figure 10.52. The overall affect of the latent heat is seen by comparing the interface position at the center and side of the crystal. The latent heat reduces the velocity of the the interface at the center on the crystal but does not change the interface location at the side of the crystal. This thus requires the interface to change in shape during crystal growth. The calculated interface shape after 70, 100 and 105 hours of growth is shown in Figure 10.53 with and without the latent heat term in Chapter 10. Modelling Results 104 0.1 10 20 30 40 50 60 70 80 90 100 110 Time (hours) Figure 10.49: Axial interface location during crystal growth for different time step sizes (Interface as measured at center of crystal) Chapter 10. Modelling Results 105 Figure 10.50: Axial interface location during crystal growth for different time step sizes (Interface as measured at side of crystal) Chapter 10. Modelling Results 106 0.1 10 20 30 40 50 60 70 80 90 100 110 Radial position (m) Figure 10.51: Axial interface location for crystal growth including latent heat and ex-cluding latent heat (Interface as measured at center of crystal) the model. The effect of the latent heat is to slow the growth at the center appreciably, and slow the growth slightly at the wall. The growth rate at the center is very markedly retarded in the top 3rd of the crystal due to latent heat. The interface velocity increases at the top of the crystal due to loss of thermal mass ahead of the interface. The release of latent heat at the interface helps to replace the loss of thermal mass thus slowing the interface movement. Chapter 10. Modelling Results 107 0.1 0 1— 1 1 1 1 1 1 i i i 10 20 30 40 50 60 70 80 90 100 110 Radial position (m) Figure 10.52: Axial interface location for crystal growth including latent heat and ex-cluding latent heat (Interface as measured at side of crystal) Chapter 10. Modelling Results Figure 10.53: Interface shapes at 70 hours, 100 hours and 105 hours from the start crystal growth, (a) including latent heat, (b) excluding latent heat Chapter 10. Modelling Results 109 10.2.4 Conductivity The effect of changing the thermal conductivity on the thermal history was examined. Using constant values for the liquid and solid thermal conductivities three cases were examined:(l) Arso;,d > knq, (2) fcso/,d = knq and (3) A:so{,d < kuq. The results are shown in Figure 10.54 in which the interface shapes are drawn after 70, 100 and 105 hours of growth. At 70 hours the interface is slightly concave facing the liquid when the solid conductivity is less than the liquid conductivity; slightly convex when the conductivites are equal; and more convex when ks < k[. At 100 and 105 hours the results are similar to that observed in a previous investigation [5], in which the interface was reported flat when ks = ki and more concave and convex than in that shown in Figure 10.54 for ks greater and less respectively than ki. The location and shape of the interface is controlled by the amount of heat that reaches it and the amount that is removed. The amount of heat supplied to the interface is determined by a number of variables which are the hot zone temperature, furnace gradient, liquid conductivity, liquid heat capacity and latent heat of solidification. By adjusting one of the variables, such as liquid conductivity in this case, the curvature of the interface can be changed. The increase in liquid conductivity increases the supply of heat to the interface thus moving it down. Decreasing the liquid conductivity to less than the solid conductivity adjusts the heat flow to the interface to be less than the removal of heat thus moving the interface upwards and giving it a convex shape. The same can be said about the removal of heat from the interface and the cold furnace temperature and solid conductivity. The value of the liquid conductivity and solid conductivity have great importance in the interface shape. The values to calculate the conductivity polyomials used in this shidy are based on 2 experimental measurements in the liquid and 5 experimental measurements Chapter 10. Modelling Results 110 Figure 10.54: Interface shapes at 70, 100 and 105 hours for constant conductivity, (a) Kolid > hiq (b) ksolid = hiq, (c) ksolid < kiq Chapter 10. Modelling Results 111 in the solid [9]. The possibility of the real values varing from the experimental values is large. Assuming that the measured conductivity may be 7 percent different from the actual values justifies examing the variation of the interface shape with variations in the thermal conductivity. The effect of changing the value of the liquid conductivity intercept with the liquidus temperature was examined. Four cases were considered: (1) kuquid > ksoud at liquidus temperature (Tmp) which is the normal case; (2) kliquid = ksolid at T m p ; (3) kHquid < ksoHd at T m p ; and (4) kuquid >> ksond at Tmp. The actual values of the solid conductivity just below the melting point is 1.87 W/m K. The equations for the liquid conductivities examined are: 1- kiiquid = 9.987 - 5.676 x 1CT3 x T = 2.24 W / m K 2. kuqutd = 9.613 - 5.676 x 10"3 x T = 1.87 W / m K 3. knquid = 9.239 - 5.676 x 10"3 x T = 1.49 W / m K 4. k[iquid = 10.735 - 5.676 x 10~3 x T = 2.99 W / m K The interface shape changes considerably between each case. By decreasing the liquid conductivity (case 1, 2 and 3) the supply of heat to the interface is lowered thus moving the interface up and making it more convex. By increasing the liquid conductivity (case 4) the supply of heat to the interface is increased, moving the interface down and giving it a concave shape. The interface changes from convex (case 1 - normal values) to concave (case 4) by changing the value of the liquid conductivity by 7 percent. The effect of changing the temperature dependence of the solid conductivity was ex-amined. Three cases were considered: (1) measured dependence, (2) linear variation of solid conductivity with temperature and (3) solid conductivity independent of temper-ature. The interface shape after 70 hours and the temperature dependence of the solid Chapter 10. Modelling Results 112 Figure 10.55: Interface shapes at 70 hours for CdTe conductivity, varying the value of the liquid conductivity intercept Chapter 10. Modelling Results 113 (a),(b),(c) I 1 10 mm Figure 10.56: Interface shapes at 70 hours varying the shape of the solid conductivity curve, (a) normal, (b) linear, (c) constant solid conductivity conductivity for the three cases considered is shown in Figure 10.56. Changing the tem-perature dependence of the solid conductivity has little effect on the interface shape or position. Varying the slope of the liquid conductivity temperature was examined. Three cases were considered:(1) Normal properties, (2) constant liquid conductivity and (3) increasing linear liquid conductivity. There is very little change in the interface shape and position for the different cases. The resulting interface position after 70 hours and a schematic of the temperature dependence of conductivity for the three cases are shown in Figure 10.57. Chapter 10. Modelling Results 114 0 0 , (b), (c) Figure 10.57: Interface stapes at 70 hours varying the shape of the liquid conductivity line, (a) normal, (b) constant, (c) increasing Chapter 10. Modelling Results 115 Adjusting the slope of the solid conductivity (Figure 10.57) or the slope of the liq-uid conductivity (Figure 10.57) relative to the melting temperature does not effect the amount of heat transfer away from the interface. In all cases examined the conductivities in the solid and liquid just below and above the melting temperature are the same. It is likely that fluid flow will be occuring in the melt during solidification. This will increase the effective thermal conductivity in the melt. In a tin system [38] multiplying the static conductivity by 9.6 was found reasonable to represent actual heat flow when a temperature gradient of 4 degrees was present across the melt. The Grashof numbers (Gr) of tin and CdTe were calculated to determine if the two systems had similar natural convective behaviour. The calculations are given in Appendix C. The two Grashof numbers are: • Grttn = 6.16 x 104 • GrcdTe = 5.88 x 104 Since both Grashof numbers are approximately the same, increasing the static thermal conductivity to include the effect of fluid flow on heat transfer is valid. The static conductivity will be increased by a factor of 7, a conservative estimate, to account for fluid flow in the melt. The calculated interface shapes with the static liquid conductivity after 70 hours of growth and the modified conductivity are shown in Figure 10.58(a) and (b) respectively. The increased liquid conductivity slows down the interface movement by 15 mm at the center of the crystal and 7 mm at the side of the crystal. The shape of the interface is convex facing the liquid using the normal liquid conductivity and concave using the inflated liquid conductivity. Increasing the bulk heat transfer to approximate fluid flow has the greatest effect on the interface position and shape. Since this is the case two Chapter 10. Modelling Results 116 •(a)-Figure 10.58: Interface shapes at 70 hours with the liquid conductivity increased by 7 times to simulate convective fluid flow, (a) static conductivity, (b) 7 times the static conductivity sets of calculations will be done, one assuming no fluid flow occurs and another using the modified thermal conductivity to account for fluid flow. 10.3 Variables Examined and Operating Ranges The model was employed to explore the effect of furnace gradient, hot and cold zone furnace temperatures and crucible velocity on the interface velocity, interface curvature and MRSS during crystal growth. These variables were examined to determine the Chapter 10. Modelling Results 117 Parameter Units Values Ampoule diameter mm 50.0 Charge height mm 90.0 Quartz wall thickness mm 2.75 Hot zone temperature Kelvin 1380.0 Furnace gradient degrees/mm 0.05, 0.1, 0.5, 1.0, 2.0, 4.0 Cold furnace temperature Kelvin 1335.0 Ampoule velocity mm/sec 1.39 x 10~4, 2.78 x 10~4, 5.56 x 10~4 Table 10.20: CdTe Model Operating Parameters for examining the effect of furnace gradient and crucible velocity in the interface curvature, interface velocity and MRSS. following: 1. Establish the lowest MRSS operating conditions, in order to produce a crystal that has the lowest possible dislocation density. 2. Establish the operating conditions which generate an interface with a planar or convex shape to increase the single crystal yield. 3. Establish how the growth velocity varies during solidification. Higher velocities increase the probability that a stray crystal or other crystal defects may form during growth. The furnace gradients used ranged between 0.05 and 4.0 degree/mm and the cru-cible velocities between 1.39 x 10~4 mm/sec and 5.56 x 10 - 4 mm/sec. The hot zone temperatures used varied between 1370 K and 1390 K while the cold zone temperatures varied between 1335 K and 1360 K. The exact crystal growth conditions used are given in Table 10.20 and Table 10.21. Chapter 10. Modelling Results 118 Parameter Units Values Ampoule diameter mm 50.0 Charge height mm 90.0 Quartz wall thickness mm 2.75 Hot zone temperature Kelvin 1390.0, 1380.0, 1370, Furnace gradient degrees/mm 1.0 Cold furnace temperature Kelvin 1335.0 , 1345.0, 1355.0, 1360.0 Ampoule velocity mm/sec 2.78 x 10- 4 Table 10.21: CdTe Model Operating Parameters for examining the effect of hot and cold zone temperatures on the interface curvature and MRSS 10.4 Model Predictions Assuming No Fluid Flow The calculated shape and distribution of the isotherms for crystal growth at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 degree/mm after 2.52 x 105 seconds is shown in Figure 10.59(a). The isotherms are slightly concave upwards in the bulk of the sample indicating radial heat loss. The solid/liquid interface is slightly convex upwards. In the liquid, the isotherms become more convex. The isotherms in the solid are close together at the solid/liquid interface and separate with increasing distance from the interface. The isotherms in the liquid vary in the same manner, they are close together at the solid/liquid interface and separate with increasing distance from interface. The change in shape and position of the solid liquid interface with time from start of solidification is shown in Figure 10.59(b). The interface is convex upwards throughout growth. At the bottom and top of the crystal the convexity is the largest. The distance along the axis beween the interface positions shown for equal time increments, increases towards the end of solidification. This indicates that the axial velocity increases at the end of the crystal. The calculated MRSS contour curves corresponding to the isotherms shown in Fig-ure 10.59(a) are shown in Figure 10.59(c). The largest values of the MRSS occur in two Chapter 10. Modelling Results 119 locations, at the interface adjacent to the quartz (2.5 MPa) and lower in the solid adja-cent to the quartz (3.0 MPa). The stress is observed to drop rapidly both with distance radially from the outside to the center of the crystal, and with distance longitudinally from the lower stress concentration. 10.4.1 Effects of Crucible Velocity and Furnace Gradient The effect of changing the furnace gradient and the growth velocity on the curvature of the solid/liquid interface and the interface velocity was examined. The curvature of the solid/liquid interface is arbitrarily defined as the interface height at the side of the crystal minus the interface height at the center of the crystal. Thus a concave interface would have a positive curvature and a convex interface would have a negative curvature. The growth velocity is calculated at the center of the crystal. The effect of temperature gradient for three crucible velocities on the interface cur-vature, at the vertical midpoint of the crystal, is shown in Figure 10.60. The curvature changes from concave (positive curvature) to convex (negative curvature) at a gradient of 1 degree/mm. Large concave curvatures are observed at gradients below 0.5 degrees/mm. Above 2 degree/mm the curvature is essentially constant at approximately -3 mm. The effect of crucible velocity on the curvature is small except at very shallow gradients where the curvature increases with increased crucible velocity. The change in interface curvature with distance from the start of solidification, at different furnace gradients is shown in Figure 10.61. For a small furnace gradient of 0.05 degrees/mm, the curvature is large (10 mm) at the start of solidification decreasing slowly and progressively to near 6 mm at the end of solidification. Increasing the gradient to 0.1 degree/mm results in the same initial curvature of 10 mm at the start, but decreases more rapidly to about 0 mm at the end of solidification. At higher gradients of 1.0 to 4.0 degrees/mm the curvature is effectively flat or slightly convex for most of the crystal Chapter 10. Modelling Results 120 (a) (b) (c) Figure 10.59: Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm. (a) Isotherm contours after 2.88 x 105 sec of growth, (b) Solid/liquid interface shape and position at the lines indicated at intervals of 3.6 x 104 sec of growth, (c) MRSS contours after 2.88 x 105 sec growth, in MPa. Chapter 10. Modelling Results 121 length. The change in interface curvature during crystal growth for furnace gradients of 0.05 and 0.1 degree/mm is related to the length of the gradient zone. Figure 10.62 shows the isotherms for a gradient of 0.1 and 1 degree/mm when the interface is 40 mm from the bottom of the crucible. Figure 10.62(a) shows that the crucible is entirely in the gradient zone. Figure 10.62(b) shows the isotherms when the bottom 10 mm is in the cold zone and the top 40 mm is in the hot zone. With the crucible entirely in the gradient zone the temperature at the bottom of the crucible will always be changing thus the amount of axial heat flow will not be constant. The changing axial heat flow will result in a variable interface curvature during crystal growth as seen in Figure 10.61 for gradients less than 0.5 degree/mm. Having the bottom of the crystal in the cold zone furnace (Figure 10.62(a)) will result in a constant axial heat flux in the crystal. The constant axial heat flux will give an interface curvature that does not change during crystal growth as seen in Figure 10.61 for gradients equal to and greater than 0.5 degree/mm. The dependence of interface velocity at the center of the crystal with distance from the start of solidification is shown in Figure 10.63, for an imposed crucible velocity of 2.78 x 10 - 4 mm/sec. The calculated results show that the center velocity is relatively constant between 2.78 x 10~4 and 3.00 x 10 - 4 mm/sec for most of the crystal. Towards the end of crystal growth the interface velocity increases rapidly, 3 to 5 times the crucible velocity. At the lower furnace gradients of 0.05 and 0.1 degrees/mm the calculated interface velocity is higher than the crucible velocity, and the increase in velocity towards the end of freezing is appreciable. At the higher furnace gradients, 1 to 4 degrees/mm, the calculated velocities are the same as the crucible velocity, and the increase in velocity near the end of freezing is small. The initial variation in interface velocity at the start of solidification is due to the oscillations discussed in the sensitivity analysis. The effect of changing the furnace gradient and growth velocity on the calculated Chapter 10. Modelling Results 122 MRSS was determined, using the largest value of the MRSS (MRSS*), located at the solid/liquid interface adjacent to the quartz and the largest value of the MRSS (MRSSt) located in the body adjacent to the quartz. The results of. the MRSS* are shown in Figure 10.64. MRSS* is observed to be very small at small furnace gradients, increasing rapidly to about 2.5 MPa at a gradient of 1 degree/mm, then increasing slowly with increasing gradient to 3 MPa. The effect of crucible velocity on the MRSS* is large, the higher velocities giving 1 MPa higher values of the MRSS* and the lower velocity values give a 1 MPa lower value of the MRSS*. The MRSSt results are shown in Figure 10.65. MRSSt is observed to be very small at small furnace gradients, increasing rapidly ap-proximately linearly to about 5 MPa at a gradient of 2 degree/mm, then increasing slowly with increasing gradient to 6 MPa. The effect of crucible velocity on the MRSSt is small, the higher velocities giving slightly higher values of the MRSSt at the larger furnace gradients. The position at which the MRSSt occurs changes with the applied gradient. Fig-ure 10.66 shows the MRSS contours for furnace gradients of 4, 1 and 0.5 degrees/mm. The location of the MRSSt is always located were the gradient furnace and cold furnace join. Thus the MRSS^ is due to the changing in thermal boundary condition from the gradient furnace to the constant cold zone temperature. 10.4.2 Effects of Cold and Hot Zone Temperature The effect of changing the cold and hot zone temperature on the curvature of the solid/liquid interface , MRSS* and MRSSt was examined using the parameters listed in Table 10.21. The furnace gradient was kept at 1 degree/mm for all model runs. The effect of cold zone temperature on the interface curvature is shown in Figure 10.67. Large curvatures are observed at a high cold zone temperature, at 1360 K the interface curvature is 9 mm. Decreasing the cold zone temperature results in a decrease in the Chapter 10. Modelling Results Figure 10.60: Dependence of interface curvature on furnace gradient and crucible velo (V). Chapter 10. Modelling Results 124 Figure 10.61: Variation on interface curvature with distance from the start of solidifica-tion for different furnace temperature gradients (G). Chapter 10. Modelling Results 125 \ -(a) (b) Figure 10.62: Isotherms for crystals grown at two different gradients with furnace lo-cations at indicated positions, V = 2.78 x 10 - 4 mm/sec (a) 7.02 x 105 sec, G = 0.1 degree/mm (b) 1.98 x 105 sec, G = 1 degree/mm Chapter 10. Modelling Results 126 Figure 10.63: Vaxiation of interface velocity with distance from the start of solidification for the furnace gradients (G) shown. Imposed crucible velocity indicated by CV. Chapter 10. Modelling Results 127 0 I 1 1 1 1 i i i 0 0.5 1 1.5 2 2.5 3 3.5 4 Furnace Gradient (degree/mm) Figure 10.64: MRSS* as a function of furnace gradient for the growth velocities (V) shown. Chapter 10. Modelling Results 128 Figure 10.65: MRSSt as a function of furnace gradient for the growth velocities (V) shown. Chapter 10. Modelling Results 129 (a) (b) (c) Figure 10.66: Crystal grown at 2.78 x 10 - 4 mm/sec with various furnace gradients. The furnace location relative the the crystals are indicated in the Figure, (a) MRSS contours after 2.16 x 105 sec of growth with a 4 dergree/mm gradient, (b) MRSS contours after 2.52 x 105 sec of growth with a 1 dergree/mm gradient, (c) MRSS contours after 3.96 x 105 sec of growth with a 0.5 dergree/mm gradient. Chapter 10. Modelling Results 130 interface curvature. The curvature was decreased to approximately 0 mm at a cold zone temperature of 1335 K. The effect of changing the cold zone temperature on the calculated MRSS was deter-mined, using the values of MRSS* and MRSSt. The results are shown in Figure 10.68. The MRSSt is observed to be constant for the cold zone temperatures examined. The MRSS* is lowest at a value of 2.5 MPa using a cold zone temperature of 1335 K. The MRSS* is constant at 3 MPa for cold zone temperatures of 1345 to 1360 K. The change in MRSS* is due to changing the cold zone temperature which changes the axial to radial heat flow at the stress location. The interface is approximately flat (approximately zero radial heat flow) for cold zone temperatures of 1335 K thus lower stresses are encountered. When the interface is curved, cold zone temperature greater than 1345 K, the MRSS* is higher. The MRSSt is constant indicating the isotherms at the stress location must have the same curvature for all cold zone temperatures. The isotherms for cold zone temperatures of 1360 K, 1345 K and 1335 K are shown in Figure 10.69. Increasing the cold zone temperature increases the number of isotherms in the solid and decreases the isotherm spacing in the solid. The isotherms in the liquid remain relatively unchanged. The decrease in curvature with increasing cold zone tem-perature shown in Figure 10.67 is also apparent as the interface changes from concave to convex with decreasing cold zone temperature. The effect of the hot zone temperature on the interface curvature is shown in Fig-ure 10.70. Slightly convex curvatures are observed at a high hot zone temperature, the interface curvature being 1 mm at 1390 K. Decreasing the hot zone temperature results in changing the interface from convex to concave. The curvature decreased to approximately -6 mm at 1370 K. The effect of changing the hot zone temperature on the calculated MRSS was deter-mined based on the MRSS* and the MRSSt. The results are shown in Figure 10.71. The Chapter 10. Modelling Results 131 1335 1340 1345 1350 1355 1360 Cold Zone Temperature (K) Figure 10.67: Dependence of interface curvature as a function of cold furnace tempera-ture. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V = 2.78xl0_4mm/sec Chapter 10. Modelling Resxdts 132 3.1 3 0 Cu 00 00 2.9 2.8 2.7 2.6 2.5 2.4 1335 1340 1345 1350 1355 Cold Zone Temperature (K) 136 Figure 10.68: MRSS* and MRSSt as a function of cold furnace temperature. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V = 2.78xl0~4mm/sec Chapter 10. Modelling Results 133 Figure 10.69: Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm, hot zone temperature of 1380 K, (a) cold zone temperature of 1360 K, (b) cold zone temperature of 1345 K, (c) cold zone temperature of 1335 K. Chapter 10. Modelling Results 134 MRSS* is observed to be highest at a lowest hot zone temperature. At a hot zone tem-perature of 1370 K MRSS* is 3 MPa, decreasing to 2.5 MPa at a hot zone temperature of 1390 K. The MRSSt behaves the opposite of the MRSS*. At a low hot zone temperature of 1370 K the MRSS^ is 2.5 MPa and increases to 3 MPa at a hot zone temperature of 1390 K. The change in MRSS* is due to changing the hot zone temperature which changes the axial to radial heat flow at the stress locations. The interface is approximately flat (approximately zero radial heat flow) for hot zone temperatures of 1380 K and 1390 K thus lower stresses are encountered. When the interface is highly curved, hot zone temperature of 1370 K, the MRSS* is higher. The change in the level of the MRSS^ is also due to the change in the axial to radial heat flow. At a high hot zone temperature the amount of radial to axial heat flow is high at the location of the MRSSt thus increasing the curvature of the isotherms. This results in a higher MRSSt level. At lower hot zone temperatures the radial heat flow is less than the axial thus resulting in flatter isotherms at the MRSSt. The flatter isotherms result in a lower MRSS. The isotherms for hot zone temperatures of 1390 K, 1380 K and 1370 K are shown in Figure 10.72. Decreasing the hot zone temperature has little effect on the temperature isotherms in the solid far from the solid liquid interface. The temperature isotherms in the liquid and solid close to the melting point get farther apart. The change from a concave interface to a convex interface with decreasing hot zone temperature as shown in Figure 10.70 can also been seen. 10.4.3 Excess Shear Stress Results The MRSS and the excess stress above the critical shear stress is compared for a Cd!_ a : Znx Te crystal grown with a thermal gradient of 1 degrees/mm, a hot zone temperature of 1380K, a cold zone temperature of 1335K, and a growth rate of 2.78 x 10~4 mm/sec Chapter 10. Modelling Results 135 Figure 10.70: Dependence of interface curvature as a function of hot furnace temperature. Cold zone temperature is 1335 K, G = 1.0 degrees/mm, V = 2.78xl0 - 4 mm/sec Chapter 10. Modelling Results Figure 10.71: MRSS* and MRSSt as a function of hot furnace temperature. Cold temperature is 1335 K, G •== 1.0 degrees/mm, V = 2.78xl0 - 4 mm/sec Chapter 10. Modelling Results 137 Figure 10.72: Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm, cold zone temperature of 1335 K, (a) hot zone temperature of 1390 K, (b) hot zone temperature of 1380 K, (c) hot zone temperature of 1370 K. Chapter 10. Modelling Results 138 at 2.52 X 105 seconds after solidification started. The results are shown in Figure 10.73, in which the MRSS contours are shown in (a) and the excess stress contours in (b) using the CRSS value as experimentally measur ed (o-crss = 0.0232exp (^2)MPa) . In (b), the excess stress is significant in the regions close to the outer wall, as with the MRSS. The MRSS is relatively low thus subtracting the CRSS produces a decrease of 0.5 MPa in the contour lines. The region close to the wall is where dislocations will be generated. In the lower part of the crystal there is lower excess stress, therefore dislocations would be generated at a lower rate in this region during growth. This would indicate that there should be a higher dislocation density in the crystal near the side, compared to the centre, providing dislocation annihilation does not occur. 10.5 Model Predictions using a Bulk Heat Transfer Coefficient Seven Times the Static Value Calculations were done using a liquid conductivity seven times the static value to take into account bulk heat transport due to fluid flow. The calculated shape and distribution of the isotherms for crystal growth at 2.78 x 10~4 mm/sec with a furnace gradient of 1 degree/mm after 2.88 x 105 seconds is shown in Figure 10.74(a). The isotherms are concave upwards in the bulk of the sample including the solid/liquid interface, indicating that radial heat loss occurs. Near the top of the sample, in the melt, the isotherms become convex upwards and are widely spaced. The vertical central gradient is relatively steep in the solid below the interface, decreasing with distance in the solid away from the interface, and shallow in the melt. The discontinuity in the isotherm at the solid CdTe/quartz interface is due to the presence of an air gap between the solid and the quartz. Chapter 10. Modelling Results 139 (a) (b) Figure 10.73: (a) MRSS contours in MPa (b) Excess stress above the CRSS (MPa). G = 0.5 degree/mm, V = 2.78 x 10~4 mm/sec, Hot Zone = 1370 K, Cold Zone = 1360 K, 3.6 x 105 sec from start of crystal growth. (acrss = 0.0232 exp MPa) Chapter 10. Modelling Results 140 The shape and position of the solid/liquid interface with time from the start of so-lidification is shown in Figure 10.74(b). The interface is observed to be concave upwards throughout growth tending to become slightly flatter towards the end of solidification. The distance along the axis between the interface positions shown for equal time incre-ments, increases towards the end of solidification. This indicates that the axial velocity increases in this region. The calculated MRSS contour curves corresponding to the isotherms shown in Fig-ure 10.74(a) are shown in Figure 10.74(c). The largest values of the MRSS (18 MPa) are observed to occur adjacent to the quartz near the solid/liquid interface. The stress is observed to drop rapidly both with distance radially from the outside to the centre of the crystal, and with distance longitudinally into the crystal from the interface. 10.5.1 Effects of Crucible Velocity and Foirnace Gradient The effect of temperature gradient for three crucible velocities on the interface curvature, at the vertical midpoint of the crystal, is shown in Figure 10.75. Large curvatures are observed at gradients below 0.5 degrees/mm. Above 0.5 degree/mm the curvature is essentially constant at near 8 mm. The effect of crucible velocity on the curvature is small except at very shallow gradients where the curvature increases with increased crucible velocity The change in interface curvature with distance from the start of solidification, at different furnace gradients is shown in Figure 10.76. For a small furnace gradient of 0.05 degrees/mm, the curvature is large (18mm) at the start of solidification decreasing slowly and progressively to near 12 mm at the end of solidification. Increasing the gradient to 0.1 degree/mm results in the same initial curvature of 18 mm at the start, but decreases more rapidly to near 8 mm at the end of solidification. At higher gradients of 0.5 to 4 degrees/mm the curvature is effectively constant at 8 mm for most of the crystal length, Chapter 10. Modelling Results 141 to (b) (c) Figure 10.74: Crystal grown at 2.78 x 10~4 mm/sec with a furnace gradient of 1 de-gree/mm using 7 x the normal liquid conductivity, (a) Isotherm contours after 2.88 x 105 sec of growth, (b) Solid/liquid interface shape and position at the lines indicated at in-tervals of 3.6 x 104 sec of growth, (c) MRSS contours after 2.88 x 105 sec growth, in MPa. Chapter 10. Modelling Results 142 dropping to 2 mm near the end of solidification. The dependence of interface velocity at the center of the crystal with distance from the start of solidification is shown in Figure 10.77, for an imposed crucible velocity of 2.78 x 10 - 4 mm/sec. The calculated results show that the center velocity is relatively constant between 2.78 x 10 - 4 and 4.17 x l O - 4 mm/sec for most of the crystal. Towards the end of crystal growth the interface velocity increases rapidly, 3 to 5 times the crucible velocity. At the lower furnace gradients of 0.05 and 0.1 degrees/mm the calculated interface velocity is higher than the crucible velocity, and the increase in velocity towards the end of freezing is appreciable. At the higher furnace gradients, 1 to 4 degrees/mm, the calculated velocities are the same as the crucible velocity, and the increase in velocity near the end of freezing is small. The effect of changing the furnace gradient and growth velocity on the calculated MRSS was determined, using the largest value of the MRSS (MRSS*), located at the solid/liquid interface adjacent to the quartz. The results are shown in Figure 10.78. MRSS* is observed to be very small at small furnace gradients, increasing rapidly and approximately linear to about 18 MPa at a gradient of 1 degree/mm, then increasing slowly with increasing gradient to 24 MPa. The effect of crucible velocity on the MRSS* is small, the higher velocities giving slightly higher values of the MRSS* at the larger furnace gradients. 10.5.2 Effects of Cold and Hot Zone Temperature The effect of changing the cold and hot zone temperature on the curvature of the solid/liquid interface and MRSS* was examined using the parameters listed in Table 10.21. The furnace gradient was kept at 1 degree/mm for all model runs. The effect of cold zone temperature on the interface curvature is shown in Figure 10.79. Large curvatures are observed at a high cold zone temperature, at 1360 K the interface curvature is 22mm. Chapter 10. Modelling Results 143 18 s G ID t-i 3 +-» c3 > t-l u o 16 14 12 10 • V = 5.56 x 10-4 mm/sec • V = 2.78 x 10~4 mm/sec • V = 1.39 x 10-4 mm/sec -o- -9 0.5 1 1.5 2 2.5 3 Furnace Gradient (degree/mm) 3.5 Figure 10.75: Dependence of interface curvature on furnace gradient and crucible velocity (V) using 7 x the normal liquid conductivity. Chapter 10. Modelling Results 144 Figure 10.76: Variation on interface curvature with distance from the start of solid-ification for different furnace temperature gradients (G) using 7 x the normal liquid conductivity. Chapter 10. Modelling Results 145 Figure 10.77: Variation of interface velocity with distance from the start of solidification for the furnace gradients (G) shown 7 x the normal liquid conductivity used. Imposed crucible velocity indicated by CV. Chapter 10. Modelling Results 146 Figure 10.78: MRSS* as a function of furnace gradient for the growth velocities (V) shown. 7 x the normal liquid conductivity. Chapter 10. Modelling Results 147 Decreasing the cold zone temperature results in a decrease in the interface curvature. The curvature was decreased to approximately 9 mm at a cold zone temperature of 1335 K. The effect of changing the cold zone temperature on the calculated MRSS was deter-mined, using the value of MRSS*, located at the solid/liquid interface adjacent to the quartz. The results are shown in Figure 10.80. The MRSS* is observed to be smallest at a high cold zone temperature. At a cold zone temperature of 1360 K the MRSS* is 4 MPa, the MRSS* increasing approximately linear to about 16 MPa at a cold zone temperature of 1335 K. The isotherms for cold zone temperatures of 1360 K, 1345 K and 1335 K are shown in Figure 10.81. Increasing the cold zone temperature increases the number of isotherms in the solid and decreases the isotherm spacing in the solid. The isotherms in the liq-uid remain relatively unchanged. The decrease in curvature with increasing cold zone temperature shown in Figure 10.79 is also apparent. The effect of the hot zone temperature on the interface curvature is shown in Fig-ure 10.82. Large curvatures are observed at a high hot zone temperature, the interface curvature being 11mm at 1390 K. Decreasing the hot zone temperature results in a decrease in the interface curvature, to approximately 5 mm at 1370 K. The effect of changing the hot zone temperature on the calculated MRSS was deter-mined based on the MRSS* located at the solid/liquid interface adjacent to the quartz. The results are shown in Figure 10.83. The MRSS* is observed to be smallest at a lowest hot zone temperature. At a hot zone temperature of 1370 K MRSS* is 12 MPa, increasing linearly to about 20 MPa at a hot zone temperature of 1390 K. The isotherms for hot zone temperatures of 1390 K, 1380 K and 1370 K are shown in Figure 10.84. Decreasing the hot zone temperature has little effect on the temperature isotherms in the solid far from the solid liquid interface. The temperature isotherms in Chapter 10. Modelling Results 148 1335 1340 1345 1350 1355 1360 Cold Zone Temperature (K) Figure 10.79: Dependence of interface curvature as a function of cold furnace temperature using 7 x the normal liquid conductivity. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V : 2.78xl0_4mm/sec Chapter 10. Modelling Results 149 Figure 10.80: MRSS* as a function of cold furnace temperature using 7 x the normal liq-uid conductivity. Hot zone temperature is 1380 K, G = 1.0 degrees/mm, V = 2.78 x l O - 4 mm/sec Chapter 10. Modelling Results 150 (a) (b) (c) Figure 10.81: Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm using 7 x the normal liquid conductivity, hot zone temperature of 1380 K,.(a) cold zone temperature of 1360 K, (b) cold zone temperature of 1345 K, (c) cold zone temperature of 1335 K. Chapter 10. Modelling Results 151 1370 1372 1374 1376 1378 1380 1382 1384 1386 1388 1390 Hot Zone Temperature (K) Figure 10.82: Dependence of interface curvature as a function of hot furnace temperature using 7 x the normal liquid conductivity. Cold zone temperature is 1335 K, G = 1.0 degrees/mm, V = 2.78 x l O - 4 mm/sec The isotherms for hot zone temperatures of 1390 K, 1380 K and 1370 K are shown in Figure 10.84. Decreasing the hot zone temperature has little effect on the temperature isotherms in the solid far from the solid liquid interface. The temperature isotherms in the liquid and solid close to the melting point get farther apart. The decrease in curvature with decreasing hot zone temperature as shown in Figure 10.82 can also been seen. A hot zone temperature of 1370 K results in an almost flat interface. Chapter 10. Modelling Results 152 1370 1372 1374 1376 1378 1380 1382 1384 1386 1388 1390 Hot Zone Temperature (K) Figure 10.83: MRSS* as a function of hot furnace temperature using 7 x the normal liquid conductivity. Cold zone temperature is 1335 K, G = 1.0 degrees/mm, V = 2.78 xlO~4 mm/sec Chapter 10. Modelling Resvdts 153 (a) (b) (c) Figure 10.84: Crystal grown at 2.78 x 10 - 4 mm/sec with a furnace gradient of 1 de-gree/mm using 7 x the normal liquid conductivity, cold zone temperature of 1335 K, (a) hot zone temperature of 1390 K, (b) hot zone temperature of 1380 K, (c) hot zone temperature of 1370 K. Chapter 10. Modelling Results 154 of 1380K, a cold zone temperature of 1335K, and a growth rate of 2.78 x 10~4 mm/sec at 3.6 x 105 seconds after solidification started. The results are shown in Figure 10.85, in which the MRSS contours are shown in (a) and the excess stress contours in (b) using the CRSS values as experimentally measured (cxcrss- = 0.0232exp ^ 4 ™ , 0 JMPa) . In (b), the excess stress is significant in the region of the interface close to the outer wall, as with the MRSS. Since the CRSS is small relative to the MRSS, subtracting its value from the MRSS produces a relatively small shift in the stress contours. This is the region where dislocations will be generated. In the lower part of the crystal there is no excess stress, therefore dislocations would not be generated in this region during growth. This would indicate that there should be a higher dislocation density in the crystal near the side, compared to the centre, providing dislocation annihilation does not occur. The high dislocation density associated with CdTe indicates there is little recovery during and after crystal growth. Chapter 10. Modelling Results 155 (a) (b) Figure 10.85: (a) MRSS contours in MPa (b) Excess stress above the CRSS (MPa) using 7 x the normal liquid conductivity. G = 0.5 degree/mm, V = 2.78 x 10~4 mm/sec, Hot Zone = 1370 K, Cold Zone = 1360 K, 3.6 x 105 sec from start of crystal growth. (acrss = 0.0232 exp ( ^ M ) MPa) Chapter 11 Discussion 11.1 Discussion of Experimental Results The critical resolved shear stress of CdTe is low, having a value of 0.43 MPa at the melting point of CdTe (1365 K). The low value indicates dislocations will be generated at low thermal stresses during crystal growth. The etch pit density as measured is approximately 105 cm - 2 . The CRSS of GaAs at its melting point is approximatley double that of CdTe (0.9 MPa) and in turn the average EPD in CdTe is an order of magnitude higher than in GaAs which is about 104 cm - 2 . The CdTe etch pit results indicate a lower dislocation density at the center of the crystal and a higher density at the outside. In GaAs a four fold S3anmetry is found for a (100) surface when the <100> direction is parallel to the crystal growth direction [4]. The degree of symmetry changes with the orientation of the single crystal. Symmetry in the etch pattern was not observed in CdTe due to the orientation of the slice relative to the crystal and the low number of samples examined. The nonuniformity of the etch pit density shown in Figure 7.18 (Slice 521sll) and Figure 7.21 (Slice 555sll) may be due some type of symmetry. The nonuniformity due to variations in the etch pit size must be controlled in any future etch pit mapping. A variation in etch pit size results in a false change in etch pit density. A pit having an area larger than average would be detected as more than one pit. 156 Chapter 11. Discussion 157 11.2 Discussion of Model Predictions Assuming no Fluid Flow The model calculations for low axial temperature gradients show that the solid/liquid interface is relatively flat in the centre of the crystal. This indicates that there is no significant radial heat loss at the interface. Increasing the furnace gradient above 1 degree/mm changes the interface to a convex shape indicating that heat flows towards the interface in the radial direction. At gradients less than 1 degree/mm the interface moves to a concave shape indicating that heat flows in the radial direction away from the interface. In the upper portion of the crystal, toward the end of solidification, the interface cur-vature tends to become more convex, the growth velocity increases, and the temperature gradient becomes smaller. These effects can be attributed to the smaller amount of heat entering into the upper part of the boule containing CdTe liquid. There are two locations in the solid were there is a maximum in the MRSS values. One location is at the interface adjacent to the quartz. The other maximum values of the MRSS occur at the outside of the crystal that is adjacent to the location where the gradient furnace and the cold furnace join. The MRSS increases with furnace gradient for gradients up to 4.0 degrees/mm. The interface curvature decreases with increasing furnace gradient in the same range. The controlling factor on the MRSS is the furnace gradient, and not the interface curvature. The shape of the interface can be manipulated by adjusting the hot and cold zone furnace temp ear atures. The hot zone is the heat source for Vertical Bridgman crystal growth while the cold zone acts as a heat sink. By decreasing the hot zone temperature the amount of heat that must be dissipated is lowered, resulting in a more convex interface for a given cold zone temperature. Similarly, the cold zone can be adjusted to control the amount of axial heat flow. The cold zone temperature is more important than the Chapter 11. Discussion 158 hot zone temperature in determining the MRSS since the cold zone temperature affects the temperature isotherms in the solid while the hot zone temperature affects the liquid temperature isotherms. The CRSS in CdTe is very low and therefore dislocations will be generated with small thermal stresses in the crystal. With no fluid flow in the melt, the results show that the levels of the MRSS are low and that subtracting the CRSS lowers the stress by 1 MPa. The optimum operating conditions for Vertical Bridgman growth of CdTe are a low cold zone furnace temperature, a low hot zone furnace temperature and a moderate gradient (approximatly 1 degree/mm) in the gradient furnace. 11.3 Discussion of Model Predictions using a Bulk Heat Transfer Coefficient Seven Times the Value with no Fluid Flow The concave upward curvature of the solid/liquid interface indicates that the radiative heat loss through the quartz near the interface is comparable to the axial heat conduc-tion through the solid. Decreasing the furnace gradient Avou ld increase the radial heat radiation and thus increase the interface curvature, in agreement with the model results. In the upper portion of the crystal, toward the end of solidification, the interface curvature tends to flatten, the growth velocity increases, and the temperature gradient becomes smaller. The effects can be attributed to the lower amount of heat entering into the upper part of the boule containing CdTe liquid. The MRSS is maximum at the interface and increases with furnace gradient for gra-dients u p to 1.5 degrees/mm. The interface curvature decreases with increasing furnace gradient in the same range. Therefore the controlling factor of the MRSS is the furnace gradient, and not the interface curvature. The hot zone is the heat source for Vertical Bridgman crystal growth while the cold Chapter 11. Discussion 159 zone acts as a heat sink. By decreasing the hot zone temperature the amount of heat that must be dissipated is lowered, resulting in a flatter interface for a given cold zone temperature. Similarly the cold zone can be adjusted to control the amount of axial heat flow. The cold zone temperature is more important than the hot zone temperature for the MRSS since the cold zone temperature affects the temperature isotherms in the solid while the hot zone temperature affects the liquid temperature isotherms. The CRSS in CdTe is very low and therefore dislocations will be generated with small thermal stresses in the crystal. It is therefore evident that gradients should be kept as small as possible. However, with low furnace gradients, high interface curvatures are obtained which could lead to the enhanced production of crystal defects such as twins, stray crystal formation, non-stoichiometry and porosity. 11.4 Comparison Between the Two Sets of Model Results The increased heat flow due to fluid flow increases the MRSS by approximatly 7 times. Fluid flow in the liquid makes the attainment of a flat or convex interface impossible. The MRSSt that occurs under no fluid flow condition does not occur when the conductivity of the liquid is increased 7 times to approximate the bulk heat transport from fluid flow. The interface is lowered due to the increased heat flow toward the solid/liquid interface, thus moving the MRSS* to a position adjacent to the location were the gradient furnace and cold furnace join. This location is the normal position of the MRSSt. In both sets of calculations the furnace gradient has a greater effect on the axial to radial heat flow ratio than the hot and cold zone temperatures. The spacing of the isotherms and MRSS levels are dramatically changed by changing the furnace gradient. The hot and cold zone temperatures change the the isotherm spacing and MRSS level by a negligable amount. Chapter 12 Summary and Conclusions The results of the present investigation can be summarized as follows: 1. A mathematical model has been developed to determine the thermal and stress fields in CdTe during Vertical Bridgman crystal growth. 2. The critical resolved shear stress has been measured in CdTe in the temperature range of 738K to 1024K. The extrapolated value of the CRSS at the melting point of CdTe is 0.43 MPa, which is low. 3. The etch pit density was found to be high at the outside and lower at the center of the samples. 4. Stresses in multigrain samples are very high close to the grain boundaries as shown by higher EPD in the grain boundary region. 5. Results Assuming No Fluid Flow Occurs (a) The interface is shown to be flat or convex to the melt and remains so during the entire growth process under the regime we have studied. The interface curvature increases with decreasing furnace gradient at low gradients. The growth velocity is constant for most of the solidification, increasing signifi-cantly towards the end of freezing. (b) The maximum shear stress in CdTe occurs at two locations. One in near the solid/liquid interface close to the outer surface of the crystal. The second is at 160 Chapter 12. Summary and Conclusions 161 the outer surface of the crystal adjacent to the portion of the furnace where the gradient furnace and cold furnace joins. The maximum stress is greater than the CRSS in this region. (c) The furnace gradient should be kept low to keep the MRSS levels low and thus keep dislocation formation during crystal growth at a minimum. (d) Changing the crucible velocity changes the MRSS* level by approximately 1 MPa when the velocity is doubled. The MRSSt level changes by 0.5 MPa when the velocity is doubled. Since the level of the MRSS is relativley low (4 MPa) these variations are significant. (e) The hot and cold zone temperature can be adjusted to reduce the interface curvature. Increasing the hot zone temperature makes the interface concave. Decreasing the cold zone temperature makes the interface convex. 6. Results for Model Predictions using a Bulk Heat Transfer Coefficient Seven Times the Value for No Fluid Flow. (a) The interface is shown to be concave to the melt and remains so during the en-tire growth process under the regime we have studied. The interface curvature increases with decreasing furnace gradient at low gradients. The growth ve-locity is constant for most of the solidification, increasing significantly towards the end of freezing. (b) The maximum shear stress in CdTe occurs near the solid/liquid interface close to the outer surface of the crystal. The maximum stress is greater than the CRSS in this region. (c) The furnace gradient must be small to keep the MRSS levels low and thus avoid excessive dislocation formation. Chapter 12. Summary and Conclusions 162 (d) Changing the crucible velocity changes the level of the MRSS* by approximaley 1.5 MPa when the velocity is doubled. High levels of the MRSS (20 MPa) make this change insignificant. (e) The hot and cold zone temperature can be adjusted to reduce the interface curvature. Increasing the hot zone temperature makes the interface more con-cave. Decreasing the cold zone temperature makes the interface less concave. 7. With fluid flow in the melt, the MRSS is estimated to be increased by a factor of seven, as compared to the case of no fluid flow, based on reported data on tin. 8. The results of the present investigation indicate that optimum growth conditions are obtained when the hot zone temperature is low, the cold temperature is low, the temperature gradient in the gradient zone is at 1 degree/mm and no fluid flow occurs. This will give a low dislocation density with a flat to convex interface. Chapter 13 Recommendations for Future Work Several more sets of experiments could be carried out. Most of the experiments would be directed toward proving the mathematical model. 1. The interface should be profiled at different times during crystal growth. Ther-mocouples should be used to record the temperatures at various radial locations during crystal growth. This would allow a model to accurately predict Vertical Bridgman Growth. 2. The etch pit density and distribution should be fully mapped in a CdTe single crystal boule. 3. A fluid flow model could be developed to establish the extent of flow in the melt during solidification. This could then be related to the heat flow through the melt in the model replacing the constant factor used in the present investigation. 163 Bibliography [1] K. Zanio, Semiconductors and Semimetals, Vol. IS, Cadmium Telluride, Academic Press (1978). [2] 0. Oda, K. Hirata, K. Matsumoto and I. Tsuboya, J. Crystal Growth, 71, 273 (1985). [3] N.R. Kyle, J. Electrochem. Soc, 118, 1790 (1971). [4] C. Schvezov, PhD Thesis, The University of British Columbia, (1986). [5] C.E. Huang, D. Elwell and R.S. Feigelson, Journal of Crystal Growth, 64, 441 (1983). [6] M.J. Crochet, F.T. Geyling and J.J. Van Schaftingen, Finite Elements in Fluids, 6, 321 (1985). [7] T. Jasinski and A.F. Witt, Journal of Crystal Growth , 71, 295 (1985). [8] Y.M. Dakhoul, R. Farmer, S.L. Lehoczky and F.R. Szofran, Journal of Crystal Growth , 86, 49 (1988). [9] S. Sen, W.H. Konkel, S.J. Tighe, L.G. Bland, S.R. Harma and R.E. Taylor, Journal of Crystal Growth , 86, 111 (1988). [10] S. Motakef, Journal of Crystal Growth , 98, 711 (1989). [11] M.J. Crochet, F. Dupret and Y. Ryckmans, Journal of Crystal Growth , 97, 173 (1989). [12] C.E. Huang, D. Elwell and R.S. Feigelson, Journal of Crystal Growth, 69, 275 (1984). [13] C. Schvezov, I.V. Samarasekera and F. Weinberg, Journal of Crystal Growth , 84, 219 (1987). [14] A. Tanaka, Y. Masa, S. Seto and T. Kawasaki, Mat. Res. Soc. Symp. Proc, 90, 111 (1987). [15] E.Y. Gutmanas, N. Travitzky, U. Plitt and P. Haasen, Scripta Metallurgica, 13, 293 (1979). [16] E.L. Hall and J.B. Vander Sande, Journal of The American Ceramic Society, 61, 418 (1978). 164 Bibliography 165 [17] D.J. Williams and A.W. Vere, Journal of Crystal Growth , 83, 341 (1987). [18] L. Li, Personal communication to Chris Parfeniuk, June 13, 1990. [19] M. Inoue, I. Teramoto, and S. Takayanagi, Journal of Applied Physics, 33, 2578 (1962). [20] Zienkiewicz, The Finite Element Method, McGraw-Hill Book Company, New York, (1977). [21] Zienkiewicz and Morgan, Finite Elements and Approximation, John Wiley & Sons, New York, (1983). [22] J.M. Bettencourt, O.C. Zienkiewicz and G. Cantin, Int. J. Num. Meth. Engng., 17, 931, (1981) [23] E.C. Lemmon Num. Methods in Heat Transfer, , 201-213 (1981). [24] R.W. Lewis and P.M. Roberts Applied Scientific Research, 44, 61-92 (1987). [25] T.W. Clyne Materials Science and Engineering,65, 111-124 (1984). [26] Huebner and Thornton, The Finite Element Method For Engineers, John Wiley &; Sons, New York, (1982). [27] F. Kreith and Z. Black, Basic Heat Transfer, Harper and Row, New York, (1980). [28] B.G. Thomas, I.V. Samarasekera and J.K. Brimacombe, Met. Trans. B, 15b, 307, (1984). [29] H.S. Carslaw and J.C. Jeager, Conduction of Heat in Solids, Oxford University Press, Glasgow, (1959). [30] Y. Nishida, W. Droste, and S. Engler, Met. Trans. B, 17B, 833 (1986). [31] E.M. Sparrow, L.U. Albers and E.R.G. Eckert, Journal of Heat Transfer, 84, (1962), 73-81 [32] F. Krieth, W. Z. Black, Basic Heat Transfer, Harper & Row, New York, (1980). [33] W.A. Brantley, J. Appl. Phys., 44, 534 (1973). [34] H. McSkimin, D. Thomas, J. Appl. Phys., 33, 56 (1961). [35] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, (1970). Bibliography 166 [36] D.D. Wagrnan et al., Selected Values of Thermodynamic Properties, National Bureau of Standards Series 270, U.S. Department of Commerce, Washington, (1968-1971). [37] Y.S. Touloukian, Series Editor, Thermophysical Properties of Matter, the TPRC data series, IFI/Plenum, New York, (1970-). [38] C. Harrison and F. Weinberg, Met. Trans. B, 16b, 355, (1985). [39] R.C. Weast, ed., CRC Handbook of Chemistry and Physics, 63r<^ Edition, CRC Press, Inc, Florida, (1982-83). [40] B.D. Cullity, Elements of X-Ray Diffraction, Addison-Wesley, Massachusetts, (1978). [41] T.W. Chapman, Am. Inst. Chem. Eng. J., 12(2), 395 (1966). Appendix A Formulation of Thermal F E M Model A . l Galerkins Method Heat transfer in a three-dimensional isotropic material Cl bounded by a surface T as shown in Figure A.86, is described by the equation: dT kV*T = P C p ^ (A.8) where k = conductivity, p = density, Cp = the specific heat, and T = temperature. Transforming the heat transfer equation to to cylindrical coordinates yields the following: k d ( dT\ k d2T 7 d2T „ dT rfr [rfr) + r W + ^ = ? C p ~d~t ( A " 9 ) where r, and <fr are the coordinates as shown in Figure A.86. The heat equation in cylindrical coordinates is very useful if the domain to be modelled is symmetric about the z axis. The axi-symmetric domain combined with axi-symmetric boundary conditions allows the assumption that the temperature in the (f> direction is constant, thus ^ § ^ 5 - = 0. Equation A.9 becomes k d ( er\ ,d'T „ dT rTr r * ) + h ^ = p C ' T t ( A 1 0 ) or k d ( dT\ d2T „ dT n . . \ r d - r \ : ^ ) + k l ^ - p C > l H = 0 ( A - U ) Solving equation A.11 using Galerkins method requires that the dependent variable, temperature (T), be approximated by a trial (shape) function (Nj). Thus at all points 167 Appendix A. Formulation of Thermal FEM Model Figure A.86: Domain for a three-dimensional field problem Appendix A. Formulation of Thermal FEM Model 169 within the domain fl the temperature is approximated by r « E ^ i ( ^ ) (A-i2) J=I Differentiating the approximation with respect to r, z, and t and substituting them into equation A.11 will not give an exact solution. The error, or residual R, that is introduced by using the approximation (trial function) is given by - g ^ ^ ^ . ^ f f ) ) - , ^ ^ } ^ ( A , 3 , The value of the residual (R) is a function of its position in the domain (fl). To reduce the residual to approximately zero over the whole domain the weighted residual method is used. The weighted residual method integrates the residual a multiple number of times. Each integral is weighted (Wi) such that the residual is reduced to zero as shown in equation A. 14. / WxRMl + I W2Rdn + . . . « 0 Jn Jn m . / WiRdn^O (A. 14) The weighting functions (Wi) are a set of independent functions that gives convergence as the number used increases (m —• oo). Galerkins method is a specific case of the weighted residual method, (equation A.13) in which the trial function (Nj) is used as the weighting function (Wi). [ N,RdD, « 0 (A.15) .=i J{> Substituting the full form of the residual into equation A.15 gives: n n -••=i;=i J q where i = 1,2,..., n and j = 1,2,..., n Note that summation signs are not displayed in any further equations for simplicity. Also, Appendix A. Formulation of Thermal FEM Model 170 dCl will be put in proper form for cylindrical coordinates dCl — dVol = dzr dr d<f). For axi-symmetry about the z axis dz r dr d<f) = 2 TT r dz dr - u-h (4 w (t)) }** -The r term from the integration with respect to 6 cancels out the £ term on the r differential simplifying the subsequent manipulation of equation A. 16. Integration by parts in two dimensions reduces equation A. 16 to a first order derivative problem. In general integration by parts is given by the Equation A.17. U/i*>—U/i^L Isa^  (A-17) r is the boundary surrounding the area described by the coordinates £ and rj. is the r] component of the normal to the surface. Equation A.16 after integration by parts takes on the form r dN- r dN -2TT / rNikTj—-*-nrdT - 2TT rNikT.-^n.dT sa 0 (A.18) Jr dr Jr dz -The surface integrals from equation A.18 represent the natural boundary conditions which are, specified heat flow at the surface, or convective heat transfer; only convective heat transfer will be examined. Applying Fouriers law to terms in the surface integrals in Equation A.18 , „dNj , m.dJV.- I m 3 / V j dr dz dn gives the convective heat transfer term. dN-kT3—^ = -hc (TjNj-T^) on where n is the normal direction to the surface and hc is the convective heat transfer coef-ficient. Radiative heat transfer is incorporated into the formulation by using a radiation Appendix A. Formulation of Thermal FEM Model 171 heat transfer coefficient. Radiation heat transfer is given by: FaaT4 = -Fcto ( T 4 - T ^ ) where a - Stefan-Boltzman constant = 5.67 x 10~8W / m 2 K 4 a = absorptivity of surface F = radiation view factor to account for the relative geometries of the emitting and absorbing surfaces The radiative heat transfer equation can be rewritten as Thus the convective heat formulation can be used for radiation. The final form of the heat transfer equation that has been solved is: Faa (T 4 - 7 £ ) = hr(T — T^) where, hr, the radiation heat transfer equation is given by: (A.19) 2TT uAiTs dN, dNj dNi dNj dr dr dz dz + PCpNiNj dTj dt dr dz (A.20) For simplicity equation A.20 is transformed into matrix form. where [K] {T} + [C] {T} = {f} [K] = [IQ + [Kh] (A.21) Appendix A. Formulation of Thermal FEM Model 172 {/} = CM + {fr} The elements of the matrices and vectors are drdz hi = I rhc NiT^dT or / rhr NiT^dT Jr Jr The conduction matrix [Kc], convection matrix [Kh], convection force vector {//,}, and capacitance matrix [C] will be formulated. A . 2 Finite Element Method The Finite Element Method is the Weighted Residual Method applied to a descrete number of elements. Each element is formulated simultaneously and since the elements are interconnected therefore giving the solution for the entire domain. Using the Finite Element Method allows an increase in accuracy while still using a small number of unique shape functions (N). The actual formulation of the equations are on the elemental level and not the whole domain. The representation of the finite element equations on the elemental level are represented by a superscript (e): A.3 Shape (Trial) Function The type of shape function used must satisfy two criteria: continuity and completeness. Continuity refers to how many times the shape function can be differentiated and still [Kp {T} + [C] ( e ) {T} = {/}(e) (A.22) Appendix A. Formulation of Thermal FEM Model 173 be continuous. To solve a differential equation that has a highest order of derivative equal to m a shape function continuous up to the m — 1 derivative is required. The notation used is C m _ 1 continuity. Completeness requires that the polynomial used as the shape function be a complete polynomial to the order m. Equation A.20 is a first order derivative (m=l). Thus the shape function used must have at minimum C° continuity and be a complete polynomial up to the order 1. Therefore a simple linear shape function could be used to solve the differential equation (ie T = a + bx + cy). The shape function used in this study is a quadratic polynomial. Using higher order polynomials as shape functions will give a more accurate solution. A quadratic polyno-mial is required for the thermal stress analysis (Chapter B) since it requires C1 continuity. The actual shape function is a function of £ and TJ and isoparametric in nature. There are eight shape functions (equation A.23) that give a four sided element (Figure A.87a). The sides of the element are curved due to the quadratic nature of the shape functions. The approximation (shape function) for temperature is: T = J2T,Ni(t,r]) t=i The corresponding values of N are: Wi = -(1 - 0(1 - + ^ + ^)/4 N2 = -(1 + 0(1 - »?)(! - e + »7)/4 N3 = -(1 + 0(1+ i/)(l-e-i/)/4 N4 = _(l_0(l + i7)(l + £-i7)/4 N5 = N6 = (WXi + 0/2 N7 = (l-n(l + >?)/2 N8 = (i-VXi-0/2 (A.23) To treat radiative/convective boundary conditions one dimensional shape functions Appendix A. Formulation of Thermal FEM Model 174 Figure A.87: Element Shape a) Quadratic isoparametric element b)Quadratic isopari-metric boundary element Appendix A. Formulation of Thermal FEM Model 175 have been employed. The one dimension element is shown in Figure A.87b and the associated shape functions are given in equation A.24. Using quadratic isoparametric shape functions complicates the equation solution technique. The domain is in terms of r and z and the element shape functions are in terms of £ and rj. This requires that the element be mapped into £ and n space for integration. The integration that is required must be done numerically. N1 = -V(l-ri)/2 N2 = T?(1 + v)/2 N3 = l-V 2 A .4 Numerical Integration and Mapping Functions Since boundaries of the element are curved the nodes in the r-z plane must be mapped to the £-77 plane were they would be isoparametric in nature (Figure A.88). The mapping function must be quadratic since the curved boundaries of the element in the r-z plane need three points for their unique location. The mapping function also must equal unity at their node and zero at all others. Since the shape function has all these characteristics it is used as the mapping function 8 8 r = X>,-TV,-(£,»;) z = Y,ZiNlU,v) i=l 1 = 1 Numerical integration is required since the integration cannot be done exactly. Since the element has been mapped into £ and r) space and the element is isoparametric the limits of integration are -1 to +1. The numerical integration method used is Gaussian Quadrature which involves evaluating the formula at a number sampling points. Each evaluation is multiplied by a weighting value and the summation of all the sampling (A.24) Appendix A. Formulation of Thermal FEM Model 176 Figure A.88: Mapping of quadratic element from r and z space into £ and r? space Appendix A. Formulation of Thermal FEM Model 177 points is equal to the integral. The gaussian quadrature formula is: i+1 f (0 = Wif (6) + w2f (6) +... wnf (60 = E wj (fc) or for two dimensions: / / / (ft) = E E WiWi/ (fcifc) «=i j = i The location of the sampling points and corresponding weight coefficients are given in Table A.22. 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 we have 2n unknowns (/,• and £,•), and thus a polynomial of degree 2n-l could be constructed and exactly integrated. The corresponding error is of the order 0 (A2"). 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 (equation A.20, pCp NiNj^f). The function to be integrated is an incomplete 4th order-polynomial. The optimum number of sampling points is 3 (ie 2n — 1 = 4.-. 4 ^ = 2.5 m 3). For the two dimension the number of integration points are 3 in the T] direction and 3 in the £ direction, thus a total of 9 sampling points are used (32 = 9). The larger the number of sampling points the longer the solving times. A reduced number of sampling points can be used and still achieve the accuracy of the larger sam-pling point integration. Using fewer than optimum sampling points is knoAvn as reduced integration [21]. Reduced integration using 4 sampling points gives comparable results to nine point integration if stability criteria are satisfied. Stability criteria for reduced inte-gration requires that more than one element be present in the domain. Three sampling points are preferred when there are temperature dependent thermo-physical properties, sampling points allow the variation in conductivity, density, etc to be better modelled over the domain. For a model were the thermo-physical properties are constant reduced integration is used to achieve faster solving times. Appendix A. Formulation of Thermal FEM Model 178 Sampling Locations Sampling Points Weight Coefficients - 0.57735 02691 + 0.57735 02691 n = 2 1.00000 00000 1.00000 00000 - 0.77459 66692 0.00000 00000 + 0.77459 66692 n = 3 0.55555 55555 0.88888 88888 0.55555 55555 Table A.22: Sampling locations and weight coefficients of the Gaussian Quadrature For-mula A.5 Conduction Evaluating the conduction matrix for an element, [A']^, requires the following integra-tion: ONidNj dNidNj' + drdz dr dr dz dz This can be simplified by using matrices, or in terms of the integration formula: (A.25) [K]{e) = 2TT J r k [B]T [B] drdz w here [B] dN-dr ' dr dN, dN; dz ' dz dNn ' dr dNn ' dz Since the element is isoparametric the integration is simplified by having the limits of the integration set to —1 and +1 in £-77 space. It is necessary to evaluate and dr ' dz drdz in terms of £ and 77. Using the chain rule to get £ and 77 in terms of r and z gives: dNi dNi dr dNi dz + di dr d£ dz dt, or in matrix form dN, 1 dr dz { Mi ' dt dr dN, dr dz dNj dr, ) L dr, ' dr, . K dz dNj dr dNj dz Appendix A. Formulation of Thermal FEM Model 179 The above matrix, [J], is the Jacobian. The Jacobian in full form is: [J] = E n (A.26) 3Af,-3r; ' dr, J The Jacobian must be inverted in order to find the first derivatives of the shape functions in terms of r and z. The matrix inversion must be done numerically. 8N{ dr dNj dz [J] - l 8Nj dt dNj dr, Il\ > -^ 22 The inverse Jacobian is used to put the [B] matrix in terms of the £-77 coordinates dNj dt dNj dr, T dN, . T dNi T dN, , j dN, 12\-^r T -122-r dN, _ L T dN: 1 1u~dt ' 12 dt ' l i dr, ' Finally drdz in £-77 coordinates is j dN, , j dN, J 2 1 " 9 f + • /22-^ i T dNg , T dNg , his? + h2-§f dr dz — det [J] d( dr] Integration of the conduction matrix changes it from [Kc]{e) = 2icJJrk [B]T [B] drdz to [Jvc](e) = 2TT ^  ^ r fc [ £ ] T [B] <fe< [J] d£ dr? (A.27) The terms of equation A.27 are determined by gaussian quadrature numerical integration which gives: [/gW = £ £ W i Wj 2*r{M k [B{iitni)]T (A.28) i = l j = l ?^  = numbering of sampling points Wi and Wj = weighting constants £i and 77; = coordinates at sampling points Appendix A. Formulation of Thermal FEM Model 180 A.6 Convection Evaluating the convective matrix and vector for an element requires integrating the following two equations: Kff. = 2TT J^rhcNiNjdF fff = 2TT / rhc NiT^dT Jr or in terms of a single element [Kh] {e) = 2TT J rhc {N} {N}T dT {h} (e) =2TT ^ rhc {N} T^ dT In this case the integrations requires that dT be put into terms of r and z before the mapping into 77 space. Since (r,z) are orthogonal coordinates: dT = Vdz2 + dr2 or for derivatives dT _ (dzV (drs \dv) \<fy, Mapping into the 77 space gives dT _ dNn Km=zl dN„ dV J \ m = l dV To simplify notation the symbol A will be used, where A \m=l dN„ dr) + E \m=l dN, drj Appendix A. Formulation of Thermal FEM Model 181 Thus dT = Adr] The integration to be conducted now becomes [Kh]{e) = 2TT f rhc {N} {N}T A drj {A}{e) =2TT J^rhc {N} Adr? Using gaussian quadrature for numerical integration gives [Kh}[e) = E ^ H ) h< K*)} {N(n,)}T A (A-29) 1=1 {/ /»} ( e ) = E Wi 2^^,.) n c { 7 V ( n 0 } Too A (A.30) A.7 Time Stepping Techniques For transient heat transfer problem, equation A.21 must be solved. [K]{T} + [C]{T} = {f} (A.31) Solving equation A.31 is generally difficult analytically and cannot be done for non-linear situations. Thus the solution technique that is used is a trial function (finite element) discretization in the time domain. This method is known as a recurrence technique since each time step is solved using the previous time steps temperatures as initial conditions. The capacitance matrix is evaluated before the time stepping technique is discussed further. Appendix A. Formulation of Thermal FEM Model 182 A.7.1 Capacitance Matrix The shape functions used are the same as the conduction matrix, equation A.23. The Jacobian, [J] is as given in equation A.26. The capacitance matrix is C\f = 2*JJrpCp NtNjdr dz [C]{e) = 2ir JJrpCp {N} {N}T dr dz or in s and t coordinates [C] ( e ) =2TT ^ J' rpCp {N} {N}T det[J]d(dri or in terms of gaussian quadrature [C]{e) = 2TTEE Wi W3 r^.) p Cp {%„.)} {%,„)}r det [J] i=i j=i A.7.2 Time Stepping Using Recurrence Technique Equation A.21 is solved using the weighted residual approach. The number of trial functions used, will affect the accuracy and stability of the solution. In general the greater the number of shape (trial) functions used the better the accuracy and stability [28]. A two point (linear) scheme was used for the initial solution with a three point (quadratic) scheme used for all subsequent calculations. A.7.3 Two Point Recurrence Scheme Temperature is also approximated by another shape function. {T} = J2Nl{Tl} i {Ti} is the temperature vector at time /. The shape function, Ni, is a function of time only. Appendix A. Formulation of Thermal FEM Model n n + 1 •{ = t/At * At Figure A.89: Linear shape function for time Appendix A. Formulation of Thermal FEM Model 184 The shape function used will be linear as shown in Figure A.89. The shape functions are o < C < i ( = i-t Nn = l-(; Nn = -tt (A-32) Nn+1 = C; Nn+i = Thus {T} = {T„} Nn + {Tn + 1} JV n + i (A.33) Substituting equation A.33 into equation A.31 will give a residual which is minimized using the weighted residual method. The weighted residual for the time domain is f1 W3 {[C] ({Tn} Nn + {Tn+1} Nn+1) + [K] ({Tn} Nn + {Tn+1} JVB + 1) - {/} } d( = 0 Jo Collecting the T„ and Tn+1 terms gives {[A-] £ Wj Nnd( + [C] £ W, Nnd(} T„+ {[K]£wjNn+1d( + [C] £WjNn+1d(} T „ + 1 - £ Wj {/} d( = 0 (A.34) Substituting the values for the shape functions from equations A.32 into equation A.34 gives {[jin^^d-cK + i c i / ; ^ ^ ) * } ^ {[K] £ W, (d( + [C] £ Wj (™) rfc} T„ + , - £ Wj {/} d( = 0 Divide all terms by / Wjd(. To simplify the equation a new variable is introduced, 0. e _ Jo W (d( So W d( Thus the equation becomes [K](l-6) + [C] (-^)}Tn+ Appendix A. Formulation of Thermal FEM Model 185 {[K] e + [C] ( ^ ) } T n + 1 - {/} = 0 (A.35) where T r y = Jo1 ^  {/} <*C The same interpolation can be applied to the vector / that has been applied to the unknown vector Tn+1. {/} = {/}n + 1 0 + { / } n U - 0 ) For the convective boundary condition the force vector does not change significant!}' during a time step thus {/}„ and {/}n+1 are equal and {/} = {/}. Using 9 = | gives the best solution [22]. Thus the final form of the matrix equation to be solved is: [K] {2r"+13+ + [C] {^^} = 17) A.7.4 Three Point Recurrence Scheme The shape function used is quadratic as shown in Figure A.90. The shape functions are o < C < i C = ii Arn_ 1 = C ( i + 0|; N»-i = {\ + ()£i Nn = (i-C)(i + C); Nn = -^ (A-36) Nn+1 = -c(i-oh Nn+1 = (-i + c) ± Thus {T} = {r„_x} Nn^ + {Tn} Nn + {Tn+1} Nn+1 (A.37) Substituting Equation A.37 into Equation A.31 will give a residual which is minimized using the weighted residual method. The weighted residual for the time domain is Wj [C] ({rn_!} Arn_! + {Tn} Nn + {Tn+1} Nn+1) Appendix A. Formulation of Thermal FEM Model 186 Figure A.90: Quadratic shape function for time Appendix A. Formulation of Thermal FEM Model 187 + [K] ({Tn} Nn + {Tn+1} Nn+1) - {/} d( = 0 (A.38) Using the same method as in deriving the two point scheme, Equation A.38 can be written as: {[/<"] At 2 Q + 3 - 7) + [C] At ( 1 - 7 ) | r„_ 1 + {[K] At2 Q - 28 + 7) + [C] At (1 - 27)j Tn+ {[K] At2/? + [C] At 7 } T „ + a - {/} At 2 = 0 (A.39) w here _ IliW (c + i) rfc /-i W {/} <*c The same interpolation can be applied to the vector / that has been applied to the unknown vector Tn+\. If) = {/}„+! 0 + (/In (l - 2/? + 7) + {/}„_, Q + 0 - 7 ) For the convective boundary condition the force vector does not change significantly during a time step thus { /}„ , {f}n_1 and {/}n+1 are equal and {/} = {/}. The value of 7 and Q used are for the Dupont form with 7 = 1 and 3 — f . The Dupont form of the equation has been successfully used in the past [28]. The final form of the equation A.39 using the Dupont variables is: [Al { 3 r" + 1 4 + r- 1} + [C] {^f^} = 17} Appendix A. Formulation of Thermal FEM Model 188 A.8 Phase Transformation Solution Heat flow involving solidification or melting is difficult to model. The problem lies in accurately accounting for the latent heat of sohdification. If the latent heat is not properly applied during any time step, the resulting temperature fields will be incorrect. Initially, the latent heat from a phase change was included in the heat capacity by in-flating its value between the solidus and liquidus temperatures. The change implemented in the specific heat at a phase change is shown in Figure A.91. Using this method to account for the latent heat release leads to a number of errors. One error is due to the discontinuous slope in the specific heat, Cp, which is difficult to numerically integrate over an element. The other source of error occurs when the solidus and liquidus are between the same node in an element. This results in the latent heat not being included during that time step. To avoid the errors, the enthalpy [21,23,24] is employed to evaluate the heat capacity of the material. The enthalpy, which is the integral of the specific heat, is smooth through the sohdification temperature range ( Figure A.91). The enthalpy of a material can be expressed as the sum of the specific and latent heats. The latent heat evolved is a function of the local fraction of liquid [25]. An alternate means of incorporating latent heat considers the segregation as well. In an alloy system that is solidifying with complete mixing in the liquid, the temperature in the transition region changes during growth due to nonlinear solute segregation, as given by the Schiel equation. This method is superior to the Enthalpy method since it allows The heat transfer equation using enthalpy is: kV2T = p OH ~dt (A.40) Appendix A. Formulation of Thermal FEM Model enthalpy / heat capacity Ts Tl Temperature Figure A.91: Variation of heat capacity and enthalpy through a phase change Appendix A. Formulation of Thermal FEM Model 190 the use of the Scheil equation which predicts the latent heat release more accurately. For this method the enthalpy is written as H= [T (Cp + f(T)*L)dT (A.41) JTr where / (T) =fraction liquid L =latent heat of solidification Equation A.41 can be rewritten to separate out the specific heat term. H = [TCpdT+ fT f (T) * LdT JTr JTr The heat transfer equation can now be rewritten as: kV2T = p^-( [TCpdT + fT f (T) * LdT Ot \JTr JTr or k v 2 T = ' c > f t + J t f j ! ( T ) * L i T ) With the equation in this form the latent heat can be added on to the normal heat transfer equation. Evaluating the latent heat term from the heat transfer equation ( g j ? r ( / ( r ) * £ d r ) ) gives: dflt Applying Galerkins method to the heat equation with the latent heat included gives an extra term that must be included into the final matrix form of the equation: [KL] = 2icJ J rp^LNiNjdrdz The same interpolation as used for the temperatures has been used for the fraction liquid. / = E Ni (r,z) ft (t) (A.42) i = l Appendix A. Formulation of Thermal FEM Model 191 = 0 ^ ^ ^ ^ ^ — \ df/dT = max > Figure A.92: Variation of fraction liquid / with temperature contours (interface shape) dT ^ Jl dT For two or more dimensions the value of is dependent on the shape of the interface (Figure A.92). The variation in fraction liquid will be zero in the direction tangent to the interface; it is maximum in the direction normal to the interface. Hence the component of the principle fraction liquid in any direction will be the projection of the maximum fraction liquidus onto the direction of interest. If r is the normal direction to the interface Appendix A. Formulation of Thermal FEM Model 192 than the variation in fraction liquid is: dT dr dT df_ jdT_ dr j dr Since r and z the are orthogonal coordinates, the enthalpy variation in the r direction can be put in terms of r and z coordinates dr \\drj \dz J The same operation that used on the enthalpy are also apply to the variation of temper-ature in the r direction. The resulting equation is: a / dT (9T\2 , /"dT\2 \[d7) + \lTz) } (A.43) Equation A.43 is very similar to equation A.25. The latent heat equation can be manip-ulated into the [B] matrix form as in section A.5. Therefore equation A.43 becomes: d£ dT {f}T[B]T[B]{f} (A.44) .{T}T[B}T[B}{T}J Where{/} are the fraction liquidus values at the nodes and {T} are the nodal tem-peratures. The latent heat matrix that must be added onto the Galerkin temperature formulation is [KL]-{f}T[B]T[B]{f} [KLp = 2TT £ £ Wi W3 r^.) p t t i f y ) {T}T[BY[B} {T}, •'=i j=i det [J] (A.45) A.9 Solidification Gap Since the volume shrinkage of Cdi-^Zn^Te on cooling is larger than quartz, by one order of magnitude between 600 K and 1364 K, the gap will form between the metal and quartz Appendix A. Formulation of Thermal FEM Model 193 on cooling. This will change the heat flux between the quartz and the metal. Before solidification there is good contact between the melt and the crucible, with conduction as the mode of heat transport. After solidification and gap formation the mode of heat exchange changes from conduction to radiation. The solidification gap is modelled as if it were a convective/radiative boundary, with the ambient temperature given by the temperature of the adjacent boundary (Figure A.93). The estimated shrinkage gap is small (/im), as a result nodes at each side of the gap occupy the same r-z space. The convective matrix and vector for the right side boundary are: The value for T/e/t is the temperature of the left side boundary. It is evaluated at the left side elements gaussian integration points that are directly opposite to the right side element gaussian integration points. The left side of the boundary is evaluated as the right side is: Before solidification the value of the heat transfer coefficient is inflated to approximate conduction. A value of 8000.0 J / m 2 s K for the heat transfer coefficient across the gap before solidification was found to approximate conduction [30]. Between the liquidus and solidus temperature the inflated heat transfer coefficient is ramped down to the normal value. Using this method the quartz and CdTe have are two seperate temperature models which are linked to each other by the solidifaction gap boundary. Appendix A. Formulation of Thermal FEM Model 194 Figure A.93: Boundary elements that are a solidification gap Appendix B Formulation of F E M for Stress Fields The stresses and strains that occur in axi-symmetric coordinates are shown in Fig-ure B.94. The stress/strain calculations are based on the displacement that occurs to the domain that is being modelled. Using the same quadratic isoparametric element as used in the temperature analysis there are eight nodes, each of which has a displacement ({a}). The total combination of possible nodal displacements can be described in terms of the r and z coordinates. The r and z direction dispacement are known as u and v displacement; thus each node has two variables. For example, the displacement ({a}) of -node i consists of the displacement in the u and v direction. Each node in the element has two variables (u and v displacement). The isoparametric node has 16 variables instead of the 8 that are used for temperature analysis. Figure B.95 shows the quadratic isoparametric element and the two variables per node. If the nodal 195 Appendix B. Formulation of FEM for Stress Fields Figure B.94: Strain and Stress involved in the analysis of axi-symmetric solids Appendix B. Formulation of FEM for Stress Fields Figure B.95: Quadratic isoparametric element for Appendix B. Formulation of FEM for Stress Fields 198 displacement is represented by the shape function approximation the vector form is: actual {N} {a}'I NlUl N2u2 N2v2 (B.46) N&us N8v8 The relationship between displacement, strain and stress will now be shown. The strain has four components and is related to the general displacement by: (4 dv dz < > = < du dr u r . 1" . du _ i _ dv { dz ^ dr . (B.47) Using the displacement function defined by equation B.46 the relationship between the strain and displacement vector is: {e} = [B] {a} in w hich dN, dz dN, dr 8N2 dr 2 , dN2 dz 0 0 dNn dr 1N, dNs ' dz 0 0 Appendix B. Formulation of FEM for Stress Fields 199 The relationship between stress and strain is given by: {*} = { } = [D] ({6} - {*}) + M OS were {eo} is the thermal (initial) strain, {cr0} is the initial stress and [D] is the material properties matrix which is given by [D] = E(l - u) 1/ 1 V 1 - 1 / ' 1 ' 1 - 1 / ' 1 - 1 / ' 1 - 1 / ' ' 0 , 0 , 0 0 0 0 l - 2 i / ' 2(1-1/) J .E is Youngs Modulas and v is Poissons ratio. Having the relationship between the dis-placement, strain and stress the general solution for stress analysis will now be considered. Solving for thermal stress requires that the Variational Principle be used. The starting point for the variation method is the total potential energy of a structure ( n p ) . The absolute minimum value of n p corresponds to the exact solution. The potential energy of a structure ( n p ) is: Hp = /fl Q (e} T [D] {e} - {ef [D] {e0} + {ef {a}) dQ ~ [ M T {F} d$l - I { « } T {<^ } dT - {a}T {P} (B.48) Jn Jr where {a} = the displacement field {e} = the strain field [D] = the material properties matrix {eo} , {o~o} — initial strains and inital stresses Appendix B. Formulation of FEM for Stress Fields 200 {F} — body forces {</>}= surface tractions {P} = loads applied to nodes by external agencies For thermal stress there are no surface tractions, initial stresses, body forces or applied loads. Thus equation B.48 becomes: Hp = / n Q Of [D] {e} - {ef [D] {£o}) <Kl (B.49) Putting equation B.49 in terms of the displacement gives: Bp = JQ (I {af [B)T [D] [B] {a} - {af [B]T [D] {e0}) dft (B.50) The displacements are nodal values thus they can be removed from the integration. n p = {af I \ [Bf [D] [B] da {a} - {af I {af [B]T [D] {e0} dtt (B.51) JQ I JQ The solution to the stuctural problem is the actual displacement {<x}actua; which makes Ftp stationary to with respect to any small variations 6 {a}actual. Thus the solution to the structure problem is the variation flip = 0 The displacement has been approximated in the potential energy by a the shape functions: 8 {^actual ~ E a>N< .'=1 The variation of lip in terms of the shape functions is d n p r d n P c d n P c oa oax oa2 This being true for any variation 8a yields a set of equations dILt (B.52) da i da = < > = 0 Appendix B. Formulation of FEM for Stress Fields 201 Thus by solving for the value of the nodal displacements can be found. Differentiating equation B.51 with respect to a gives: I k = L \[B]T [D] [B] d^{a}+{a}T L \[B]T [D] [B] dn~ I [B]T[D] {e0}dn = 0 (B.53) Jn Using the following definition / \[B}T[D][B]d£l{a}^{a}T ! \[B]T [D][B)d£l Jn l Jn I allows equation B.53 to become ^ = / [Bf [D] [B] dn {a} - j [B]T [D] {e0} dQ. = 0 (B.54) oa Jn Jn Thus to determine the thermal stresses equation B.54 must be solved. Putting the integrals in terms of the radial coordinates gives 2TT J J r [B)T [D] [B] dr dz {a} - 2TT J J r [B}T [D] {e0} dr dz = 0 (B.55) Equation B.55 can be put in terms of a matrix equation [K] {a} - {/} = 0 where [K] = 2TT J J r [B]T [D] [B] dr dz and {/} = 2TT fJr[B]T [D}{e0}drdz As in section A.5 the integration must be done numerically and the element must be mapped into £ and 77 space. The resulting equation that is nummerically integrated and solved is: [K] = 2TT £ £ WiWjr^ [D] [B^)} det [J] (B.56) i=i j=i Appendix B. Formulation of FEM for Stress Fields 202 and [/] = 2TT E E WiWjr^ [%.try)]T [D] {e0(^.)} det [J] (B.57) where [B] in £ and 77 space is given as iVi r a/v, , T a/v, ^21-9^ + ^22-9^ J2i-gf + J22-9+ 0 0 Jii-flF- -r J12-7- a/v. , r a/Vg , Jll 9^  + J 12-9^ T dNg , 7 a/v* , i21 9^  + J 2 2-^f 0 0 7 a/v* , r a/Vo , -'li-aT + J i 2 - 9 ^ d'Z 1 ± r i dr, ' ••• ' J21-9£- T *22-QZ > ^dZ where n is the number of sampling points used for gaussian integration and W are the weightings applied at the sampling locations (£;,r?,-). The equation for the Jacobian [J] is given is section A.5. The thermal strain for an isotropic material are given as > = < e<?o 0 V J where a is the thermal expansion coefficient and A T is the temperature change that is causing the thermal stress. The temperature change are calculated at each gaussian integration point. Appendix C Grashof Number Calculation The Grashof number (Gr) is the ratio of the buoyancy force to shear force. Buoj'ancy is the driving force for natural convective heat flow. The Grashof number is: p2 g (3 ATP3 Gr = p? p = fluid density p, = fluid viscosity (3 = fluid coefficient of thermal expansion — 1 dp ~ p dT g = acceleration of gravity D — characteristic dimension of body A T = temperature difference between body and ambient fluid C l Grashof Number for Tin A temperature change (AT) of 4 degrees across the fluid was used to calculate the Gr number. A fluid temperature of 280°C was used as the mean fluid temperature to get the thermophysical properties [39]. The thermophysical properties are: (C.58) 203 Appendix C. Grashof Number Calculation 204 p = 7.00 g/cm3 = 1.678 cp at 280°C p = 9.02 x 10-5 1/deg: 9 = 981 cm/s2 D = 1 cm AT = 4.0 degrees The Gr calcualtion is: G r _ (7.00)2 g 2/cm 6 981 cm/s2 9.02 x 10~5 1/degree 4 degrees (l) 3 cm3 (0.01678)2 g 2/cm 2 sec2 The grashof number number for tin at 280°C with a 4 degree temperature difference in the fluid is 6.16 x 104. Since different systems are being examined it is better to have both at the same homologous temperature (TJJ). The Tjj for the Gr calculation is: TH = T/Tm _ 280 231 = 1.212 C . 2 Grashof Number for CdTe No CdTe data readily exists, therefore the Cd and Te values will be averaged using the same source for the thermophysical properties as for tin [39]. The mean fluid temperature of the values is determined from the homologous temperature: T = THxTm = 1.212 x 1364K = 1653 K The thermophysical properties are: Appendix C. Grashof Number Calculation 205 p = \ (8.0 + 5.8) g/cm3 p = the liquid conductivity of Te cannot be found /? = 1.31 x 10~4 1/degree g = 981 cm/s2 D = 1 cm A T = 4.0 degrees Since the viscosity of Te cannot be found the theoretical viscosity of CdTe will be calculated [41]. C.2.1 Theoretical CdTe viscosity calculation In order to calculate the theoretical viscosity of a liquid metal the following properties must be known: Appendix C Grashof Number Calculation 206 M — molecular weight = 120 g/mole 8 = atomic radius 8 = 1(3.2 + 3.08) °A = 3.14 A Tm = melting point = 1364 K T = fluid temperature = 1653 K n = particle number density = 6.9g/cm3 x x 6.02 x 1023 atoms/mole x (lO"8)3 cm 3/A 3 = 0.0346 atoms/A3 R = gas constant = 8.28 x 107g cm2 / s2 mole K Determine energy parameter The energy parameter (|) must be calculated, f = 5.20 xTm = 5.20 x 1364 = 7092.8 Determine reduced temperature(T*) Using the mean fluid temperature the reduced temperature is: Appendix C. Grashof Number Calculation 207 Figure C.96: The reduced viscosities of liquid metals and their dependence on reduced temperature T* = f x r _ 1653 7092.8 = 0.233 K Thus using 1/T* and Figure C.96 gives a reduced viscosity function (u* (V*)2) of 2.60 Appendix C. Grashof Number Calculation 208 Calculate reduced volume (V*) The reduced volume is: V* V* _n 6s 1 0.0346 atoms/A3 x (3.14A) The reduced volume is V* = 0.9331. Using the reduced volume and the reduced viscosity function gives: . 2.6 (0.9331)2 H* = 2.99 Theretical Viscosity The theoretical viscosity as related to the reduced viscosity is: _ u* x (M x R x T)* M ~ S2 x N0 _ 2-99 1/atoms x ( 120 g/mole x 8.28 x 107g cm2/s2 mole K x 1653 K y (3.14A)2 (l0- 8cm/A) 2 x 6.02 x 1023atoms/mole Thus the viscosity is 2.04 x 10 - 2 g/cm sec or 2.4 centipoise for CdTe at 1653 K. C.2.2 Grashof number _ (6-9)2 g 2/cm 6 981 cm/s2 1.31 x 10~4 1/degree 4 degrees (l) 3 cm3 (0.0204)2 g 2/cm 2 sec2 The Gr number for CdTe at 1563 K with a temperature change of 4 degrees is 5.88 x 104. Appendix D Relative Thermal Resitance The thermal circuit of Vertical Bridgman growth of CdTe is shown in Figure D.97. There is radiation that is transmitted through the quartz to the CdTe (Ri), radiation that is absorbed by the quarts (R2) and forced convection along the quartz ( A 3 ) . The thermal resistances for the boundary conditions of CdTe are worked out assuming a furnace temperature of 1380 K and a quartz and CdTe temperature of 1380 K. Resistances are as follows: radiation convection conduction were the radiative heat transfer coefficient with a temperature of 1380 K and 1370 K is: hr = a e (r2 + T22) (2\ + T2) = 5.67 x 10~8 e (l3 802 + 13 702) (1380 + 13 70) = e590W//m2A' R R R hrA 1 hc A L The resistances were calculated with the following values: quartz: e = 0.2 209 Appendix D. Relative Thermal Resitance 210 CdTe R4 A A / V (radiation to CdTe) R4 R l A A / V (radiation transmitted) (radiation to CdTe) R3 R2 AAA/ V v V (conduction) R3 (conduction) (radiation to quatz) R5 Radiation Transmitted Radiation A A A / W V / V W -Radiation Absorbed Forced (convection to quartz) Convection Figure D.97: Thermal circuit of Vertical Bridgman boundary conditions. Appendix D. Relative Thermal Resitance 211 cross section forced convection k = 2.13W/Km L = 0.00275m A = lm2 hc = 50W/m2K The calculated resistance are: hr A 0.8 x 590 ' R2 = = n 0 * n = Q.QQ85K/W hr A 0.2 x 590 ' ft = O = S i = a 0 0 2 3 / ^ Figure D.98 gives the thermal circuit with the resistances. Examining Figure D.98 it is evident that the radiation from the furnace that is transmitted through the quartz will dominate all other boundary conditions. The resistance of this mode of heat transfer (0.0021 K/W) is an order of magnitude lower than the other resistnces (0.0115 K/W). Appendix D. Relative Thermal Resitance 212 0.0021 K/W - A A A r CdTe 0.0115 K/W — W V 0.023 K/W Radiation Transmitted Radiation Absorbed Radiation Forced Convection Figure D.98: Thermal circuit of Vertical Bridgman boundary conditions with resistances. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0078458/manifest

Comment

Related Items