Mathematical modeling of thermal processes during “cassette” crystallization of chalcogenides

A new, relatively simple and highly efficient modification of the directional melt crystallization method in the form of a multi-cassette process has been considered. This study is based on Russian Patents and technological studies conducted at National Research and Technological University MISiS. As a result, mathematical models of a multi-cassette method have been developed for 3D radiation and conduction analysis of thermal processes in the entire volume of the heating unit and 2D analysis of convection and conduction heat exchange in a separate cassette. Parameters have been calculated on the basis of these mathematical models for clarifying the effect of heating unit component arrangement and dimensions on the formation of thermal fields in cassette units, the effect of vertical homogeneity of heat supply to the cassette unit and heating power reduction rate during crystallization on the shape of the crystallization front, as well as the effect of small asymmetry in cassette design and violation of cassette bottom cooling homogeneity on convection and asymmetrical heat transfer. Application of the conductive and radiative heat exchange model to the entire heating unit has allowed us to calculate process parameters on the basis of which we have analyzed the effect of heating unit components, their arrangement and temperature on the heat exchange conditions at the cassette unit boundaries. Application of the convective and conductive model to one growth cassette has shown that asymmetrical design and boundary thermal conditions as well as unstable vertical temperature gradient lead to the formation of convection vortices and substantial crystallization front deviation from planar shape. Calculations on the basis of the convective mass exchange model have shown that an increase in the crystallization rate by one order of magnitude greatly increases the tellurium flow into the crystal thus substantially altering the melt composition in the vicinity of the crystallization front and hence serving as a potential origin of dendrite growth. The authenticity of the calculation results has been verified in a number of tests aimed at analyzing the effect of heat and mass transport on crystallization front shape for cassette cooling rates that are typical of polycrystalline bismuth telluride growth processes.


Introduction
Nowadays a number of bulk thermoelectric material synthesis methods have been developed. Most of them include powder compaction using various methods e.g. synthesis of ultrafine powders by gas condensation in inert gas atmospheres [1] or using a plasma-chemical method [2], chemical synthesis [3] followed by spark plasma sintering [4] as well as powder grinding in a ball mill [5]. Still there are problems in the development of these methods originating from residual porosity after compaction, specimen contamination during powder preparation or consolidation and fabricating specimens of greater dimensions.
The objective of intense plastic deformation methods [6] is to produce nanostructures in bulk specimens and workpieces by refining their microstructure to a nanoscale level. A wafer casting method has been suggested [7] wherein the opposite faces of the as-cast wafer of an A 5 B 6 material are parallel. These wafers have a layered structure forming at least two matrices of cleavage planes that are misoriented relative to each other in a way that the cleavage planes of one matrix are tilted relative to the cleavage planes of the other matrix and to the base wafer surfaces. The presence of two mutually misoriented cleavage plane matrices in the as-cast wafer structure raises problems with wafer cutting into rectangular dice since one does not know the actual cutting plane orientation relative to the two cleavage plane matrices and to the base wafer surfaces.
The aim of this work is to modify the Bridgman method by introducing cassette wafer crystallization. This directional crystallization method has already been described [8] for the growth of crystalline wafers of bismuth telluride solid solutions. As-grown wafers are then cut into dice in a direction perpendicular to their base surfaces. The process provides for high mechanical strength but produces considerable misorientation between cleavage planes. This method has been further developed [9] to by directional crystallization growth of crystalline planes with a more perfect crystalline structure with lower cleavage plane misorientation angles due to a more efficient cleavage plane orientation control at crystal nucleation stage and during further growth.
The applied aspects of the study were as follows: -search for science-based approaches to growth unit design improvement and choice of higher cooling rates providing for dendrite-free crystallization of bismuth telluride melt; -identification of potential negative thermal factors associated with radiation and conduction heat exchange in the growth unit volume and convection and conduction heat exchange in separate growth cassettes that affect crystallization, thermoelectric properties and strength of as-grown bismuth telluride wafers.
Empirical studies of wafer growth, structure, thermoelectric properties and strength were carried out earlier [10]. Mathematical modeling [11] was used for compara-tive analysis of two basic technological approaches to the synthesis of bismuth telluride based thermoelectric materials: equal channel angular extrusion [12] and Bridgman cassette crystallization. The effect of design features and temperature growth mode on plastic formation and crystallization was analyzed.
Mathematical modeling of crystallization was carried out [13] taking into account of phase transformations occurring in accordance with the phase diagram of the multicomponent material in its liquid and solid phases. It was assumed that the permeability of the two-phase zone is isotropic and depends not only on the porosity but also on the geometry of the porous material. For example, introduction of an empirical parameter allowed presetting the distance between dendrite branches and the well-known Kozeny-Carman formula was used based on similarities between the porous material and a system of parallel permeable capillary channels [14].
The mathematical approach used in this work is based on the mathematical modeling of thermal processes occurring during cassette crystallization. Mathematical models were developed for this purpose. One of them takes into account the specific features of thermal processes in the entire volume of the cassette unit while others incorporate thermal and mass exchange processes in separate cassettes.
Schematic of the heating unit for the cassette growth method [8,9] is shown in Fig. 1. There the water-cooled chamber 1 contains longitudinally arranged resistive heater plates 2. The top chamber 3 located between the heater plates is used for melting and convective homogenization of the charge mixture and is connected with the chamber 4 containing the crystallization cassettes. The melt is poured a single time through the tubes 5 from the melting chamber into the crystallization cassettes installed on the massive cooled metallic plate 6 resting on the bottom of the water-cooled chamber 1. After melt is poured a vertical temperature gradient is produced in the crystallization cassettes and crystallization starts in the cassettes due to resistive heater power reduction. The heater power reduction rate determines the crystallization rate and is to be optimized for the required crystal wafer quality. This mathematical model which takes into account the complete geometry of the heating unit is useful for the calculation of thermal boundary conditions that can be preset in another simpler geometrical model for each separate cassette. Schematic of the model for thermal and mass exchange processes in separate cassettes is presented in Fig. 2. The cassette comprises the graphite case I having a number of cavities: the main cavity containing the crystallizing melt II, the crystal seed cavity III, the melt pouring channels IV and the mounting holes for cassette assembly into units V. The overall dimensions of a single cassette (x × y × z) are 4.4 × 6 × 1.3 cm 3 .
Initially the melt fills a cassette between two narrow graphite plates along which a temperature gradient is generated. The advantage of this method is the possibility of simultaneously providing a large number of gaps between the graphite plates so as to obtain a large number of thermoelectric material wafers in a single crystallization process. The polycrystals grown using this method have higher mechanical strength than single crystals, their electrophysical properties being close to those of single crystals [9]. The process provides for crystallization rate control through temperature gradient adjustment which in turn determines the shape of the crystallization front. Experiments have shown that for a flat crystallization front the growing texture is the most efficient from the viewpoint of bismuth telluride structural anisotropy. However significant distortion of the crystallization front causes large misorientation between as-grown polycrystal grains.
We present below our results of development of a theoretical and methodical approach to analysis of heat and mass transport during cassette crystallization. This approach was implemented in three mathematical models used in conjugation for the key input and output parameters: conduction and radiation heat exchange for the entire heating unit, conduction and convection heat exchange in the growth cassette and convection mass exchange in the two-component bismuth and tellurium melt. These mathematical models were software-implemented on the basis of the finite difference method and the finite element method on the CrystmoNet software platform [15].

Modeling of thermal processes in heating unit
The 3D model of the growth unit was developed on the basis of the finite element method. Schematic and the main structural components of the growth unit are shown in Fig. 3.
The arrangement of the heating unit structural components is quite complex. Radiation heating of the cassette unit by the resistive heater located in the heating unit is of great importance. The heating unit typically contains a number of materials having different thermophysical properties (Table 1). Data on the thermophysical properties were borrowed from earlier works [16][17][18]. The solidus T s = 863 K and liquidus T l = 865 K temperatures were chosen on the basis of the Bi 2 Te 3 -Sb 2 Te 3 phase diagram [19].
The heater plates 3 radiate heat toward the cassette unit 2. The heat flows down conductively toward the metallic plate under the cassettes. This heating unit is located in the closed chamber 1 whose outer case surface is maintained at room temperature through permanent water circulation in the case.  The overall dimensions of the heating unit components (x × y × z) are as follows: chamber 27 × 13.5 × 14 cm 3 ; heater 1 × 9 × 9 cm 3 ; cassette unit 4.4 × 6 × 10 cm 3 ; steel plate 8 × 0.5 × 14 cm 3 .
Numerical modeling of this design amounts to solving the heat transfer equation in the complex-shaped heating unit zone with allowance for different thermophysical properties of heating unit components and in conjugation with calculation of radiation heat exchange and crystallization processes.
The heat transfer equation has the following form: (1) Here the index i shows the respective component of the heat unit ( Fig. 3 and calculation parameters in the Table 1). It is hypothesized that the melt crystallizes in the range between the liquidus T l and solidus T s temperatures. In this range (ΔТ = T l -T s ) the crystallization front is a two-phase zone containing the crystalline phase and the melt. This is modeled by introducing an additional heat source in the crystallizing material along with the volume heat source in the resistive heater Q H into the heat conductivity equation Eq. (1): Here L = 1.35 × 10 5 J/kg is the latent melting heat and φ is the volume fraction of the solid phase in the two-phase zone which is determined by the crystal growth kinetics. However, melt undercooling in this case is but moderate and φ can be determined from the equilibrium phase diagram of the melt with allowance for the identity transformation: Practical calculations are carried out using The boundary conditions at the solid surfaces are set as follows. It is hypothesized that the outer surface of the heating unit case meets the condition T = T c where T c is the temperature of the water-cooling circuit which is usually maintained at room temperature 300 K. The common contact surface of two materials including those with different heat conductivities meets the heat flux balance condition.
The open inner surfaces of the heating unit are incorporated into radiation cuvettes for which the following heat balance condition is met: where it is assumed that two heat exchange mechanisms are implemented on the surface k having the temperature T k and belonging to the material with the heat conductivity λ k : -convection heat transfer with the surface heat exchange coefficient α from the environment having the temperature T c ; -radiation heat exchange of this surface k with the "visible" surrounding surfaces of the radiation cuvette.
Here ε k is the thermal emissivity of the surface k, σ is the Stefan-Boltzmann constant and T e,k is the effective radiation environment temperature for this surface.
The calculation results allowed us to identify specific features of the thermal processes, e.g. to evaluate the effect of design factors (mutual arrangement and dimensions of heating unit components). Now we consider specific features of the thermal field in the beginning and in the end of growth.
In accordance with the schematic illustrated in Fig. 2 the melt initially fills the cavities II, III and IV under the conditions of the preset vertical temperature gradient with the lowest temperature at the bottom and the highest temperature at the top of the cassette. The cassette unit was cooled down by reducing the power of the heater. Figure 4 shows the thermal field in the cassette unit. The cassette unit was arranged so heat radiation from the heater was directed toward its side surface (y × z). The butt surfaces (x × y) were oriented toward the water-cooled chamber wall. Therefore, the central part of the cassette unit was heated to a lower degree and this can be seen from the pronounced upward bent of the isothermal curves in the center. It can also be seen that this arrangement of the cassette unit develops inhomogeneous thermal conditions for the cassettes in the center and at the periphery of the cassette unit. This is confirmed in practice by different quality of wafers grown in central and peripheral cassettes [20].
The process starts with pouring the melt into the cassettes at 900 K. The initial heater temperature is 1300 K corresponding to the heater power Q H = 9 kW. Figure 5a shows isothermal curves in the heating unit and in the cassette unit. Analysis of these isothermal curves shows that the central part of the cassette is more heated compared with the peripheral and bottom parts. Cassette bottom Notations: ρ is density; λ is heat conductivity; С p is heat capacity; ε is thermal emissivity; 1-4 are heating unit component numbers as per Fig. 3.
contact with the massive metallic plate resting on the water-cooled chamber bottom maintains the crystallization front (865 K) almost flat. This cassette cooling pattern does not allow melting of the seed in the plane III (Fig. 2). At the end of the process the heater temperature decreases to 1250 K corresponding to the heater power Q H = 7 kW. The 865 K isothermal curve corresponding to the crystallization front has a W-shape with a small downward bent at the crystal location (Fig. 5b).

Modeling of thermal processes in separate cassette
For a separate cassette as shown in Fig. 2 we assumed that by analogy with the heating unit calculation there is the intermediate phase between the crystal (solid phase) and the melt (liquid phase), i.e., the crystallization zone at the temperature T which is above the solidus temperature T s = 863 and below the liquidus temperature T l = 865 K. The heat conductivity equation (Eq. (1)) also takes into account latent crystallization heat release , where the volume fraction of the solid phase in the two-phase zone is set by the following linear relationship: The heat transfer equation for the melt is written with allowance for thermal convection in the following form: To determine the velocity vector V and the pressure p in the bismuth melt (Bi) we solve the Navier-Stokes equation and the continuity equation with allowance for gravity heat convection in the Boussinesq approximation: where g is the gravity vector, β Т is the thermal expansion coefficient and µ is the dynamical viscosity coefficient of the melt. We solve the tellurium (Te) transport equation jointly with Eqs. (5) and (6): where M = ρ l C is the Te concentration in the melt (kg/ m 3 ) and C is the relative mass of Te per 1 kg of melt. The following boundary conditions are imposed in the calculation upon the sought velocity and concentration distributions: at the top cassette boundary the melt motion is absent but the Te concentration is set as C = C eo ; at the side cassette surfaces the melt motion is absent and a zero Te flux is set; at the crystallization front the melt motion is absent and Te transport is taken into account in the following mass balance relationship: where R is the crystallization front velocity relative to the normal n to the crystallization front. Parameters for mass exchange calculation in the crystal-melt system are given below.

Mass exchange parameters for Bi-Te melt
Dynamical  ....................................................................... The problems of increasing the melt crystallization rate in the cassette are important since this is a key task for increasing the growth process output. However practical attempts to significantly accelerate the process lead to violation of the thermal balance in the growth cassette and hence large crystallization front distortion and dendrite growth. Ideally, cassette cooling should always maintain a stable vertical temperature gradient with a gradual temperature decline. This is required for avoiding intense convection in the melt and hence maintaining a flat or close to flat crystallization front.
The convection model was implemented for a 2D case. The boundary conditions were borrowed from the 3D radiation and conduction model. For the preset cassette unit heating conditions (Fig. 4) the thermal fields in a separate cassette were calculated for the beginning and the end of crystallization (Fig. 6). Figure 6a shows a prominent convexity of the crystallization front toward the melt which is to a large extent smoothened in the course of the process (Fig. 6b). A stable vertical temperature gradient is maintained throughout all the crystallization process stages (Fig. 2, temperature profiles 1 and 2). Thus, for a stable vertical temperature gradient and absence of convection and thermal asymmetry throughout the growth process, melt convection is suppressed almost completely and the isothermal curve distribution corresponds to the pattern for conduction heat exchange in the melt.
The situation however changes dramatically if there is some thermal asymmetry caused by significantly inhomogeneous cooling of the massive plate on which the cassette is mounted. Mathematical modeling showed the specific features of this thermal situation. If cooling occurs under these conditions the vortex shapes become different, with the greatest vortex circumventing the crystallization front at a relatively large velocity of ~0.034 cm/s (Fig. 7a) distorting the crystallization front (Fig. 7b, isothermal curve T = 865 K).
Another factor leading to convection fluxes in the melt is associated with inhomogeneous vertical heating of the cassette. Slow cooling corresponding to the crystallization front velocity R = 0.15 mm/min produces symmetrical vortex structures implementing weak convection mode (Fig. 8a) which provides for slight crystallization front convexity toward the melt (Fig. 8b).
However, there are attempts in the industry to accelerate cassette cooling by drastically reducing the heater power. In practice this causes considerable crystallization front distortion and dendrite growth. Calculations showed that this is entailed by an unstable vertical temperature gradient generating asymmetrical vortices. Upon fast   cooling (R = 1.2 mm/min) a small asymmetry of cassette design leads to asymmetrical vortices and isothermal curves (Fig. 9). Vortex asymmetry is caused by small design differences between the right-hand and the left-hand parts of the cassette. The emerging intense convection significantly changes the thermal field in the melt bulk and in the vicinity of the crystallization front.
At R = 6 mm/min melt convection makes the volume distribution of Te inhomogeneous and increases its flow into the crystal (Fig. 10). At T = 865 K the decrease in the Te concentration becomes noticeable and violates the melt composition required for the crystallization of Bi 2 Te 3 as per the phase diagram [19].
For comparison the radial Te distributions are shown at the crystallization front (see Fig. 11a). Their analysis suggest that only low crystallization rate e.g. R = 0.3 mm/ min provide for for a homogeneous radial Te distribution.
If R is far greater e.g. by one order of magnitude the radial variation in the Te concentration becomes tangible. Growth of the relative radial inhomogeneity of the Te concentration with an increase in R is shown in Fig. 11b. Thus, convection produces a large radial inhomogeneity in the tellurium distribution. Violation of the required melt composition in the vicinity of the crystallization front may lead to unstable formation of the crystalline bismuth telluride phase.

Conclusion
Mathematical models for modification of the Bridgman method in the form of cassette crystallization were developed.
Application of the conduction and radiation heat exchange model to the entire heating unit allowed us to calculate process parameters on the basis of which we analyzed the effect of heating unit components, their arrangement and temperature on the heat exchange conditions at the cassette unit boundaries. Application of the convection and conduction model to one growth cassette b Figure 9. (a) Vortex structures and (b) isothermal curves in the beginning of crystallization for unstable temperature stratification leading to intense asymmetrical vortices in melt during rapid cassette cooling at V = 1.2 mm/min. Figure 10. C/C eo Te concentration isocontours in cassette volume for stable temperature stratification and high crystallization rate R 3 = 6 mm/min. b Figure 11. Te concentration profiling: (a) radial C/C eo profiles along crystallization front for different crystallization rates (R 1 = 0.3, R 2 = 3 and R 3 = 6 mm/min); (b) maximum radial ∆C/C max inhomogeneity at crystallization front as a function of crystallization rate R.
showed that asymmetrical design and boundary thermal conditions as well as unstable vertical temperature gradient lead to the formation of convection vortices and substantial crystallization front deviation from planar shape.
Calculations on the basis of the convective mass exchange model showed that an increase in crystallization rate by one order of magnitude greatly increases the tellurium flow into the crystal thus substantially altering the melt composition in the vicinity of the crystallization front and hence serving as a potential origin of dendrite growth.
The authenticity of the calculation results was verified in a number of tests aimed at analyzing the effect of heat and mass transport on crystallization front shape for cassette cooling rates that are typical of polycrystalline bismuth telluride growth processes.