- Email: [email protected]
Phase field modeling of grain structure evolution during directional solidification of multi-crystalline silicon sheet
Journal of Crystal Growth 475 (2017) 150–157
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/crys
Phase field modeling of grain structure evolution during directional solidification of multi-crystalline silicon sheet H.K. Lin, C.W. Lan ⇑ Department of Chemical Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 7 March 2017 Received in revised form 26 May 2017 Accepted 16 June 2017 Available online 19 June 2017 Communicated by M. Plapp Keywords: A1. Computer simulation A1. Directional solidification A1. Planar defects A1. Nucleation B2. Semiconducting silicon
a b s t r a c t Evolution of grain structures and grain boundaries (GBs), especially the coincident site lattice GBs, during directional solidification of multi-crystalline silicon sheet are simulated by using a phase field model for the first time. Since the coincident site lattice GBs having lower mobility, tend to follow their own crystallographic directions despite thermal gradients, the anisotropic energy and mobility of GBs are considered in the model. Three basic interactions of GBs during solidification are examined and they are consistent with experiments. The twinning process for new grain formation is further added in the simulation by considering twin nucleation. The effect of initial distribution of GB types and grain orientations is also investigated for the twinning frequency and the evolution of grain size and GB types. Ó 2017 Published by Elsevier B.V.
1. Introduction Grain boundaries (GBs) of multi-crystalline silicon (mc-Si), grown by directional solidification, play a crucial role on the electrical properties of solar cells. Controlling GBs is often necessary for improving the ingot quality. For example, the dendritic casting method [1] is to initiate many highly symmetric coincident site lattice (CSL) GBs, e.g., R3 GBs, through grain growth at high undercooling generated along the crucible bottom [1], and they are useful for the high lifetime in the grown ingot. In ribbon growth, the twin grains have better lifetime and their formation have also been discussed by Stockmeier et al. [2] based on thermodynamics arguments. Nevertheless, the detailed GB interactions and how to control them are still not quite clear. Recently, the ingot growth has made a significant progress by initiating uniform and small grains at the crucible bottom, i.e., the so-called high-performance mc-Si (HPMC) [3]. The large amount of non-R GBs, generated at the beginning of solidification, would relax the thermal stress and thus reduce the generation of the dislocations; they also terminate the propagation of the dislocations during ingot growth [3–5]. However, during ingot growth, more R3 GBs appear while the proportion of non-R GBs decreases, and the ingot quality near the top is thus deteriorated. Therefore, understanding the evolution of grain structures is very important to further control the directional solidification of Si ingots or ribbons. ⇑ Corresponding author. E-mail address: [email protected] (C.W. Lan). http://dx.doi.org/10.1016/j.jcrysgro.2017.06.016 0022-0248/Ó 2017 Published by Elsevier B.V.
GB evolution during directional solidification of silicon has been paid much attention in recent years. Fujiwara et al. [6] found that the low-energy grain, such as h1 1 1i, tends to overgrow the others at low undercooling. Chen et al. [7] were the first to simulate this phenomenon by using a two-dimensional phase field model considering anisotropic interfacial energy and kinetic coefficient with some success. Duffar and Nadri [8] further described the twinning mechanism for the nucleation on the {1 1 1} facet plane near GBs. However, the critical undercooling proposed to form a twinned nucleus in their mechanism was 9 K or larger. Their model was recently modified by Lin and Lan [9], who considered the interactions between the nucleus and the neighbor grain at the groove. A more reasonable undercooling for the twin nucleation on the {1 1 1} facet plane was derived, i.e., less than 1 K. This value was more consistent with the experimental observation [10]. To confirm the twin formation, in-situ X-ray synchrotron imaging was used to investigate silicon crystal growth by Tandjaoui et al. [10]. They characterized the birth of new grains and showed that twins nucleated exactly on {1 1 1} facets at the GB grooves [11]. Moreover, Wong et al. [12] investigated the evolution of grain orientations and GBs by analyzing wafers at different heights from directional solidification ingots seeded with small and randomly oriented silicon beads. Indeed, twin nucleation was found at the tri-junctions. In addition to these research, there have been several works related to the interactions between R3n GBs [13–16]. Recently, the interactions between R3n and non-R GBs were further investigated by Lin et al. [17] for the directional solidification of mc-Si sheets at different speeds. At lower speeds, they found
151
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
that R3 GBs decreased due to the interactions between R3 GBs and non-R GBs, while twinning was a key mechanism for the increase of R3 GBs at higher speeds. Their study using mc-Si sheets unfolded the detailed mechanisms of GB interactions in silicon. However, up till now these have not yet been simulated successfully. Moreover, the anisotropic properties of GBs in silicon play a crucial role in the GB interactions, and thus in their evolution. The classic Read and Shockley dislocation model [18] describes the energy in low-angle non-R GBs. For high-angle GBs, the wellknown Brandon criteria [19] are used to calculate the maximum deviation angle from CSL GBs, which have relatively lower energy comparing with other non-R GBs [20]. In addition, CSL GBs have a much lower mobility as well [21]. Therefore, several simulations for annealing have been reported by taking these anisotropic energy and mobility of GBs into account [22–26]. Twodimensional (2D) grain growth during coarsening has also been simulated and discussed for the evolution of grain orientation distribution [24–26]. Nevertheless, the geometric relation and interactions between CSL GBs have not been applied to the simulation of mc-Si directional solidification. Furthermore, there were few studies on the simulation of twinning nucleation. Pohl et al. [27] found that a stable twin only existed at the three-phase boundary by using molecular dynamic simulation. Nadri et al. [28] also added twinning to their model to simulate the grain structure in mc-Si ingot by using a geometric model. The results looked close to the experiments, but the model did not consider GB types and their interactions. In this paper, we simulate the grain structure and the evolution of GBs during directional solidification of mc-Si sheet by phase field modeling. The geometric and energetic properties of R3n GBs are considered. In the next section, the phase field model used in the simulation is briefly described. Section 3 is devoted to results and discussion followed by the conclusion in Section 4. 2. Phase field modeling The phase field model (PFM) used here is based on the thininterface model proposed by Karma and Rappel [29]. The phase field variable / is set to 1 in the solid, -1 in the melt, and 0 at the interface. To represent the model in dimensionless form, the length is rescaled by W0, which characterizes the interface thickness, and the time t is rescaled by s0, which characterizes the atomic movement. The dimensionless variables are denoted by a superscript asterisk, unless otherwise stated. The dimensionless phase field equation could be written as follows:
s ðnÞ
@/ ¼ r ½W ðnÞ2 r / þ ½/ kc uð1 /2 Þð1 /2 Þ @t 1þ/ @W ðnÞ H jDhj þ @ x jr /j2 W ðnÞ 2 @ð@ x /Þ @W ðnÞ ; þ @ y jr /j2 W ðnÞ @ð@ y /Þ
where
h
i
s ðnÞ ¼ a2s ðnÞ þ a1ba02DWm 0 as ðnÞak ðnÞ
ð1Þ
for the thin-interface
model, where n is the normal unit vector at the interface and as(n) is the anisotropy function for the interfacial free energy, b0 is the kinetic coefficient, Dm is the mean thermal diffusivity, a1 and a2 are constants, and ak ðnÞ is the anisotropy function for the kinetic coefficient. Moreover, W ⁄ (n) = as(n) and kc is a coupling constant between the phase field and the temperature field. In addition, u is the dimensionless temperature, i.e., u = Cp,l (T Tm)/DH, where T is temperature, Cp,l is the specific heat of the liquid, and DH is the heat of fusion. H is a parameter related to grain boundary energy, and |Dh| is the angular difference between two grain orien-
tations, as the coupling between orientation and phase fields. The relationship between H, |Dh| and grain boundary energy cgb can be derived based on the same procedure in [30]. The highly anisotropic interfacial free energy as(n) and kinetic ak(n) functions are selected from the previous work for the facet formation of silicon [31], which showed that the wavelength of facets was affected by as(n) and ak(n) influenced the facet tips. To model polycrystalline materials, especially for describing the CSL GBs, the normal vector of each grain needs to be calculated individually and correctly. The concepts of orientation-field [32] and multi-phase-field models [33] are adopted. We introduce N crystalline variables ui to the orientation-field equation for specifying the grains, where i = 1, N. For the crystalline variables, ui Z is set to 1 for grain i, and 0 for the others. This crystalline field resemble the concept of the orientation field. The real orientation of a grain relative to a reference frame, h, is assigned to each grain when ui exceeds a threshold, e.g., 0.9 is used here. More detail derivation and procedure has been discussed in [34]. Then, the crystalline evolution equations could be derived as follows:
" # a @ ui 1þ/ r ui ¼ M ui r h PN ; 2 @t i¼1 jr ui j
i ¼ 1; 2; . . . ; N
ð2Þ
2
where Mui ¼ ð1 ð1þ/ Þ Þ aui is the mobility of each grain. More2 over, aui is the anisotropic function for grain boundary mobility; h and a are the interpolation parameters, which are set to 5 and 100 in our simulation for eliminating the crystalline field diffusion into the liquid phase. Larger a makes the boundary layer of crystalline field smaller, and this makes the determination of the orientation easier. For example, if we have four seeds, the initial profiles of the grain crystalline variables ui (i = 1, 4) can be set, as shown in Fig. 1(a). The crystalline variable ui is set to 1 for grain i, and 0 for the rest. The orientation field variable h is assigned to each grain, as shown in Fig. 1(b). In this study, we take h as the tilt angle 1i and from h0 1 1i. When h is 0 , the x-axis and y-axis are h0 1 h1 0 0i, respectively. To simplify our calculation, the frozen temperature approximation (FTA) is also adopted. The dimensionless temperature distribution is given by:
u ¼ G ðy V t Þ
ð3Þ
where G and V are the dimensionless temperature gradient and the drift velocity. A positive temperature gradient (G) of 200 K/ mm and the drift velocity (V) of 2 mm/min are used in our simulation, which are close to our previous experimental parameters [17]. The temperature profile is shown in Fig. 1(c). The above equations including one phase equation in Eq. (1) and N crystalline equations in Eq. (2) are solved by an adaptive finite volume method [35]. The interface thickness W0 is taken to be 2.5 lm, and the corresponding s0 is 9.5 103 s. The domain size and an example of the corresponding adaptive mesh are shown in Fig. 1(d). Detail numerical implementation can be found elsewhere [35]. 3. Results and discussion 3.1. GB properties Typical CSL GBs for Si h0 1 1i are R3, R9, and R27, which have relatively lower GB energy. The Read-Shockley model [18] is considered for small angle GBs. The GB energy is controlled by the parameter H and is extended from [25] as follows:
(
D h H ¼ 12:5 1 ln 15 for Dh < 15 15 ; n H ¼ 12:5 Dh ð1 A cos ðDh hCSL ÞÞ for Dh > 15
ð4Þ
152
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
Fig. 1. (a) The initial profiles of the identity variables, u1 to u4 , for the seeds; (b) the corresponding orientation variable h ¼ 70 ; 0 ; 70 ; and 70 for the four seeds from left to right, respectively, where h is the tilt angle around h0 1 1i; (c) the temperature field; (d) the adaptive mesh and domain size.
where hCSL is the misorientation for h0 1 1i CSL GBs, for example 70.5° and 109.5° for R3, 38.9° and 141.1° for R9, and 31.6° and 148.4° for R27, respectively [20]. The parameter A and n are depth and the width for the cusp function. These values are 0.95 and 128 for R3, 0.4 and 256 for R9, and 0.08 and 256 for R27, respectively. This corresponding GB energy with misorientation is shown in Fig. 2(a). As shown, relatively lower GB energy is given for the R3 GBs. Shallower cusps are also assigned for the R9 and R27 GBs. Here we assume that all the CSL GBs have the same energy. In addition to the energy, it’s also important to consider the anisotropic mobility for GBs. For CSL GBs, much lower mobility could be found if the CSL GB is a coherent twin boundary [36], so we further consider the GB plane for each CSL GB. For example, a coherent twin boundary plane for a R3 GB is {1 1 1}, and we must find the {1 1 1} twin plane of this R3 coherent twin boundary. This is much simpler in 2D simulations. The tilt angle of this {1 1 1} twin plane is an average of grain orientations of two adjacent grains. The anisotropic mobility function is given as follows [25]:
( au ¼
"
2 1
64 #)10 cosð2ðUGB Ucoherent ÞÞ þ 1 ; 2
ð5Þ
where UGB is the tilt angle of the grain boundary, which is defined as u UGB ¼ 90 U, and U ¼ sin1 jr ruj nx ; nx is the unit vector of xcoordinate. With this function, some specific tilt angles from the tilt GB direction, i.e., coherent twin boundaries, are assigned in the term Ucoherent . For example, if the grain orientations of two adjacent seeds are 0° and 70°, the R3 GB has a relatively lower mobility at Ucoherent ¼ 35 , as shown in Fig. 2(b) plotted from Eq. (5). 3.2. GB interactions To illustrate the GB interactions, three typical cases are considered: CSL with CSL, CSL with non-R, and non-R with non-R. In Fig. 3(a), the general relations between two CSL GBs can be seen from our simulation, such as R3 + R9 ? R3 and R3 + R3 ? R9. These GB types change if we change the grain orientations, as shown in Fig. 3(b). The interaction becomes R3 + non-R ? nonR. This indicates that the non-R blocked R3 GBs, which was found in the previous experiments [17]. In the third interaction, a non-R GB remains when two non-R GBs meet, as shown in Fig. 3(c). Apparently, non-R GBs are obtained most easily due to the second and third interactions, and it looks quite similar to the experiments
Fig. 2. (a) GB energy with different misorientations; (b) the specific mobility for a certain CSL coherent GB, e.g., a
P 3 coherent GB here.
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
153
P P P Fig. 3. Three typical types of GB interactions: (a) CSL with CSL; (b) CSL with non- ; (c) non- with non- . Also, the evolution of GBs in a convex thermal gradient: (d) CSL P GBs and (e) non- GBs.
Fig. 4. (a) Flow chart of nucleation algorithm in the simulation; (b) an example of twinning process.
at the lower drift velocity [17]. However, the GB interactions are more complicated in reality and a three-dimension simulation is necessary to describe the GB more accurately.
To consider the effect of thermal gradients, a temperature profile for a convex growth front is considered. The results for the convex FTA are shown in Fig. 3(d) and (e). As shown, the CSL GBs are
154
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
not affected by the temperature gradients, and they grow with their own direction. This is similar to the experiments [17], because the CSL GBs have a lower mobility. On the other hand, the non-R GBs follow the direction of thermal gradients. This also resembles the experimental observations that most of the non-R GBs tend to align with the direction of thermal gradients. 3.3. Grain growth including twinning Twin nucleation occur rather frequently during mc-Si crystal growth, especially at a higher growth speed [17]. The twinning mechanism for grain nucleation in mc-Si was first proposed by
Duffar and Nadri [8]. It was modified recently by Lin and Lan [9] by considering the neighbor grain and the GB type. Here we incorporate this modified model in our simulation for twinning. The flow chart of the algorithm is illustrated in Fig. 4(a). As shown, after each time step snuc, we first find all tri-junctions (solid-solidliquid) from the intersection of the solid-liquid interface and GBs. These tri-junctions can be viewed as the nucleation sites. Then we calculate their twinning probabilities PT based on the formula derived by Lin and Lan [9]. In each nucleation site, there exists two twinning probabilities because the tri-junction consists of two facet planes. We then calculate these twinning probabilities and pick the one with the highest PT value, which is presumed as
Fig. 5. (a) Grain structure with different initial distributions of GB types; (b) the corresponding GB types; the he evolution of the number of twin nucleation (c), average grain size (d) and GB types (e).
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
the criterion for twin nucleation in this nucleation site. A random number P between 0 to 1 is generated by computer and a nucleus is introduced in this nucleation site if P < PT . An example of twinning process is given in Fig. 4(b). As shown, a new orange grain appear at a certain time, and a tilted R3 GB appears after the nucleation succeeds. By considering a more realistic simulation, we put 20 seeds initially and simulate the evolution of the grain orientation, type of GBs, grain size, and the twin nucleation number for some parameters. First, we study the effect of initial GB types. Here we simulate two cases; the first case has seeds with 70% R3 GBs and the second has seeds with 70% non-R GBs initially. As shown in Fig. 5(a), we can see the evolution of the grain structure. The black dash line
155
is the initial interface position, and the grains grow upward. For the case of seeds with 70% R3 GBs, most twin nucleation occur at the left region due to the initial random orientation distribution. Grains at the central region grow straightly along the thermal gradient. In contrast, for the case of seeds with 70% non-R GBs, twin nucleation appears everywhere. The corresponding GB types for these two cases are shown in Fig. 5(b). Tilted R3 GBs generate after twin nucleation, but they are usually blocked by non-R GBs. Moreover, the evolution of the twin nucleation number is shown in Fig. 5(c). Higher frequency of the occurrence of twin nucleus can be observed for the case of seeds with 70% non-R GBs. Nevertheless, the nucleation rate gradually decreases as grains grow in both cases due to the reduction of nucleation sites
Fig. 6. (a) Grain structure with different initial distributions of grain orientations; (b) the corresponding GB types. The evolution of the number of twin nucleation (c), average grain size (d), and GB types (e).
156
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
while grains coarsen. Moreover, the evolution of average grain size is shown in Fig. 5(d); the average grain size is estimated by the domain width divided by the number of grains. Smaller average grain size for the case of seeds with 70% non-R GBs can be seen due to more nucleated grains during twinning. In addition, we could also trace the variation of GB types. As shown in Fig. 5(e), for the case of seeds with 70% R3 GBs, the proportion of R3 GBs decreases with the increasing R9 GBs at early stage of solidification owing to R3 + R3 ? R9. As grains grow further, the proportion of R3 GBs gradually decreases with the increasing non-R GBs. The reason is that R3 GBs are blocked by non-R GBs; therefore, the proportion of non-R GBs is higher than that of R3 GBs at the final stage. On the other hand, for the case of seeds with 70% non-R GBs, because of the twin nucleation, the proportion of non-R GBs gradually decreases with the increasing R3 GBs. However, the proportion of non-R GBs is still higher than R3 GBs at the final stage. In general, near the end the ratio of non-R GBs to R3 GBs becomes similar, which is similar to the experiments by Wong et al. [12]. This indicates that the frequencies of appearance and disappearance of R3 GBs are comparable. Secondly, we discuss the effect of the initial distribution of grain orientation. When the initial orientation h is concentrated between 20° and 20°, we refer it to the case of 70% {1 0 0} seeds initially. In contrast, the case of 70% {1 1 1} seeds initially can be represented by the initial orientation h concentrated between 65° to 45° and 45° to 65°, respectively. As shown in Fig. 6(a) for the evolution of grain structure, more small grains emerge during the growth from 70% {1 0 0} seeds. On the other hand, most grains grow straightly and not so many grains nucleate in the case of 70% {1 1 1} seeds. The corresponding evolution of GB types is also shown in Fig. 6(b). Obviously, there are many tilted R3 GBs in the case of 70% {1 0 0} seeds. For the case of 70% {1 1 1} seeds, most
non-R GBs grow upward due to the large thermal gradient and less R3 GBs appear during solidification. The twin nucleation number is shown in Fig. 6(c). The twining frequency is higher in the case of 70% {1 0 0} seeds. However, as shown in Fig. 6(d), the difference in the average grain size evolution is small for both cases. This is because the R3 GBs could be easily blocked by non-R GBs. In other words, the competition between the forming and blocking for R3 GBs results in the comparable grain size for both cases. Again, the proportion of GB types is also analyzed in Fig. 6(e). In both cases, we set over 80% non-R GBs initially, and they all decrease during solidification. The proportion of R3 GBs increases more in the case of 70% {1 0 0} seeds due to more twin nucleation, but the blocking of R3 GBs from non-R GBs inhibits the increase of the proportion of R3 GBs. Less twinning for the case of 70% {1 1 1} seeds and thus lower proportion of R3 GBs can be seen clearly. To summarize the simulation results of the present section, the effect of GB types on the twinning can be resulted from the undercooling at the GB groove. Larger undercooling induces more twin nucleation, and the undercooling at the GB grooves with different GB types is extracted from our simulation, as shown in Fig. 7(a). There is a large difference in the undercooling between non-R and R3 GBs. Deeper grooves for non-R GBs lead to a larger undercooling than that for R3 GBs, which form small or invisible grooves from the experiments of Fujiwara et al. [37]. Besides, we can also see more twin nucleation occur in the case of 70% {1 0 0} seeds. As shown in Fig. 7(b), the twinning probability changes with the angle m between the facet and the GB plane [9]. A schematic diagram for the facet-facet groove is also depicted. As the GB groove consists of two {1 1 1} grains, twinning probability is rather low. This is the reason why much less twin nucleation is observed in the solidification from 70% {1 1 1} seeds than that from 70% {1 0 0} seeds. As a result, more {1 1 1} oriented seeds seem to be favored by the HP-mc Si because fewer R3 GBs are generated from twin nucleation according to our simulation. On the other hand, if the growth speed is reduced, the undercooling in the tri-junction groove is reduced, and this reduces twinning probability and the generation of R3 GBs.
4. Conclusions
Fig. 7. (a) Undercooling at the grooves of different GB types from simulation; the groove shapes of non-R and R3 GBs are shown. (b) The relation between twinning probability and the facet angle m. The schematic illustrations of the facet-facet groove for {1 1 1} and {1 0 0} grains are shown.
In this paper, we use the phase field and crystalline field models to simulate the development of grain structures and GBs during the solidification of mc-Si sheets. The CSL GBs are described successfully by the anisotropic grain boundary energy and mobility. Typical interactions between CSL and non-R GBs are examined and they are consistent with the experiments. The temperature gradient direction only affects the growth direction for non-R GBs, but not for CSL GBs. By applying the twinning model, the effects of initial distribution of GB types and grain orientations are also studied. Obviously, the frequency of the occurrence of twin nucleation is higher for non-R GBs due to the larger undercooling at their tri-junction grooves. The twinning probability is also higher for the solidification from {1 0 0} oriented seeds owing to the larger angle between the facet and the GB plane. Moreover, the tilted R3 GBs generated from twinning are supposed to decrease the average grain size during solidification. However, R3 GBs are easily blocked by the non-R GBs. Furthermore, we find that near the end of solidification, the proportion of GB types is independent of its initial state. The appearance of R3 GBs from twinning becomes comparable to the disappearance of R3 GBs due to the blocking by non-R GBs after some time. Hence, the ratio of non-R GBs to R3 GBs remains about the same till the end of solidification. Of course, to better simulate the experimental results in a realistic manner, 3D simulation would be necessary, but this is beyond the scope of this study.
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
Acknowledgement This work was supported by the Ministry of Science and Technology of Taiwan.
References [1] K. Fujiwara, W. Pan, K. Sawada, M. Tokairin, N. Usami, Y. Nose, A. Nomura, T. Shishido, K. Nakajima, Directional growth method to obtain high quality polycrystalline silicon from its melt, J. Cryst. Growth 292 (2) (2006) 282–285. [2] L. Stockmeier, G. Müller, A. Seidl, T. Lehmann, C. Reimann, J. Friedrich, Preferred grain orientations in silicon ribbons grown by the string ribbon and the edge-defined film-fed growth methods, J. Cryst. Growth 395 (2014) 74–79. [3] Y.M. Yang, A. Yu, B. Hsu, W.C. Hsu, A. Yang, C.W. Lan, Development of highperformance multicrystalline silicon for photovoltaic industry, Prog. Photovolt. Res. Appl. 23 (3) (2015) 340–351. [4] G. Stokkan, Y. Hu, Ø. Mjøs, M. Juel, Study of evolution of dislocation clusters in high performance multicrystalline silicon, Sol. Energ. Mater. Sol. C. 130 (2014) 679–685. [5] C.C. Hsieh, Y.C. Wu, A. Lan, H.P. Hsu, C. Hsu, C.W. Lan, Comparison of defect formations in solar silicon growth from small random and large oriented seeds, J. Cryst. Growth 419 (2015) 1–6. [6] K. Fujiwara, Y. Obinata, T. Ujihara, N. Usami, G. Sazaki, K. Nakajima, Grain growth behaviors of polycrystalline silicon during melt growth processes, J. Cryst. Growth 266 (4) (2004) 441–448. [7] P. Chen, Y.L. Tsai, C.W. Lan, Phase field modeling of growth competition of silicon grains, Acta Mater. 56 (15) (2008) 4114–4122. [8] T. Duffar, A. Nadri, On the twinning occurrence in bulk semiconductor crystal growth, Scripta Mater. 62 (12) (2010) 955–960. [9] H.K. Lin, C.W. Lan, Revisiting the twinning mechanism in directional solidification of multi-cyrstalline silicon sheet, Acta Mater. 131 (2017) 1–10. [10] A. Tandjaoui, N. Mangelinck-Noel, G. Reinhart, B. Billia, X. Guichard, Twinning occurrence and grain competition in multi-crystalline silicon during solidification, C. R. Phys. 14 (2–3) (2013) 141–148. [11] T. Riberi-Béridot, N. Mangelinck-Noël, A. Tandjaoui, G. Reinhart, B. Billia, T. Lafford, J. Baruchel, L. Barrallier, On the impact of twinning on the formation of the grain structure of multi-crystalline silicon for photovoltaic applications during directional solidification, J. Cryst. Growth 418 (2015) 38–44. [12] Y.T. Wong, C. Hsu, C.W. Lan, Development of grain structures of multicrystalline silicon from randomly orientated seeds in directional solidification, J. Cryst. Growth 387 (2014) 10–15. [13] V. Randle, Twinning-related grain boundary engineering, Acta Mater. 52 (14) (2004) 4067–4081. [14] B.W. Reed, M. Kumar, Mathematical methods for analyzing highly-twinned grain boundary networks, Scripta Mater. 54 (6) (2006) 1029–1033. [15] V. Randle, Y. Hu, The role of vicinal Sigma 3 boundaries and Sigma 9 boundaries in grain boundary engineering, J. Mater. Sci. 40 (12) (2005) 3243– 3246. [16] M. Trempa, I. Kupka, C. Kranert, T. Lehmann, C. Reimann, J. Friedrich, Evolution of grain structure and recombination active dislocations in extraordinary tall conventional and high performance multi-crystalline silicon ingots, J. Cryst. Growth 459 (2017) 67–75.
157
[17] H.K. Lin, M.C. Wu, C.C. Chen, C.W. Lan, Evolution of grain structures during directional solidification of silicon wafers, J. Cryst. Growth 439 (2016) 40–46. [18] W.T. Read, W. Shockley, Dislocation Models of Crystal Grain Boundaries, Phys. Rev. 78 (3) (1950) 275–289. [19] D.G. Brandon, The structure of high-angle grain boundaries, Acta Metall. 14 (11) (1966) 1479–1484. [20] M. Kohyama, R. Yamamoto, M. Doyama, Structures and energies of symmetrical <011> tilt grain boundaries in silicon, Phys. Status Solidi B 137 (1) (1986) 11–20. [21] K.G.F. Janssens, D. Olmsted, E.A. Holm, S.M. Foiles, S.J. Plimpton, P.M. Derlet, Computing the mobility of grain boundaries, Nat. Mater. 5 (2) (2006) 124–127. [22] A.Y. Badmos, H.J. Frost, I. Baker, Simulation of microstructural evolution during directional annealing with variable boundary energy and mobility, Acta Mater. 51 (10) (2003) 2755–2764. [23] H. Håkan, Influence of anisotropic grain boundary properties on the evolution of grain boundary character distribution during grain growth—a 2D level set study, Model. Simul. Mater. Sci. Eng. 22 (8) (2014) 085005. [24] N. Moelans, B. Blanpain, P. Wollants, Quantitative phase-field approach for simulating grain growth in anisotropic systems with arbitrary inclination and misorientation dependence, Phys. Rev. Lett. 101 (2) (2008) 025502. [25] T. Hirouchi, T. Tsuru, Y. Shibutani, Grain growth prediction with inclination dependence of h110i tilt grain boundary using multi-phase-field model with penalty for multiple junctions, Comput. Mater. Sci. 53 (1) (2012) 474–482. [26] T. Hirouchi, Y. Shibutani, Explicit distinctions between 2D MPF grain growth simulations and EBSD analyses to determine driving mechanism of grain growth, Mater. Trans. 54 (10) (2013) 1884–1893. [27] J. Pohl, M. Müller, A. Seidl, K. Albe, Formation of parallel (1 1 1) twin boundaries in silicon growth from the melt explained by molecular dynamics simulations, J. Cryst. Growth 312 (8) (2010) 1411–1415. [28] A. Nadri, Y. Duterrail-Couvat, T. Duffar, Two-dimensional numerical modeling of grain structure in multi-crystalline silicon ingot, J. Cryst. Growth 385 (2014) 16–21. [29] A. Karma, W.-J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. 57 (4) (1998) 4323–4349. [30] J.A. Warren, R. Kobayashi, A.E. Lobovsky, W.C. Carter, Extending phase field models of solidification to polycrystalline materials, Acta Mater. 51 (20) (2003) 6035–6058. [31] H.K. Lin, H.Y. Chen, C.W. Lan, Phase field modeling of facet formation during directional solidification of silicon film, J. Cryst. Growth 385 (2014) 134–139. [32] R. Kobayashi, J.A. Warren, W.C. Carter, Vector-valued phase field model for crystallization and grain boundary formation, Physica D 119 (1998) 415–423. [33] I. Steinbach, F. Pezzolla, B. Nestler, M. Seesselberg, R. Prieler, G.J. Schmitz, J.L.L. Rezende, A phase field concept for multiphase systems, Physica D 94 (1996) 135–147. [34] G.Y. Chen, H.K. Lin, C.W. Lan, Phase-field modeling of twin-related faceted dendrite growth of silicon, Acta Mater. 115 (2016) 324–332. [35] C.W. Lan, C.C. Liu, C.M. Hsu, An adaptive finite volume method for incompressible heat flow problems in solidification, J. Comput. Phys. 178 (2) (2002) 464–497. [36] D.L. Olmsted, E.A. Holm, S.M. Foiles, Survey of computed grain boundary properties in face-centered cubic metals-II: grain boundary mobility, Acta Mater. 57 (13) (2009) 3704–3713. [37] K. Fujiwara, M. Ishii, K. Maeda, H. Koizumi, J. Nozawa, S. Uda, The effect of grain boundary characteristics on the morphology of the crystal/melt interface of multicrystalline silicon, Scripta Mater. 69 (3) (2013) 266–269.
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/crys
Phase field modeling of grain structure evolution during directional solidification of multi-crystalline silicon sheet H.K. Lin, C.W. Lan ⇑ Department of Chemical Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 7 March 2017 Received in revised form 26 May 2017 Accepted 16 June 2017 Available online 19 June 2017 Communicated by M. Plapp Keywords: A1. Computer simulation A1. Directional solidification A1. Planar defects A1. Nucleation B2. Semiconducting silicon
a b s t r a c t Evolution of grain structures and grain boundaries (GBs), especially the coincident site lattice GBs, during directional solidification of multi-crystalline silicon sheet are simulated by using a phase field model for the first time. Since the coincident site lattice GBs having lower mobility, tend to follow their own crystallographic directions despite thermal gradients, the anisotropic energy and mobility of GBs are considered in the model. Three basic interactions of GBs during solidification are examined and they are consistent with experiments. The twinning process for new grain formation is further added in the simulation by considering twin nucleation. The effect of initial distribution of GB types and grain orientations is also investigated for the twinning frequency and the evolution of grain size and GB types. Ó 2017 Published by Elsevier B.V.
1. Introduction Grain boundaries (GBs) of multi-crystalline silicon (mc-Si), grown by directional solidification, play a crucial role on the electrical properties of solar cells. Controlling GBs is often necessary for improving the ingot quality. For example, the dendritic casting method [1] is to initiate many highly symmetric coincident site lattice (CSL) GBs, e.g., R3 GBs, through grain growth at high undercooling generated along the crucible bottom [1], and they are useful for the high lifetime in the grown ingot. In ribbon growth, the twin grains have better lifetime and their formation have also been discussed by Stockmeier et al. [2] based on thermodynamics arguments. Nevertheless, the detailed GB interactions and how to control them are still not quite clear. Recently, the ingot growth has made a significant progress by initiating uniform and small grains at the crucible bottom, i.e., the so-called high-performance mc-Si (HPMC) [3]. The large amount of non-R GBs, generated at the beginning of solidification, would relax the thermal stress and thus reduce the generation of the dislocations; they also terminate the propagation of the dislocations during ingot growth [3–5]. However, during ingot growth, more R3 GBs appear while the proportion of non-R GBs decreases, and the ingot quality near the top is thus deteriorated. Therefore, understanding the evolution of grain structures is very important to further control the directional solidification of Si ingots or ribbons. ⇑ Corresponding author. E-mail address: [email protected] (C.W. Lan). http://dx.doi.org/10.1016/j.jcrysgro.2017.06.016 0022-0248/Ó 2017 Published by Elsevier B.V.
GB evolution during directional solidification of silicon has been paid much attention in recent years. Fujiwara et al. [6] found that the low-energy grain, such as h1 1 1i, tends to overgrow the others at low undercooling. Chen et al. [7] were the first to simulate this phenomenon by using a two-dimensional phase field model considering anisotropic interfacial energy and kinetic coefficient with some success. Duffar and Nadri [8] further described the twinning mechanism for the nucleation on the {1 1 1} facet plane near GBs. However, the critical undercooling proposed to form a twinned nucleus in their mechanism was 9 K or larger. Their model was recently modified by Lin and Lan [9], who considered the interactions between the nucleus and the neighbor grain at the groove. A more reasonable undercooling for the twin nucleation on the {1 1 1} facet plane was derived, i.e., less than 1 K. This value was more consistent with the experimental observation [10]. To confirm the twin formation, in-situ X-ray synchrotron imaging was used to investigate silicon crystal growth by Tandjaoui et al. [10]. They characterized the birth of new grains and showed that twins nucleated exactly on {1 1 1} facets at the GB grooves [11]. Moreover, Wong et al. [12] investigated the evolution of grain orientations and GBs by analyzing wafers at different heights from directional solidification ingots seeded with small and randomly oriented silicon beads. Indeed, twin nucleation was found at the tri-junctions. In addition to these research, there have been several works related to the interactions between R3n GBs [13–16]. Recently, the interactions between R3n and non-R GBs were further investigated by Lin et al. [17] for the directional solidification of mc-Si sheets at different speeds. At lower speeds, they found
151
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
that R3 GBs decreased due to the interactions between R3 GBs and non-R GBs, while twinning was a key mechanism for the increase of R3 GBs at higher speeds. Their study using mc-Si sheets unfolded the detailed mechanisms of GB interactions in silicon. However, up till now these have not yet been simulated successfully. Moreover, the anisotropic properties of GBs in silicon play a crucial role in the GB interactions, and thus in their evolution. The classic Read and Shockley dislocation model [18] describes the energy in low-angle non-R GBs. For high-angle GBs, the wellknown Brandon criteria [19] are used to calculate the maximum deviation angle from CSL GBs, which have relatively lower energy comparing with other non-R GBs [20]. In addition, CSL GBs have a much lower mobility as well [21]. Therefore, several simulations for annealing have been reported by taking these anisotropic energy and mobility of GBs into account [22–26]. Twodimensional (2D) grain growth during coarsening has also been simulated and discussed for the evolution of grain orientation distribution [24–26]. Nevertheless, the geometric relation and interactions between CSL GBs have not been applied to the simulation of mc-Si directional solidification. Furthermore, there were few studies on the simulation of twinning nucleation. Pohl et al. [27] found that a stable twin only existed at the three-phase boundary by using molecular dynamic simulation. Nadri et al. [28] also added twinning to their model to simulate the grain structure in mc-Si ingot by using a geometric model. The results looked close to the experiments, but the model did not consider GB types and their interactions. In this paper, we simulate the grain structure and the evolution of GBs during directional solidification of mc-Si sheet by phase field modeling. The geometric and energetic properties of R3n GBs are considered. In the next section, the phase field model used in the simulation is briefly described. Section 3 is devoted to results and discussion followed by the conclusion in Section 4. 2. Phase field modeling The phase field model (PFM) used here is based on the thininterface model proposed by Karma and Rappel [29]. The phase field variable / is set to 1 in the solid, -1 in the melt, and 0 at the interface. To represent the model in dimensionless form, the length is rescaled by W0, which characterizes the interface thickness, and the time t is rescaled by s0, which characterizes the atomic movement. The dimensionless variables are denoted by a superscript asterisk, unless otherwise stated. The dimensionless phase field equation could be written as follows:
s ðnÞ
@/ ¼ r ½W ðnÞ2 r / þ ½/ kc uð1 /2 Þð1 /2 Þ @t 1þ/ @W ðnÞ H jDhj þ @ x jr /j2 W ðnÞ 2 @ð@ x /Þ @W ðnÞ ; þ @ y jr /j2 W ðnÞ @ð@ y /Þ
where
h
i
s ðnÞ ¼ a2s ðnÞ þ a1ba02DWm 0 as ðnÞak ðnÞ
ð1Þ
for the thin-interface
model, where n is the normal unit vector at the interface and as(n) is the anisotropy function for the interfacial free energy, b0 is the kinetic coefficient, Dm is the mean thermal diffusivity, a1 and a2 are constants, and ak ðnÞ is the anisotropy function for the kinetic coefficient. Moreover, W ⁄ (n) = as(n) and kc is a coupling constant between the phase field and the temperature field. In addition, u is the dimensionless temperature, i.e., u = Cp,l (T Tm)/DH, where T is temperature, Cp,l is the specific heat of the liquid, and DH is the heat of fusion. H is a parameter related to grain boundary energy, and |Dh| is the angular difference between two grain orien-
tations, as the coupling between orientation and phase fields. The relationship between H, |Dh| and grain boundary energy cgb can be derived based on the same procedure in [30]. The highly anisotropic interfacial free energy as(n) and kinetic ak(n) functions are selected from the previous work for the facet formation of silicon [31], which showed that the wavelength of facets was affected by as(n) and ak(n) influenced the facet tips. To model polycrystalline materials, especially for describing the CSL GBs, the normal vector of each grain needs to be calculated individually and correctly. The concepts of orientation-field [32] and multi-phase-field models [33] are adopted. We introduce N crystalline variables ui to the orientation-field equation for specifying the grains, where i = 1, N. For the crystalline variables, ui Z is set to 1 for grain i, and 0 for the others. This crystalline field resemble the concept of the orientation field. The real orientation of a grain relative to a reference frame, h, is assigned to each grain when ui exceeds a threshold, e.g., 0.9 is used here. More detail derivation and procedure has been discussed in [34]. Then, the crystalline evolution equations could be derived as follows:
" # a @ ui 1þ/ r ui ¼ M ui r h PN ; 2 @t i¼1 jr ui j
i ¼ 1; 2; . . . ; N
ð2Þ
2
where Mui ¼ ð1 ð1þ/ Þ Þ aui is the mobility of each grain. More2 over, aui is the anisotropic function for grain boundary mobility; h and a are the interpolation parameters, which are set to 5 and 100 in our simulation for eliminating the crystalline field diffusion into the liquid phase. Larger a makes the boundary layer of crystalline field smaller, and this makes the determination of the orientation easier. For example, if we have four seeds, the initial profiles of the grain crystalline variables ui (i = 1, 4) can be set, as shown in Fig. 1(a). The crystalline variable ui is set to 1 for grain i, and 0 for the rest. The orientation field variable h is assigned to each grain, as shown in Fig. 1(b). In this study, we take h as the tilt angle 1i and from h0 1 1i. When h is 0 , the x-axis and y-axis are h0 1 h1 0 0i, respectively. To simplify our calculation, the frozen temperature approximation (FTA) is also adopted. The dimensionless temperature distribution is given by:
u ¼ G ðy V t Þ
ð3Þ
where G and V are the dimensionless temperature gradient and the drift velocity. A positive temperature gradient (G) of 200 K/ mm and the drift velocity (V) of 2 mm/min are used in our simulation, which are close to our previous experimental parameters [17]. The temperature profile is shown in Fig. 1(c). The above equations including one phase equation in Eq. (1) and N crystalline equations in Eq. (2) are solved by an adaptive finite volume method [35]. The interface thickness W0 is taken to be 2.5 lm, and the corresponding s0 is 9.5 103 s. The domain size and an example of the corresponding adaptive mesh are shown in Fig. 1(d). Detail numerical implementation can be found elsewhere [35]. 3. Results and discussion 3.1. GB properties Typical CSL GBs for Si h0 1 1i are R3, R9, and R27, which have relatively lower GB energy. The Read-Shockley model [18] is considered for small angle GBs. The GB energy is controlled by the parameter H and is extended from [25] as follows:
(
D h H ¼ 12:5 1 ln 15 for Dh < 15 15 ; n H ¼ 12:5 Dh ð1 A cos ðDh hCSL ÞÞ for Dh > 15
ð4Þ
152
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
Fig. 1. (a) The initial profiles of the identity variables, u1 to u4 , for the seeds; (b) the corresponding orientation variable h ¼ 70 ; 0 ; 70 ; and 70 for the four seeds from left to right, respectively, where h is the tilt angle around h0 1 1i; (c) the temperature field; (d) the adaptive mesh and domain size.
where hCSL is the misorientation for h0 1 1i CSL GBs, for example 70.5° and 109.5° for R3, 38.9° and 141.1° for R9, and 31.6° and 148.4° for R27, respectively [20]. The parameter A and n are depth and the width for the cusp function. These values are 0.95 and 128 for R3, 0.4 and 256 for R9, and 0.08 and 256 for R27, respectively. This corresponding GB energy with misorientation is shown in Fig. 2(a). As shown, relatively lower GB energy is given for the R3 GBs. Shallower cusps are also assigned for the R9 and R27 GBs. Here we assume that all the CSL GBs have the same energy. In addition to the energy, it’s also important to consider the anisotropic mobility for GBs. For CSL GBs, much lower mobility could be found if the CSL GB is a coherent twin boundary [36], so we further consider the GB plane for each CSL GB. For example, a coherent twin boundary plane for a R3 GB is {1 1 1}, and we must find the {1 1 1} twin plane of this R3 coherent twin boundary. This is much simpler in 2D simulations. The tilt angle of this {1 1 1} twin plane is an average of grain orientations of two adjacent grains. The anisotropic mobility function is given as follows [25]:
( au ¼
"
2 1
64 #)10 cosð2ðUGB Ucoherent ÞÞ þ 1 ; 2
ð5Þ
where UGB is the tilt angle of the grain boundary, which is defined as u UGB ¼ 90 U, and U ¼ sin1 jr ruj nx ; nx is the unit vector of xcoordinate. With this function, some specific tilt angles from the tilt GB direction, i.e., coherent twin boundaries, are assigned in the term Ucoherent . For example, if the grain orientations of two adjacent seeds are 0° and 70°, the R3 GB has a relatively lower mobility at Ucoherent ¼ 35 , as shown in Fig. 2(b) plotted from Eq. (5). 3.2. GB interactions To illustrate the GB interactions, three typical cases are considered: CSL with CSL, CSL with non-R, and non-R with non-R. In Fig. 3(a), the general relations between two CSL GBs can be seen from our simulation, such as R3 + R9 ? R3 and R3 + R3 ? R9. These GB types change if we change the grain orientations, as shown in Fig. 3(b). The interaction becomes R3 + non-R ? nonR. This indicates that the non-R blocked R3 GBs, which was found in the previous experiments [17]. In the third interaction, a non-R GB remains when two non-R GBs meet, as shown in Fig. 3(c). Apparently, non-R GBs are obtained most easily due to the second and third interactions, and it looks quite similar to the experiments
Fig. 2. (a) GB energy with different misorientations; (b) the specific mobility for a certain CSL coherent GB, e.g., a
P 3 coherent GB here.
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
153
P P P Fig. 3. Three typical types of GB interactions: (a) CSL with CSL; (b) CSL with non- ; (c) non- with non- . Also, the evolution of GBs in a convex thermal gradient: (d) CSL P GBs and (e) non- GBs.
Fig. 4. (a) Flow chart of nucleation algorithm in the simulation; (b) an example of twinning process.
at the lower drift velocity [17]. However, the GB interactions are more complicated in reality and a three-dimension simulation is necessary to describe the GB more accurately.
To consider the effect of thermal gradients, a temperature profile for a convex growth front is considered. The results for the convex FTA are shown in Fig. 3(d) and (e). As shown, the CSL GBs are
154
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
not affected by the temperature gradients, and they grow with their own direction. This is similar to the experiments [17], because the CSL GBs have a lower mobility. On the other hand, the non-R GBs follow the direction of thermal gradients. This also resembles the experimental observations that most of the non-R GBs tend to align with the direction of thermal gradients. 3.3. Grain growth including twinning Twin nucleation occur rather frequently during mc-Si crystal growth, especially at a higher growth speed [17]. The twinning mechanism for grain nucleation in mc-Si was first proposed by
Duffar and Nadri [8]. It was modified recently by Lin and Lan [9] by considering the neighbor grain and the GB type. Here we incorporate this modified model in our simulation for twinning. The flow chart of the algorithm is illustrated in Fig. 4(a). As shown, after each time step snuc, we first find all tri-junctions (solid-solidliquid) from the intersection of the solid-liquid interface and GBs. These tri-junctions can be viewed as the nucleation sites. Then we calculate their twinning probabilities PT based on the formula derived by Lin and Lan [9]. In each nucleation site, there exists two twinning probabilities because the tri-junction consists of two facet planes. We then calculate these twinning probabilities and pick the one with the highest PT value, which is presumed as
Fig. 5. (a) Grain structure with different initial distributions of GB types; (b) the corresponding GB types; the he evolution of the number of twin nucleation (c), average grain size (d) and GB types (e).
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
the criterion for twin nucleation in this nucleation site. A random number P between 0 to 1 is generated by computer and a nucleus is introduced in this nucleation site if P < PT . An example of twinning process is given in Fig. 4(b). As shown, a new orange grain appear at a certain time, and a tilted R3 GB appears after the nucleation succeeds. By considering a more realistic simulation, we put 20 seeds initially and simulate the evolution of the grain orientation, type of GBs, grain size, and the twin nucleation number for some parameters. First, we study the effect of initial GB types. Here we simulate two cases; the first case has seeds with 70% R3 GBs and the second has seeds with 70% non-R GBs initially. As shown in Fig. 5(a), we can see the evolution of the grain structure. The black dash line
155
is the initial interface position, and the grains grow upward. For the case of seeds with 70% R3 GBs, most twin nucleation occur at the left region due to the initial random orientation distribution. Grains at the central region grow straightly along the thermal gradient. In contrast, for the case of seeds with 70% non-R GBs, twin nucleation appears everywhere. The corresponding GB types for these two cases are shown in Fig. 5(b). Tilted R3 GBs generate after twin nucleation, but they are usually blocked by non-R GBs. Moreover, the evolution of the twin nucleation number is shown in Fig. 5(c). Higher frequency of the occurrence of twin nucleus can be observed for the case of seeds with 70% non-R GBs. Nevertheless, the nucleation rate gradually decreases as grains grow in both cases due to the reduction of nucleation sites
Fig. 6. (a) Grain structure with different initial distributions of grain orientations; (b) the corresponding GB types. The evolution of the number of twin nucleation (c), average grain size (d), and GB types (e).
156
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
while grains coarsen. Moreover, the evolution of average grain size is shown in Fig. 5(d); the average grain size is estimated by the domain width divided by the number of grains. Smaller average grain size for the case of seeds with 70% non-R GBs can be seen due to more nucleated grains during twinning. In addition, we could also trace the variation of GB types. As shown in Fig. 5(e), for the case of seeds with 70% R3 GBs, the proportion of R3 GBs decreases with the increasing R9 GBs at early stage of solidification owing to R3 + R3 ? R9. As grains grow further, the proportion of R3 GBs gradually decreases with the increasing non-R GBs. The reason is that R3 GBs are blocked by non-R GBs; therefore, the proportion of non-R GBs is higher than that of R3 GBs at the final stage. On the other hand, for the case of seeds with 70% non-R GBs, because of the twin nucleation, the proportion of non-R GBs gradually decreases with the increasing R3 GBs. However, the proportion of non-R GBs is still higher than R3 GBs at the final stage. In general, near the end the ratio of non-R GBs to R3 GBs becomes similar, which is similar to the experiments by Wong et al. [12]. This indicates that the frequencies of appearance and disappearance of R3 GBs are comparable. Secondly, we discuss the effect of the initial distribution of grain orientation. When the initial orientation h is concentrated between 20° and 20°, we refer it to the case of 70% {1 0 0} seeds initially. In contrast, the case of 70% {1 1 1} seeds initially can be represented by the initial orientation h concentrated between 65° to 45° and 45° to 65°, respectively. As shown in Fig. 6(a) for the evolution of grain structure, more small grains emerge during the growth from 70% {1 0 0} seeds. On the other hand, most grains grow straightly and not so many grains nucleate in the case of 70% {1 1 1} seeds. The corresponding evolution of GB types is also shown in Fig. 6(b). Obviously, there are many tilted R3 GBs in the case of 70% {1 0 0} seeds. For the case of 70% {1 1 1} seeds, most
non-R GBs grow upward due to the large thermal gradient and less R3 GBs appear during solidification. The twin nucleation number is shown in Fig. 6(c). The twining frequency is higher in the case of 70% {1 0 0} seeds. However, as shown in Fig. 6(d), the difference in the average grain size evolution is small for both cases. This is because the R3 GBs could be easily blocked by non-R GBs. In other words, the competition between the forming and blocking for R3 GBs results in the comparable grain size for both cases. Again, the proportion of GB types is also analyzed in Fig. 6(e). In both cases, we set over 80% non-R GBs initially, and they all decrease during solidification. The proportion of R3 GBs increases more in the case of 70% {1 0 0} seeds due to more twin nucleation, but the blocking of R3 GBs from non-R GBs inhibits the increase of the proportion of R3 GBs. Less twinning for the case of 70% {1 1 1} seeds and thus lower proportion of R3 GBs can be seen clearly. To summarize the simulation results of the present section, the effect of GB types on the twinning can be resulted from the undercooling at the GB groove. Larger undercooling induces more twin nucleation, and the undercooling at the GB grooves with different GB types is extracted from our simulation, as shown in Fig. 7(a). There is a large difference in the undercooling between non-R and R3 GBs. Deeper grooves for non-R GBs lead to a larger undercooling than that for R3 GBs, which form small or invisible grooves from the experiments of Fujiwara et al. [37]. Besides, we can also see more twin nucleation occur in the case of 70% {1 0 0} seeds. As shown in Fig. 7(b), the twinning probability changes with the angle m between the facet and the GB plane [9]. A schematic diagram for the facet-facet groove is also depicted. As the GB groove consists of two {1 1 1} grains, twinning probability is rather low. This is the reason why much less twin nucleation is observed in the solidification from 70% {1 1 1} seeds than that from 70% {1 0 0} seeds. As a result, more {1 1 1} oriented seeds seem to be favored by the HP-mc Si because fewer R3 GBs are generated from twin nucleation according to our simulation. On the other hand, if the growth speed is reduced, the undercooling in the tri-junction groove is reduced, and this reduces twinning probability and the generation of R3 GBs.
4. Conclusions
Fig. 7. (a) Undercooling at the grooves of different GB types from simulation; the groove shapes of non-R and R3 GBs are shown. (b) The relation between twinning probability and the facet angle m. The schematic illustrations of the facet-facet groove for {1 1 1} and {1 0 0} grains are shown.
In this paper, we use the phase field and crystalline field models to simulate the development of grain structures and GBs during the solidification of mc-Si sheets. The CSL GBs are described successfully by the anisotropic grain boundary energy and mobility. Typical interactions between CSL and non-R GBs are examined and they are consistent with the experiments. The temperature gradient direction only affects the growth direction for non-R GBs, but not for CSL GBs. By applying the twinning model, the effects of initial distribution of GB types and grain orientations are also studied. Obviously, the frequency of the occurrence of twin nucleation is higher for non-R GBs due to the larger undercooling at their tri-junction grooves. The twinning probability is also higher for the solidification from {1 0 0} oriented seeds owing to the larger angle between the facet and the GB plane. Moreover, the tilted R3 GBs generated from twinning are supposed to decrease the average grain size during solidification. However, R3 GBs are easily blocked by the non-R GBs. Furthermore, we find that near the end of solidification, the proportion of GB types is independent of its initial state. The appearance of R3 GBs from twinning becomes comparable to the disappearance of R3 GBs due to the blocking by non-R GBs after some time. Hence, the ratio of non-R GBs to R3 GBs remains about the same till the end of solidification. Of course, to better simulate the experimental results in a realistic manner, 3D simulation would be necessary, but this is beyond the scope of this study.
H.K. Lin, C.W. Lan / Journal of Crystal Growth 475 (2017) 150–157
Acknowledgement This work was supported by the Ministry of Science and Technology of Taiwan.
References [1] K. Fujiwara, W. Pan, K. Sawada, M. Tokairin, N. Usami, Y. Nose, A. Nomura, T. Shishido, K. Nakajima, Directional growth method to obtain high quality polycrystalline silicon from its melt, J. Cryst. Growth 292 (2) (2006) 282–285. [2] L. Stockmeier, G. Müller, A. Seidl, T. Lehmann, C. Reimann, J. Friedrich, Preferred grain orientations in silicon ribbons grown by the string ribbon and the edge-defined film-fed growth methods, J. Cryst. Growth 395 (2014) 74–79. [3] Y.M. Yang, A. Yu, B. Hsu, W.C. Hsu, A. Yang, C.W. Lan, Development of highperformance multicrystalline silicon for photovoltaic industry, Prog. Photovolt. Res. Appl. 23 (3) (2015) 340–351. [4] G. Stokkan, Y. Hu, Ø. Mjøs, M. Juel, Study of evolution of dislocation clusters in high performance multicrystalline silicon, Sol. Energ. Mater. Sol. C. 130 (2014) 679–685. [5] C.C. Hsieh, Y.C. Wu, A. Lan, H.P. Hsu, C. Hsu, C.W. Lan, Comparison of defect formations in solar silicon growth from small random and large oriented seeds, J. Cryst. Growth 419 (2015) 1–6. [6] K. Fujiwara, Y. Obinata, T. Ujihara, N. Usami, G. Sazaki, K. Nakajima, Grain growth behaviors of polycrystalline silicon during melt growth processes, J. Cryst. Growth 266 (4) (2004) 441–448. [7] P. Chen, Y.L. Tsai, C.W. Lan, Phase field modeling of growth competition of silicon grains, Acta Mater. 56 (15) (2008) 4114–4122. [8] T. Duffar, A. Nadri, On the twinning occurrence in bulk semiconductor crystal growth, Scripta Mater. 62 (12) (2010) 955–960. [9] H.K. Lin, C.W. Lan, Revisiting the twinning mechanism in directional solidification of multi-cyrstalline silicon sheet, Acta Mater. 131 (2017) 1–10. [10] A. Tandjaoui, N. Mangelinck-Noel, G. Reinhart, B. Billia, X. Guichard, Twinning occurrence and grain competition in multi-crystalline silicon during solidification, C. R. Phys. 14 (2–3) (2013) 141–148. [11] T. Riberi-Béridot, N. Mangelinck-Noël, A. Tandjaoui, G. Reinhart, B. Billia, T. Lafford, J. Baruchel, L. Barrallier, On the impact of twinning on the formation of the grain structure of multi-crystalline silicon for photovoltaic applications during directional solidification, J. Cryst. Growth 418 (2015) 38–44. [12] Y.T. Wong, C. Hsu, C.W. Lan, Development of grain structures of multicrystalline silicon from randomly orientated seeds in directional solidification, J. Cryst. Growth 387 (2014) 10–15. [13] V. Randle, Twinning-related grain boundary engineering, Acta Mater. 52 (14) (2004) 4067–4081. [14] B.W. Reed, M. Kumar, Mathematical methods for analyzing highly-twinned grain boundary networks, Scripta Mater. 54 (6) (2006) 1029–1033. [15] V. Randle, Y. Hu, The role of vicinal Sigma 3 boundaries and Sigma 9 boundaries in grain boundary engineering, J. Mater. Sci. 40 (12) (2005) 3243– 3246. [16] M. Trempa, I. Kupka, C. Kranert, T. Lehmann, C. Reimann, J. Friedrich, Evolution of grain structure and recombination active dislocations in extraordinary tall conventional and high performance multi-crystalline silicon ingots, J. Cryst. Growth 459 (2017) 67–75.
157
[17] H.K. Lin, M.C. Wu, C.C. Chen, C.W. Lan, Evolution of grain structures during directional solidification of silicon wafers, J. Cryst. Growth 439 (2016) 40–46. [18] W.T. Read, W. Shockley, Dislocation Models of Crystal Grain Boundaries, Phys. Rev. 78 (3) (1950) 275–289. [19] D.G. Brandon, The structure of high-angle grain boundaries, Acta Metall. 14 (11) (1966) 1479–1484. [20] M. Kohyama, R. Yamamoto, M. Doyama, Structures and energies of symmetrical <011> tilt grain boundaries in silicon, Phys. Status Solidi B 137 (1) (1986) 11–20. [21] K.G.F. Janssens, D. Olmsted, E.A. Holm, S.M. Foiles, S.J. Plimpton, P.M. Derlet, Computing the mobility of grain boundaries, Nat. Mater. 5 (2) (2006) 124–127. [22] A.Y. Badmos, H.J. Frost, I. Baker, Simulation of microstructural evolution during directional annealing with variable boundary energy and mobility, Acta Mater. 51 (10) (2003) 2755–2764. [23] H. Håkan, Influence of anisotropic grain boundary properties on the evolution of grain boundary character distribution during grain growth—a 2D level set study, Model. Simul. Mater. Sci. Eng. 22 (8) (2014) 085005. [24] N. Moelans, B. Blanpain, P. Wollants, Quantitative phase-field approach for simulating grain growth in anisotropic systems with arbitrary inclination and misorientation dependence, Phys. Rev. Lett. 101 (2) (2008) 025502. [25] T. Hirouchi, T. Tsuru, Y. Shibutani, Grain growth prediction with inclination dependence of h110i tilt grain boundary using multi-phase-field model with penalty for multiple junctions, Comput. Mater. Sci. 53 (1) (2012) 474–482. [26] T. Hirouchi, Y. Shibutani, Explicit distinctions between 2D MPF grain growth simulations and EBSD analyses to determine driving mechanism of grain growth, Mater. Trans. 54 (10) (2013) 1884–1893. [27] J. Pohl, M. Müller, A. Seidl, K. Albe, Formation of parallel (1 1 1) twin boundaries in silicon growth from the melt explained by molecular dynamics simulations, J. Cryst. Growth 312 (8) (2010) 1411–1415. [28] A. Nadri, Y. Duterrail-Couvat, T. Duffar, Two-dimensional numerical modeling of grain structure in multi-crystalline silicon ingot, J. Cryst. Growth 385 (2014) 16–21. [29] A. Karma, W.-J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. 57 (4) (1998) 4323–4349. [30] J.A. Warren, R. Kobayashi, A.E. Lobovsky, W.C. Carter, Extending phase field models of solidification to polycrystalline materials, Acta Mater. 51 (20) (2003) 6035–6058. [31] H.K. Lin, H.Y. Chen, C.W. Lan, Phase field modeling of facet formation during directional solidification of silicon film, J. Cryst. Growth 385 (2014) 134–139. [32] R. Kobayashi, J.A. Warren, W.C. Carter, Vector-valued phase field model for crystallization and grain boundary formation, Physica D 119 (1998) 415–423. [33] I. Steinbach, F. Pezzolla, B. Nestler, M. Seesselberg, R. Prieler, G.J. Schmitz, J.L.L. Rezende, A phase field concept for multiphase systems, Physica D 94 (1996) 135–147. [34] G.Y. Chen, H.K. Lin, C.W. Lan, Phase-field modeling of twin-related faceted dendrite growth of silicon, Acta Mater. 115 (2016) 324–332. [35] C.W. Lan, C.C. Liu, C.M. Hsu, An adaptive finite volume method for incompressible heat flow problems in solidification, J. Comput. Phys. 178 (2) (2002) 464–497. [36] D.L. Olmsted, E.A. Holm, S.M. Foiles, Survey of computed grain boundary properties in face-centered cubic metals-II: grain boundary mobility, Acta Mater. 57 (13) (2009) 3704–3713. [37] K. Fujiwara, M. Ishii, K. Maeda, H. Koizumi, J. Nozawa, S. Uda, The effect of grain boundary characteristics on the morphology of the crystal/melt interface of multicrystalline silicon, Scripta Mater. 69 (3) (2013) 266–269.