Phase field modeling of growth competition of silicon grains

(1)

Phase ﬁeld modeling of growth competition of silicon grains

P. Chen, Y.L. Tsai, C.W. Lan

*

Department of Chemical Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan, ROC Received 13 December 2007; received in revised form 18 April 2008; accepted 21 April 2008

Available online 5 June 2008

Abstract

The growth competition of two grains, with orientations (1 1 0) and (1 0 0), during directional solidification of silicon is simulated by using a phase field model. The two-dimensional simulations show two distinct competition mechanisms, with either interfacial or kinetic dominance, depending on the undercooling. At low undercooling, the interfacial effect is dominant so that (1 1 0) grain grows laterally, expelling the other grain. On the other hand, at high undercooling, the grain competition follows the same pattern at the beginning, but the (1 0 0) grain eventually becomes dominant, expanding its domain. In addition, the facet vanishing process and the dihedral angle evo-lution are discussed. The simulated results and phenomena are consistent with the experimental observations of Fujiwara et al. [Fujiwara K, Obinata Y, Ujihara T, Usami N, Sazaki G, Nakajima K. J Cryst Growth 2004;266:441] and the analytical predictions of Atwater et al. [Atwater HA, Thompson CV, Smith HI. J Mater Res 1988;6:1232].

Keywords: Phase ﬁeld model; Polycrystalline silicon; Solidiﬁcation; Grain selection

1. Introduction

The crystal growth of polycrystalline silicon (poly-Si) from the melt has become an important subject in the pho-tovoltaic industry because of the increasing market share of poly-Si solar cells. Because the grain size and orientation dramatically affect the minority carrier lifetime and thus the conversion efficiency, control of the grain growth behav-ior, especially grain competition, is critical[1–4]. Fujiwara et al.[1]observed the silicon melt growth starting from a bicrystal seed in (1 0 0) and (1 1 1) orientations. They found that the cooling rate affected the grain growth behaviors sig-nificantly. At low cooling rate (1 K min1), the planar interface was observed both on (1 0 0) and (1 1 1) grains. The (1 1 1) grain grew laterally to (1 0 0) diminishing the (1 0 0) grain during the solidification. However, at a much higher cooling rate (30 K min1), the (1 0 0) grain grew faster than the (1 1 1) grain and gradually eliminated the (1 1 1) grain. Fujiwara et al.[1]proposed that the behavior at the lower cooling case was dominated by the force

balance at the grain boundary tri-junction; the grain with a lower surface energy will be dragged by the higher one. For the higher cooling case, Atwater et al. [2] proposed the relative growth rates of the grains as the main factor determining growth behavior; the faster grain wins the com-petition and occupies the space during growth. In addition to the grain competition, in the experiments by Fujiwara et al., macro-steps also appeared from the (1 0 0) interface at high cooling rates, while the (1 1 1) grain remained planar. These macro-steps on the (1 0 0) grain, having various sizes and derivates from the (1 1 1) face at unsteady state, grow with time, and only those composed of the (1 1 1) face remained at steady state [3]. Simulation of these experi-ments could be useful, and provide a foundation for the large-scale modeling of poly-Si solidiﬁcation. However, the simulation of grain growth remains a great challenge in solidiﬁcation, and no such a simulation showing the key experimental features has been carried out so far.

The phase ﬁeld model (PFM) has been used extensively to study the microstructural dynamics of solidiﬁcation (e.g.

[5,6]), but generally without considering the grain bound-ary eﬀects. The modeling of polycrystalline material requires the treatment of crystallographic orientations

*

Corresponding author. Fax: +886 2 2363 3917. E-mail address:cwlan@ntu.edu.tw(C.W. Lan).

www.elsevier.com/locate/actamat

Available online at www.sciencedirect.com

(2)

and interfacial energies, but very few models have been proposed. Indeed, in polycrystalline phase field modeling, the complicated physical effect of crystalline orientation on solidification and grain interactions needs to be consid-ered. Morrin et al.[7]proposed a free energy density func-tion with n wells for n different orientafunc-tions. Steinbach et al.[8]developed a multiphase field model and succeeded in modeling the solidification of eutectic and peritectic materials. However, their models were not able to keep the rotation symmetry of the free energy potential in the modeling of grain boundary formation and dynamics. By proposing an orientation field h, which describes the local grain orientation, Kobayashi et al.[9,10]made a successful attempt, based on a rotation-invariant free energy func-tion, at the phase field modeling of grain boundaries. Warren et al.[11]further extended the model to the solid-ification of polycrystalline materials. Recently, nucleation phenomena have also been considered by Granasy et al.

[12,13]. For a material having a highly anisotropic interfa-cial energy and a tendency to form facets, such as silicon, special attention is needed. By using a regularization method, Eggleston et al. [14] have simulated successfully the crystal shape with corners and edges. Although suitable PFMs have been proposed for the solidification of poly-crystalline materials, there are some computational issues. One of the most critical problems in phase field simulation arises from the diffusive interface thickness d. The choice of dneeds to be small enough so that the sharp-interface limit can be approached. However, such a value, which is usu-ally of the order of the microscopic capillary length d0, is

too small from the simulation point of view. With the use of adaptive meshes [15,16], the computational eﬃciency can be greatly improved, and this can be particularly useful for polycrystalline solidiﬁcation having grain boundaries.

In this paper, we present adaptive phase field modeling of poly-Si solidification. During simulation, the mesh is adaptively refined or coarsened with the phase field vari-able and varivari-able gradients. The directional solidification of silicon from a bicrystal seed is focused, and the effect of undercooling is discussed. The PFM and the adaptive numerical scheme used are described briefly in Section 2. In Section3, the validation of the present PFM by an ana-lytical solution is presented through grain boundary groove simulations. Section4 is devoted to the simulations of the competitive melt growth of silicon, where the analysis and discussions of the growth behaviors are also given. Finally, short conclusions are given in Section5.

2. Phase ﬁeld model

The present PFM is developed following the work of Warren et al.[11]. The free energy function F in a constant volume V is postulated as F¼ Z V fð/; T Þ þe 2 / 2jr/j 2 þ s/2jrhj þe 2 h 2/ 2 jrhj2 " # dV ; ð1Þ

where / is a phase ﬁeld variable representing the crystalline state as / = 1 and the melt state as / = 0, and the param-eter e/is used to measure the anisotropic eﬀect. The

orien-tation ﬁeld variable h measures the local orienorien-tation of the crystal with respect to the ﬁxed coordinate system, and lies in the domain of p 6 h 6 p. The Helmholtz free energy density f(/, T) takes the form

fð/; T Þ ¼ pð/ÞT Tm Tm

þ W /2_{ð1 /Þ}2

; ð2Þ

where T is the temperature, Tmthe melting point, and W is

the height of the free energy potential barrier between bulk phases and p(/) = /3(10 15/ + 6/2). The parameter s measures the energy cost due to the misorientation at the grain boundary.

To model the anisotropic crystal shape, the gradient energy coeﬃcient e/needs to be further modiﬁed. As

pro-posed by Kobayashi [17], as a function of the crystalline orientation w, which is deﬁned as the angle between the interface normal and the reference axis, the coeﬃcient e/

or a crystal with 4-fold symmetry can be written as e/ðwÞ ¼ e/gðwÞ ¼ e/½1 þ e cos 4ðw w0Þ: ð3Þ

The crystalline orientation w is calculated by w= tan1(/y//x), where the subscript in the phase ﬁeld

variable denotes the first derivative with respect to the spa-tial coordinate. However, when the anisotropic strength e is greater than 1/15, missing orientations occur and the crys-tal shape obtained from the anisotropic function is not realistic where the reciprocal c-plot becomes concave. To overcome this, and at the same time to reveal the edges and corners for silicon crystals, we adopt the modified gra-dient energy coefficient e/ within the missing orientations

as proposed by Eggleston et al. [14]:

e/¼

e/ðwmÞ

cosðwmÞ

cos w; ð4Þ

where wmis the ﬁrst missing orientation, and is calculated

from the following relationship:

gðwmÞ sin wmþ

ogðwmÞ

ow cos wm¼ 0: ð5Þ Finally, through the variation principle that minimizes the total free energy, the evolution of the phase ﬁeld, orien-tation ﬁeld and temperature variables can be derived:

o/ ot ¼ M/ ðge/Þ 2 r2 /þ e2 / o ox g og ow o/ oy e2 / o oy g og ow o/ ox 2W /ð1 /Þð1 2/Þ ð6Þ 30M/L/2ð1 /Þ2 T Tm Tm M/2s/jrhj; /2oh ot¼ Mhr / 2 s jrhjrh Mhe2/g og ohj ~r/j 2 ; ð7Þ oT ot ¼ r ðarT Þ 30/ 2 ð1 /Þ2 L Cp o/ ot: ð8Þ

(3)

In the above equations, M/and Mhare the mobility of the

phase ﬁeld and the orientation ﬁeld variables, respectively. L is the latent heat, Cpthe heat capacity and a the thermal

diffusivity. For simplicity, we also assume equal thermal diffusivity in all phases in the simulation. To link the PFM parameters with the realistic physical properties, the parameters used in the WBM model[5]are taken: e/¼ ffiffiffiffiffiffiffiffi 6rd p ; ð9Þ W ¼ 3r d; ð10Þ M/¼ l_kTm 6dL ; ð11Þ

where r is the solid–liquid surface energy.

To facilitate the simulation, the governing equations are further represented in dimensionless form, where the length is rescaled by the characteristic length l0 and the time t

rescaled by l2₀=a. The dimensionless temperature u is deﬁned by (TTm)/Tm unless otherwise stated. Fig. 1a

shows the computational domain of the directional solidi-fication of bicrystalline silicon. A bicrystal seed with differ-ent crystalline oridiffer-entations is placed in front of the undercooled silicon melt. Because the simulation is two-dimensional, it is assumed that the (1 1 1) and (1 1 0) orien-tations are the same. The symmetry boundary condition is applied to the both sides of the domain, while at the top the temperature is given; the lower thermal boundary is set to be adiabatic. The physical properties of silicon used in the simulation are listed inTable 1. The interface thickness is taken as 4.1 108m in all simulations. According to Wheeler et al.[5], the interface thickness has to be of the

order of capillary length (of the order of 109m) to approach the sharp-interface solution. However, this will greatly increase the computational cost. Hence, a much thicker interface is adopted, and the comparisons between theory and the simulations thus are merely qualitative. Nevertheless, we anticipate that certain physical features of the experiments could be revealed through the simula-tion. The governing equations with the boundary condi-tions are solved by an adaptive ﬁnite volume method

[15,16]. A sample adaptive mesh during simulation is shown inFig. 1b. As shown, the mesh is refined to give bet-ter numerical resolution near the inbet-terfaces, while grids are coarser in the area without significant physical features. The criterion for the mesh refinement is based on the phase field variable in the range of 0.05 < / < 0.95, and the min-imum dimensionless mesh size Dx = 0.0625 is used. Before the simulation for silicon melt growth, which will be dis-cussed in Section 4, is presented, we shall first simulate grain boundary grooving problems. Through the compari-son with the analytical solutions, the adaptive PFM simu-lation can be validated.

3. Grain boundary grooving

The grain boundary groove is formed during solidiﬁca-tion between two adjacent grains. The groove shape at the junction is a force balance of the interfacial tensions. To obtain a static groove shape, a linear temperature pro-ﬁle having a gradient G of 1.25 105

K m1 is given; a domain 16 lm 8 lm is used. Because once a steady state has been reached, only the static shape is of interest, the release of latent heat does not play a role. The characteris-tic length l0is 2 lm here, which is calculated from T_GLmr

1=2

[18]. Two examples, with isotropic and anisotropic interfa-cial energies, are discussed in the following sections.

3.1. Dihedral angle

Again, as just mentioned, the dihedral angle is a result of the force balance at the tri-junction. For simplicity, we have ﬁrst assumed that the interfacial energy is isotropic. Therefore, the analytical solution of the dihedral angle hd

could be easily obtained from Young’s law as

Fig. 1. (a) A schematic directional solidiﬁcation of silicon bicrystal in undercooled melt. (b) A sample adaptive mesh in silicon growth simulation.

Table 1

Physical properties of silicon

Properties Name Value References

a Thermal diﬀusivity (m2s1) 2.8 105 _[17]

r Surface energy (J m2) 0.59a Estimated L Latent heat (J m3) 4.15 109

[17]

lj Kinetics coeﬃcient (m K1s1) 0.02b Estimated

e Anisotropy strength 0.25 Estimated

a _{The value of the surface energy is estimated from the equation}

~c¼ rðV2 mÞ

1=3_=L_{where ~}_c

Si¼0.525[15].

b_{The kinetic coeﬃcient we used is simply chosen by assumption;}

how-ever, the order of the value is close to those evaluated from the experi-ments and the simulations[16–18].

(4)

hd ¼ 2 cos1

cgb

2csl

; ð12Þ

where cgbis the grain boundary energy and cslis the solid–

liquid surface energy, i.e. r. cgbcan be calculated from the

derivation in Ref.[11]: cgb¼ ffiffiffiffiffiffiffiffiffi We/ p 3 ð1 / 3 minÞ; ð13Þ /min¼ 1 1þ s e/Dh : ð14Þ

Three cases with diﬀerent grain boundary energies are sim-ulated: (a) cgb= csl, (b) cgb= 1.5csland (c) cgb= 1.95csl. In

addition, to determine the dihedral angle from simulated results, the position of / = 0.5 is taken for measuring the angle; the dihedral angle is obtained from the intersection of the grains.Fig. 2 shows the simulation results, and the

comparison of the calculated dihedral angles with the ana-lytical solutions is shown inTable 2. As shown, they are in good agreement; the trivial diﬀerence is due to the diﬀusive interface assumed in the PFM.

3.2. Anisotropic grain boundary groove

We have obtained a good agreement between the simu-lated results and the analytical ones for the isotropic grain boundary energy. We now need to further check the agree-ment for the case with highly anisotropic surface energy. The analytical solution for the grain boundary groove with highly anisotropic surface energy is also available [18] as the following: x¼Tm GL Z w hd cslþ @2c_sl ow2 _{sin w} y dw; ð15Þ y¼ 2Tm GL cslsin wþ oc_sl ow cos w csl w¼p 2 !1=2 : ð16Þ

The above equations are solved in two intervals, i.e. hd6w 6 h1and h26w 6p₂, where h1and h2 are

orienta-tions very close to the kink orientation. In Eq. (15), the dihedral angle hdcan be obtained from Herring’s relation

[19]:

c_gb 2ocsl

ow sin w 2cslcos w¼ 0: ð17Þ In order to prevent the tilting of the grain boundary due to unequal horizontal forces at the junction, the crystal orien-tations of the grains are chosen to have the same value but with an opposite sign. Three cases with diﬀerent crystal ori-entation w0 are simulated: p/6, p/4 and p/3. The grain

boundary energy is ﬁxed at 1.99 csl, so that only the surface

energy is varied for all three cases. The anisotropic strength e= 0.08 is adopted. The reason for not using a higher e as listed inTable 1is to avoid the oscillating wave, so that a static shape can be obtained by simulation for comparison.

Fig. 3shows a comparison of the simulated results with the analytical solutions. Again, the calculated groove shapes are in good agreement with the analytical ones. This good agreement reinforces our conﬁdence in the simulation of bicrystal problems.

4. The competitive growth of silicon grains

The simulation of the competitive growth of silicon grains, as illustrated inFig. 1a, is presented in this section. The following values are chosen: coupling constant

Fig. 2. Simulations of isotropic grain boundary groove: (a) cgb= csl, (b)

cgb= 1.5csl, (c) cgb= 1.95csl. In each ﬁgure, melt is at the upper half and

the crystals are at the lower half. As predicted by Young’s law, the dihedral angle will decrease as the grain boundary energy increases.

Table 2

Comparison of calculated dihedral angles with analytical solutions Dihedral angle comparison

cgb/csl 1 1.5 1.95

Analytical 118° 78° 25°

(5)

s = 2 J m2and characteristic length l0= 2.1 106m; the

domain size is taken as 42 lm 16.8 lm. In addition, the orientation ﬁeld mobility Mh is set to 0.05 M/, because

the time scale of solidification is much smaller than the grain boundary dynamics. In other words, evolution of the crystal orientation is much slower than solidification. It should be pointed out that the domain size and the sim-ulation time are three orders of magnitude smaller than those in the real experiments [1]. In view of the diffusive

interface kinetics, the characteristic time for silicon growth is about 107s; therefore it is too costly to simulate the real experiments, which lasted a few minutes. In fact, since the physical features are the focus of this study, simulation using smaller length and time scales should reveal the same physical phenomena. Two undercoolings, D = 1.685 K/Tm

and D = 6.74 K/Tm, are simulated, and the simulated

results for the evolution of grain growth are shown inFigs. 4 and 5, respectively. As shown, the interface morphology and the grain competition are consistent with the observed ones[1]; the (1 1 0) grain prevails at low undercooling, while the (1 0 0) grain becomes dominant at high undercooling. Further discussion of the morphology and grain competi-tion is given in the following seccompeti-tions.

4.1. Morphology

The evolution of morphology for D = 1.685 K/Tm and

D= 6.74 K/Tm extracted from Figs. 4 and 5 are put

together for comparison inFig. 6a and b, respectively. As shown, with a smaller undercooling at D = 1.685 K/Tm,

the interface of the (1 1 0) grain remains planar, but small edges emerge from the (1 0 0) interface. The diﬀerence in

Fig. 3. Comparison of calculated anisotropic grain boundary groove shape with analytical solution at e = 0.08.

(6)

the morphology can be explained by the grain surface energy, which stabilizes the solid–liquid interface. Appar-ently, the (1 1 0) surface energy plays a stronger role in

stabilizing the interface. These edges appearing in the (1 0 0) surface are smaller than 1 lm, and the sizes remain about the same during solidiﬁcation at D = 1.685 K/Tm.

However, the situation is slightly diﬀerent for D= 6.74 K/Tmin that these edges coarsen into larger

fac-ets. Furthermore, the intersection angle of these edges is also changing during solidiﬁcation. As the edges become lar-ger, only those with an angle close to 70° remain, as shown in

Fig. 7a. This could be explained by the minimization of the surface energy. When the missing orientations disappear, the equilibrium edge angle is obtained. This can be further explained by the equilibrium crystal shape grown from a seed in an undercooled melt; the anisotropic strength e is greater than 1/15, so that the edge and angle appear in the range of the missing orientations[14].Fig. 7b shows a sample simula-tion of the equilibrium crystal shape. As shown, its equilib-rium edge angle is about 70°.

Interestingly, the vanishing of the facet, leading to big (1 1 0) facets, or the coarsening process, is suppressed by the surface energy at low undercooling, but becomes signif-icant at high undercooling. This is consistent with the experimental observations. Nevertheless, it should be pointed out that the crystal shape with edges obtained

Fig. 5. Simulation of Si melt growth at D = 6.74 K/Tm.

Fig. 6. Evolution of the morphology of competitive growth at diﬀerent undercoolings: (a) 1.685 K/Tmand (b) 6.74 K/Tm.

(7)

using the present anisotropic function, though similar to the experimental results, is not the case for the faceted mor-phology; they are not totally ﬂat surfaces. According to Wulﬀ’s theorem, the faceted portion of the crystal corre-sponds to the cusps in the polar plot of the surface energy, but the function we use here is a smooth function. In other words, the simulation based on the smooth function leads to curved edges instead of facets. However, the curved sur-face at high anisotropic strength looks quite similar to a facet.

4.2. Force balancing at the groove root

According to Fujiwara et al.[1], growth behavior at low undercooling is dominated by the force balance at the grain boundary root. In other word, the grains with lower sur-face energies remain. To check this, we apply Herring’s relation[19]to analyze the force balance in the groove dur-ing growth at D = 1.685 K/Tm. Herring’s relation describes

the interfacial tension balance at the junction for aniso-tropic interfacial energies:

c_gbcos b¼X 2 i¼1 c_isin w_ioci owi cos w_i ; ð18Þ cgbsin b¼ X2 i¼1 c_icos w_iþoci owi sin w_i : ð19Þ

In the above equations, the grain boundary energy is as-sumed to be isotropic, where the angle b is deﬁned as the tilt angle of the grain boundary with the x-axis. It should be pointed out that by deﬁnition the interfacial tension is equivalent to the surface energy. For the horizontal part of the force, we have taken the force as positive when it is pointing to the y-axis, which here is the direction of

the (1 0 0) interfacial tension. For the vertical part, the force is positive when it is towards the positive x-axis. The over-all interfacial tension is defined by jRF j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiF2_hþ F2

v

q

, where Fh and Fv are the overall horizontal and vertical

interfacial tensions, respectively. The decrease in the over-all interfacial tension, in other words, corresponds to the minimization of the surface energy.

Fig. 8shows the overall horizontal and vertical interfa-cial tension evolutions at D = 1.685 K/Tm; the total

interfa-cial tension is also included. As shown, the overall horizontal tension is much greater than zero, indicating that the interfacial tension on the (1 0 0) face is significantly larger than that on the (1 1 1) face. The value is decreasing with time as force relaxes, consistent with the experimental observation[1]. For the overall vertical interfacial tension, the value reduced to about0.6 J m2 (in the negative x-direction). This can be easily understood in terms of the reluctance of the grains to grow, so that the interfacial ten-sion is a suppressing force. However, because the underco-oling is sufficiently large, solidification still proceeds [14]. Nevertheless, from the plot of the overall interfacial tension shown inFig. 8, the value is decreasing with time. Since the thermal driving force remains approximately the same in the simulation, the horizontal interfacial is indeed the key factor in determining the trend of the tension evolution. In other word, the decrease in the interfacial tension indi-cates the minimization of the surface energy. Eventually, the gain having a lower surface energy becomes dominant in the growth, while the groove shape reaches an equilibrium.

A similar interfacial force balance is also carried for high undercooling at D = 6.74 K/Tm, and the interfacial tension

balance results are given inFig. 9, again for the component and overall interfacial tensions. As shown in Fig. 9, the evolution of the interfacial force balances is quite diﬀerent from the previous one. InFig. 9, the horizontal and overall interfacial tensions increase at the beginning, and then gradually decrease. According to the discussion in Ref.

[2], instead of the interfacial tension balance at the groove, the growth at high undercooling is dominated by the rela-tive growth rates among the grains. Hence, we should ﬁrst

Fig. 7. (a) Coarsening of edges at D = 6.74 K/Tm. (b) Equilibrium shape

of the crystal at e = 0.25.

(8)

address the growth rate anisotropy in the present PFM. According to the work by Kobayashi [17], modifying the gradient energy coeﬃcient as an orientation function in the PFM also leads to an anisotropic kinetic coeﬃcient, i.e.

lk ¼ 6M/e/L Tm ffiffiffiffiffiffiffi 2W p gðwÞ: ð20Þ

For a crystal with 4-fold symmetry in two-dimensions, the anisotropic kinetic coeﬃcient is maximum in the (1 0 0) direction (h = 0°) and minimum in the (1 1 0) direction (h = 45°). Atwater et al. [2] further explained that the growth direction of the grain boundary could be found from the anisotropic growth rates:

sin a¼tslower grain tfaster grain

; ð21Þ

where a is the angle between the grain boundary and the slower growing grain surface. For our case here, initially the growth is affected by the interfacial tension balance so that the grain boundary is pulled toward the right as in the case of low undercooling. However, since the effect of kinetics becomes significant at high undercooling, the faster growing grain (h1 0 0i) forces the dihedral angle to decrease. As a result, the resultant interfacial tension in-creases, instead of relaxing, until the grain boundary is tilted toward the (1 1 0) grain. The tilting to the left approx-imately follows the relationship in Eq.(21). Once the grain boundary tilting angle a is established, the interfacial ten-sion start to relax and eventually an equilibrium groove will form.

We also confirm this via force balance analysis of the experimental results in Ref.[1]. For the growth experiment at a cooling rate of 30 K min1, the growth rate of each grain was approximately constant, so that the growth direction of the grain boundary remained about the same. By applying Herring’s relation to the experimental results, we find that the overall interfacial tension decreases from about 0.7 to 0.15 J m2 within the growth period. There-fore, we can conclude that the interfacial tension balance at D = 6.74 K/Tm, initially is destroyed by the kinetic effect

and the interfacial tension balancing restarts when the kinetic eﬀect no longer dominates the growth behavior.

Furthermore, the evolution of the dihedral angle can also be examined. As we just illustrated, the interface energy eventually decreases with time in both cases. In other words, the dihedral angle between two grains should also reach an equilibrium. Since, we have no actual balance angle, we take the ﬁnal one at the end of simulation as a reference. Then, the evolutions of the angle deviating from the ﬁnal one can be examined as shown inFig. 10a and b for low and high undercoolings respectively. For low und-ercooling (D = 1.685 K/Tm), since the total process is

dom-inated by the interfacial energy, the same trend is found for the angle deviating from the ﬁnal one as shown inFig. 10a. This is also consistent with the experimental observation in Ref.[1]. On the other hand, for the high undercooling case (D = 6.74 K/Tm) as shown inFig. 10b, at the beginning, the

balance is destroyed by the kinetic eﬀect. However, once the grain competition has reached to a certain extent, the angular deviation from the ﬁnal state starts to decrease, showing the same trend as the total interfacial energy after the peak.

5. Conclusions

In this study, phase field simulation is carried out for the competitive growth of silicon grains. The orientation field variable is used to treat polycrystalline solidification, and a high anisotropic function is adopted to simulate the edges. The present model is first validated through grain boundary groove simulations. The force balance test and the equilib-rium groove shape from the simulation are in good

agree-Fig. 9. Plot of the overall force component vs. time at D = 6.74 K/Tm.

Fig. 10. Evolutions of the dihedral angle deviation from the ﬁnal state: (a) low undercooling at 1.685 K/Tmand (b) high undercooling at 6.74 K/Tm.

(9)

ment with the analytical results. Furthermore, competitive growth of silicon grains is simulated, and the competition between the interfacial and kinetic eﬀects is illustrated. The undercooling is found to be the most important factor for the competition and the simulation results, which are consis-tent with the experimental observation and model predictions. Although the present phase ﬁeld simulation is primitive and qualitative, it indeed provides deeper physical insights into the polycrystalline crystal growth and grain control.

References

[1] Fujiwara K, Obinata Y, Ujihara T, Usami N, Sazaki G, Nakajima K. J Cryst Growth 2004;266:441.

[2] Atwater HA, Thompson CV, Smith HI. J Mater Res 1988;6:1232. [3] Fujiwara K, Nakajima K, Ujihara T, Usami N, Sazaki G, Hasegawa

H, et al. J Cryst Growth 2002;243:275.

[4] Aoyama T, Kuribayashi K. Acta Mater 2000;48:3739.

[5] Wheeler AA, Boettinger WJ, MacFadden GB. Phys Rev A 1992;45:7424.

[6] Loginova I, Amberg G, Argen J. Acta Mater 2001;49:573.

[7] Morin B, Elder KR, Sutton M, Grant M. Phys Rev Lett 1995;75:2156.

[8] Steinbach L, Pezzola F, Nestler B, Sesselberg M, Prieler R, Schmitz GJ, et al. Physica D 1996;94:135.

[9] Kobayashi R, Warren JA, Carter WC. Physica D 1998;119:415. [10] Kobayashi R, Warren JA, Carter WC. Physica D 2000;140:141. [11] Warren JA, Kobayashi R, Lobkovsky AE, Carter WC. Acta Mater

2003;51:6035.

[12] Granasy L, Borzsonyi T, Pusztai T. Phys Rev Lett 2002;88:206105. [13] Granasy L, Borzsonyi T, Pusztai T. Phys Rev Lett 2007;98:035703. [14] Eggleston JJ, McFadden GB, Voorhees PW. Physica D 2001;150:91. [15] Lan CW, Liu CC, Hsu CM. J Comp Phys 2002;178:464.

[16] Lan CW, Hsu CM, Chang YC. Phys Rev E 2002;65:61601. [17] Kobayashi R. Physica D 1993;63:410.

[18] Voorhees PW, Coriell SR, McFadden GB. J Cryst Growth 1984;67:425.

[19] Herring C. The physics of powder metallurgy. New York: McGraw-Hill; 1951. p. 143.