Cellular Automaton Modeling of the Transition of Multi-Crystalline Silicon from a Planar Faceted Front to Equiaxed Growth

A modeling approach combining the lattice Boltzmann (LB) method and the cellular automaton (CA) technique are developed to simulate the faceted front to equiaxed structure transition (FET) of directional solidification of multi-crystalline silicon. The LB method is used for the coupled calculation of velocity, temperature and solute content field, while the CA method is used to compute the nucleation at the silicon-crucible interface and on SiC particles, as well as the mechanism of growth and capturing. For silicon, the interface kinetic coefficient is rather low, which means that the kinetic undercooling can be large. A strong anisotropy in the surface tension and interfacial kinetics are considered in the model. A faceted front in conjunction with a sufficiently high carbon content can lead to equiaxed growth by nucleation on SiC. The temperature gradient in Si melt at the interface is negative, which leads to the occurrence of a faceted interface. The higher the absolute value of thermal gradients, the faster the growth velocity. Due to differences in the degree of undercooling, there will be the unification of facets in front of the solid-liquid interface. Transitions from faceted front to thermal equiaxed dendrites or faceted equiaxed grains are observed with smaller or larger impurity contents, respectively.


Introduction
Multi-crystalline solidification of silicon is a low-cost way of producing photovoltaic cells although it is less efficient than single crystalline silicon [1].The solidification process is an important step in manufacturing components.During solidification, the shape of the solid-liquid interface plays an important role in the formation of defects and the selection of grain.For silicon, the Jackson factor is greater than two [2].Therefore, the interface tends to form facets.The formation mechanism of a faceted crystal-melt interface was investigated by in situ observation [3][4][5].It has been shown that when a wavy perturbation is introduced into a planar crystal-melt interface, the perturbation would result in zigzag facets.Lan et al. [6,7] extended phase field simulation of facet formation from two to three dimensions for directional solidification of the silicon thin film.Chen et al. [8] investigated the growth of a twin-related silicon dendrite through a novel phase field model.Miller et al. [9] applied a phase field model to compute the evolution of two grains during the solidification of silicon.However, phase field cannot calculate such a large-scale microstructure growth as the transition from faceted front to equiaxed dendrites in the process of the solidification of multi-crystalline silicon ingot, yet cellular automaton (CA) method is an alternative.In this paper, therefore, we use cellular automaton (CA) method to calculate the facet growth of silicon.
The growth behavior of columnar to equiaxed transition (CET) of multi-crystalline silicon is an important phenomenon that affects the final microstructure of a multi-crystalline silicon ingot.X-ray examination have shown that a transition from planar front to equiaxed growth is likely to occur in case of faceted interfaces, observed in multi-crystalline silicon ingots [10].A few models and mechanisms have been proposed to understand the columnar to equiaxed transition (CET).Hwang [11] developed a three-dimensional cellular automaton-finite element (CA-FE) model to simulate the microstructural evolution of multi-crystalline silicon ingots.In addition, this model can represent the occurrence of equiaxed grains observed ahead of a planar faceted interface due to carbon segregation during solidification.Although cellular automaton-finite element (CA-FE) model can calculate the columnar to equiaxed transition (CET) process, it is difficult to calculate the flow field, and the model is also very complicated.However, cellular automaton-lattice Boltzmann (CA-LB) method can be used to calculate the complex melt convection, especially the fluid flow between the gaps of different grains which is difficult for finite element (FE) to compute and the model is simple.The cellular automaton-lattice Boltzmann (CA-LB) method has been widely used in the field of simulated crystal growth [12][13][14].
In general, macroscopic simulation can improve the production process, yet it cannot simulate crystal growth mechanism accurately.Microscopic simulation has a solid theoretical basis, it is subjected to the efficiency of computation and cannot simulate large-scale of crystal growth.In mesoscopic scale, the cellular automaton combining lattice Boltzmann method is a good means to simulate crystal growth.Zhao et al. [15] set up a three-dimensional cellular automaton (CA) model to simulate the microstructure change of iron carbon alloy in the solidification process.This model took into account the solid-liquid interface curvature, surface energy and dendrites anisotropy.Meanwhile, it analyzed the grain and dendrite growth under different undercooling and thermal conductivity parameter.Huang et al. [16] solved the solid-liquid phase change problems by the lattice Boltzmann (LB) method and developed a new approach to treat the latent heat source term by modifying the equilibrium distribution function for the temperature.
The coupled cellular automaton-lattice Boltzmann (CA-LB) model retains both advantages of cellular automaton (CA) and lattice Boltzmann (LB), and the model is simple and effective.Compared with traditional model, it breaks the limitations of the continuum mode.Therefore, in this work, a new cellular automaton-lattice Boltzmann (CA-LB) model was established for calculating the transition of the planar faceted front of columnar grains to equiaxed growth in the directional solidification of a multi-crystalline silicon ingot, which is on the basis of exploring the calculation of faceted interface of multi-crystalline silicon with cellular automaton (CA) method.

Description of the Model and Numerical Calculation Method
In the condition of small Peclet number and low Reynols number, a model of coupled CA and LB method was built which was used to simulate the directional solidification of multi-crystalline silicon ingot.The LB method was used for the coupled calculation of velocity, temperature and solute content field, while the CA method was used to compute the liquid/solid phase change.

LBM Model
A single relaxation time lattice Bhatnagar-Gross-Krook (LBGK) method, the D2Q9 model, were employed [17].In consideration of latent heat source term and concentration source term, the particles' space-time evolution equations used to calculate velocity, temperature and solute content field's distribution function are given as follows [18].
Appl.Sci.2017, 7, 1251 3 of 13 In above equation, r and t represent position and time, respectively.∆t is time step.e i is the speed of i direction.τ f , τ t and τ c are the relaxation time of velocity, temperature and solute content field respectively, then the relaxation time shows the speed level of distribution functions of each particles ( f i , h i and g i ) close to equilibrium distribution functions of particles ( f eq i , h eq i and g eq i ).f i , h i and g i are the distribution functions of particles of velocity, temperature and solute content field respectively.f eq i , h eq i and g eq i are the equilibrium distribution functions of velocity, temperature and solute content field, respectively.F i , H i and G i are the source terms of external force, heat and concentration, respectively.The source term of external force F i shows interaction force on particles.The heat source term H i and concentration source term G i are the latent heat release and solute redistribution during the growth of multi-crystalline silicon.
In the D2Q9 format, the discrete particle velocity e i ( i = 0 ∼ 8) in lattice point can be expressed as In above equation, c = ∆x/∆t is the grid velocity, ∆x is the space step.The macroscopic physical quantity of the fluid, such as fluid density ρ, speed u, concentration C l and temperature T can be obtained by the corresponding distribution functions of particles (Equation ( 5)-( 8)) where F is the buoyancy of liquid.According to the approximation of Boussinesq, the density in the buoyancy term is the linear function of the temperature gradient and the concentration gradient.Hence, the buoyancy term approximates to the equation of ( 9) where ρ 0 is the density when the temperature is T 0 and the concentration is C 0 .g is the acceleration of gravity.β T and β C are the expansion coefficients of temperature and concentration respectively.In Equation ( 1)-(3), the external force term F i , H i and G i can be expressed as the equation of ( 10)-( 12) where w i is the weight coefficient, ∆ f s is the increment of solid fraction, ∆H f is the molar latent heat, C p is the molar liquid heat capacity, C l is the liquid phase solute component in the unit, k is the solute partition coefficient.The value of w i is the following Equation (13).
According to the analysis of Chapman-Enskog, continuity equation, N-S equation and convection diffusion equation can be obtained from the above mentioned LB equation.Fluid dynamics viscosity ν, thermal diffusivity α and solute diffusion coefficient D are connected with fluxion relaxation time τ f , temperature relaxation time τ t and concentration relaxation time τ c respectively (Equations ( 17)-( 19)) [14].

CA Model
The CA model for crystal growth includes the geometry of a cell, the state of a cell, the neighborhood state and the transition rules during one time step.In this work, a two dimensional (2-D) CA model is described.In order to observe the transition from planar front to equiaxed growth, the computational domain size, 5 cm × 5 cm, was divided into 500 × 500 uniform square cells with a length of 100 µm .The governing equations and numerical algorithms for calculating the crystal growth of crystalline silicon are described in details below.

Nucleation
This work uses a continuous nucleation model based on Gaussian distribution [19].The density of nucleation n(∆T) can be expressed by In the equation, ∆T is undercooling, dn/d(∆T) represents an increase in nucleation density dn when the undercooling d(∆T) is increased, and the calculation of dn/d(∆T) under the Gaussian normal distribution is given by where n max is the maximum nucleation density, ∆T σ is the standard deviation and ∆T mea is the mean undercooling.In addition, then the above function is neatened as follows: where ∆T max is maximum nucleation undercooling.We have considered two different models to find nucleation points in our system.One is nucleation at silicon-crucible interface as a result of the interfacial interaction between the molten silicon and crucible.The other is nucleation on SiC particles as shown in Figure 1.In fact, during the solidification process, carbon rejection at the solid-liquid interface leads to an increase of the carbon concentration in the silicon liquid.When carbon concentration reaches the maximum solubility limit, the silicon reacts with carbon to form SiC, allowing SiC particles to precipitate into the Si-melt (i.e., the supersaturation with respect to SiC) [20].
The solubility limit of carbon in Si-melt C L (T) is expressed by a polynomial function [21]: n max is generally given in a three-dimensional form, therefore in this paper we convert it into two-dimensional: where n * S and n S are the maximum nucleation density of silicon-crucible interface for the two-dimensional and three-dimensional, respectively, n * V and n V are the maximum nucleation density of impurities sites for the two-dimensional and three-dimensional, respectively.Consequently, the probability of nucleation is written as: where δN is the total number of new grains, N is the total number of cells, V C is a cell volume.The nucleation probability P of any cell cannot be too large, its range is [0, 1).Compared probability nucleation P with a random number r(0 < r < 1), and if P > r, the cell will be nucleated, whereas the cell cannot be nucleated [19].n and V n are the maximum nucleation density of impurities sites for the two-dimensional and three-dimensional, respectively.Consequently, the probability of nucleation is written as: where N  is the total number of new grains, N is the total number of cells, C V is a cell volume.The nucleation probability P of any cell cannot be too large, its range is   1 , 0 .Compared probability nucleation P with a random number   , and if r P  , the cell will be nucleated, whereas the cell cannot be nucleated [19].

Model for the Growth of the Columnar Faceted Interface
The undercooling of the front of the solid-liquid interface consists of solute diffusion undercooling, thermodynamic undercooling and curvature undercooling.
In above equation, m is the liquidus slope, l C is the actual solute concentration of the liquid  The undercooling of the front of the solid-liquid interface consists of solute diffusion undercooling, thermodynamic undercooling and curvature undercooling.
In above equation, m is the liquidus slope, C l is the actual solute concentration of the liquid phase in the cell, C 0 is the initial concentration, T L is the liquidus temperature, T is the node temperature of the cell, Γ is Gibbs-Thomson coefficient given by γ SL ∆S f , K is the interface curvature given by the following Equation ( 28) [22].
In above equation, ε is the anisotropy strength, φ is the angle between the local growth direction and the X axis direction, θ is the crystallographic orientation.
The relationship between the growth rate and the undercooling can be written as where v e is the growth rate, ∆T is the kinetic undercooling.µ(φ) is the kinetic coefficient equation obtained by the equation of (31) [24] µ(φ) = µq(φ) (31) where µ is the angular average of the kinetic coefficient, q(φ) is the anisotropy function obtained by the equation of (32).
Here, δ and λ are parameters.
In the case of the presence of undercooling, once the cell has nucleated, it will grow with a preferential direction relying on its crystallographic orientation.The growth rate is calculated by Equation (30), and the horizontal projection length of crystal growth l(t) can be given by (33) [25] where δt is time step.The solid fraction f s (t) can be obtained with where δ x is the distance between two adjacent cells.When f s (t) ≥ 1, which means the growth front of the solid cell touches the center of its neighboring liquid cell, the neighboring cell transforms its state from liquid to solid and gets the same orientation index as this solid cell [26].
Fujiwara [27] proposed that the undercooling at the growth front is the key parameter to divide growth behaviors.The grain growth behavior is expected to be changed between ∆T = 3.8 and 6.3 K.When the moving velocity of solid-liquid growth interface is slow and the interfacial morphology is flat, the (111) grain extended to lateral direction during growth.In this simulation, the solid-liquid interface front-side undercooling fluctuates between 1.5 and 0.4 K. Its value is less than 3.8 K and the moving velocity of solid-liquid growth interface is slow and the interfacial morphology is flat.Therefore, this work mainly simulates the growth morphology of the dominant grain orientation < 111 >.
An integer in the range of [1 − q] is randomly selected and assigned to each cell after nucleation, which indicates the crystallographic index (i.e., grain orientation) of that point.Sixty crystallographic indexes (q = 60) have been considered in our model.The orientation θ i of the grain growth direction in Equation (33) varies with the angle between ±π/4, which can be calculated using the state of the cell i as follows: In order to take into account the dominant grain orientation < 111 >, each grain is allocated a misorientation θ i in 2D space where θ i is zero for < 111 > grain.Three types of cells are considered in the CA model: solid, liquid, and interface.Each cell is characterized by temperature, solute concentration, crystallographic orientation, and solid fraction.When the solid fraction of interface cell goes to 1, the attribute of this cell changes into solid state.The solid cells will capture the nearest surrounding liquid cells making them into interface cells.This research's capturing way is the 4-neighbor rule of Von-Neumann [28].

Model for the Growth of the Thermal Equiaxed Dendrites
The growth rate of the thermal equiaxed dendrites can be written as [10]: (36) where α is the thermal diffusivity, ∆S f is the molar melting entropy, ∆H f is the molar latent heat, σ is classically taken as 1/4π 2 , γ SL is the free interfacial energy, Ω is the molar volume, C p is the molar liquid heat capacity.
When SiC precipitates in the melt and the kinetic undercooling exceeds nucleation undercooling, nucleation and growth of equiaxed grains in the melt become possible.The radius of the equiaxed grains after time δt is expressed by [10] where v e is the equiaxed growth rate in Equation (36), which leads to

Model for the Growth of the Faceted Equiaxed Dendrites
In that case and to simplify the model, we use the same kinetic law for the faceted equiaxed grains as for the columnar front in Equation (30).

Boundary Conditions
In a 2D {(x, y)|0 ≤ x ≤ L, 0 ≤ y ≤ L} computational domain, the unknown distribution function on the four boundaries and solid-liquid boundary have a strong impact on the calculation of the flow field, temperature field and concentration field.Under the condition of natural convection, the boundary conditions are treated as follows.
(1) In the velocity field, four regional boundaries and solid-liquid boundary are disposed by the form of non-slipping rebound boundary [18].
(2) In the temperature field, solid and liquid have the same heat transfer and the left, right and top boundary of the field is adiabatic.The heat dissipation of the bottom is in constant heat flux.The temperature distribution functions of the four boundaries are calculated by non-equilibrium extrapolation method.(3) In calculating solute field, the solute diffusion of the solid phase is ignored.The four regional boundaries are calculated by no diffusion boundary conditions and the solid-liquid boundary is calculated by rebound form.

Validation of the Simulation of Faceted Growth by CA
In order to verify the correctness of the CA model for faceted growth, this paper simulates the faceted growth for strong anisotropy of kinetic coefficient and compares the simulation results with the literature and experimental observations.Initially, a small circular solid nucleus is placed in an undercooled liquid.The temperature is initially set as T = 1683 K in the nucleus and T = 1680 K in the surrounding undercooled liquid.Adiabatic boundary conditions are applied to temperature field.As shown in Figure 2, a facet is formed having a normal vector in the 45 • direction for q(φ) with λ = 3 and δ = 0.01.The facet forms from an initially circular crystal, and the shape is maintained while growing.In order to observe the morphological transformation from planar to zigzag facets in detail we attempted to use the computational domain size of 613 µm × 500 µm.Figure 3 shows the shape of the Si crystal-melt interfaces with λ = 3, δ = 0.01 and ε = 0.15.A wavy perturbation is introduced into a planar crystal-melt interface, and the perturbation is amplified and results in zigzag facets.The simulation results and phenomena are consistent with the literature of Takuya [24] and experimental observations [3,[34][35][36].Through the comparison with the literature and experimental observations, the CA simulation can be validated. .A wavy perturbation is introduced into a planar crystal-melt interface, and the perturbation is amplified and results in zigzag facets.The simulation results and phenomena are consistent with the literature of Takuya [24] and experimental observations [3,[34][35][36].Through the comparison with the literature and experimental observations, the CA simulation can be validated.

The Growth of the Faceted Interface
To examine the thermal gradient in silicon, we compare our simulated temperature distributions using CA-LBM with those obtained by the analytical solution in [5] for different growth velocity.Figure 4 shows the thermal fields of the Si crystal and melt during crystal growth for e v = 106, 153, and 218 s μm .The growth velocity is controlled by changing the initial undercooling.The initial undercooling of the three different growth velocity is 0.8 K, 1 K and 1.2 K respectively.As shown in Figure 4, as the growth velocity increases, the absolute value of thermal gradients at both sides of the interface increase as well.The temperature gradient in Si melt is negative and it is necessary to produce faceted interface.When growth velocity is slow, the temperature gradient is positive.Hence, the morphology of the interface remains planar.On the other hand, when growth velocity is high, the temperature gradient in the Si melt becomes negative and the perturbation introduced into the planar interface is amplified leading to interface instability and forms zigzag facets.

The Growth of the Faceted Interface
To examine the thermal gradient in silicon, we compare our simulated temperature distributions using CA-LBM with those obtained by the analytical solution in [5] for different growth velocity.Figure 4 shows the thermal fields of the Si crystal and melt during crystal growth for v e = 106, 153, and 218 µm/s .The growth velocity is controlled by changing the initial undercooling.The initial undercooling of the three different growth velocity is 0.8 K, 1 K and 1.2 K respectively.As shown in Figure 4, as the growth velocity increases, the absolute value of thermal gradients at both sides of the interface increase as well.The temperature gradient in Si melt is negative and it is necessary to produce faceted interface.When growth velocity is slow, the temperature gradient is positive.Hence, the morphology of the interface remains planar.On the other hand, when growth velocity is high, the temperature gradient in the Si melt becomes negative and the perturbation introduced into the planar interface is amplified leading to interface instability and forms zigzag facets.
In order to observe the morphological transformation from planar to zigzag facets in detail we attempted to use the computational domain size of .A wavy perturbation is introduced into a planar crystal-melt interface, and the perturbation is amplified and results in zigzag facets.The simulation results and phenomena are consistent with the literature of Takuya [24] and experimental observations [3,[34][35][36].Through the comparison with the literature and experimental observations, the CA simulation can be validated.

The Growth of the Faceted Interface
To examine the thermal gradient in silicon, we compare our simulated temperature distributions using CA-LBM with those obtained by the analytical solution in [5] for different growth velocity.Figure 4 shows the thermal fields of the Si crystal and melt during crystal growth for e v = 106, 153, and 218 s μm .The growth velocity is controlled by changing the initial undercooling.The initial undercooling of the three different growth velocity is 0.8 K, 1 K and 1.2 K respectively.As shown in Figure 4, as the growth velocity increases, the absolute value of thermal gradients at both sides of the interface increase as well.The temperature gradient in Si melt is negative and it is necessary to produce faceted interface.When growth velocity is slow, the temperature gradient is positive.Hence, the morphology of the interface remains planar.On the other hand, when growth velocity is high, the temperature gradient in the Si melt becomes negative and the perturbation introduced into the planar interface is amplified leading to interface instability and forms zigzag facets.Furthermore, to examine the facet formation mechanism in more detail, we extracted a slice from the columnar crystals for discussion.Figure 5 shows the simulation results for q(φ) with λ = 3 and δ = 0.2 in the domain of 613 µm × 500 µm .Figure 5a represents the position and shape of the solid-liquid interface at every 0.2 s.The growth velocity of the solid-liquid interface is 153 µm/s.A wavy perturbation was introduced into the interface and the perturbation grew and resulted in zigzag facets.The interface is still rounded up to 0.6 s, and facets become clear after 0.6 s.Two small facets merged into a bigger one.The undercooling of facet B is higher than that of facet A, leading to the faster growth velocity of facet B than facet A. Hence facet B can catch up with facet A and unite B to form facet C which is shown inside the small box in Figure 5a.We can further examine the temperature field at 3 s from Figure 5b.As shown, temperature rises at the solid-liquid interface due to the latent heat.The latent heat diffuses mainly toward undercooling liquid.Undercooling is higher at the peaks than that at the concave.The greater interface undercooling, the faster growth velocity of faceted interface.The interface undercooling is about 0.58 K.
Furthermore, to examine the facet formation mechanism in more detail, we extracted a slice from the columnar crystals for discussion.Figure 5 shows the simulation results for   wavy perturbation was introduced into the interface and the perturbation grew and resulted in zigzag facets.The interface is still rounded up to 0.6 s, and facets become clear after 0.6 s.Two small facets merged into a bigger one.The undercooling of facet B is higher than that of facet A, leading to the faster growth velocity of facet B than facet A. Hence facet B can catch up with facet A and unite B to form facet C which is shown inside the small box in Figure 5a.We can further examine the temperature field at 3 s from Figure 5b.As shown, temperature rises at the solid-liquid interface due to the latent heat.The latent heat diffuses mainly toward undercooling liquid.Undercooling is higher at the peaks than that at the concave.The greater interface undercooling, the faster growth velocity of faceted interface.The interface undercooling is about 0.58 K.

The Transition from Columnar Faceted Front to Equiaxed Growth
In this work, we mainly discussed different contents of carbon impurity for the faceted front to equiaxed structure transition of multi-crystalline silicon.Mangelinck [10] proposed two theoretical equiaxed crystal mode (thermal equiaxed dendrites and faceted equiaxed grains) and the reason for the two morphologies of equiaxed crystal is their low undercooling.Therefore, in this section, we changed the constitutional undercooling by different impurity contents.
The growth mode of the equiaxed Si grains was discovered.They can be faceted or rough and dendritic due to their low contents of impurity.The morphological evolution of the thermal equiaxed dendrites and faceted equiaxed grains during the process of FET are shown in Figure 6.When the contents of carbon 0 C increases from wt.% 003 .0 to wt.% 04 .0 , the morphology of equiaxed crystal changes from rough equiaxed dendrites to faceted equiaxed grains.This is because that the more contents of carbon impurity, the lower constitutional undercooling.Therefore, the impurity content is an important factor for affecting the equiaxed crystal morphology.
The transition occurs at the same position regardless of the initial carbon concentration.This is due to the fact that nucleation on SiC particles begins at a undercooling of about K 04 .0 , regardless of the initial contents of carbon.In the phase diagram [33], when the initial temperature is determined, the undercooling degree of nucleation on SiC particles is determined.

The Transition from Columnar Faceted Front to Equiaxed Growth
In this work, we mainly discussed different contents of carbon impurity for the faceted front to equiaxed structure transition of multi-crystalline silicon.Mangelinck [10] proposed two theoretical equiaxed crystal mode (thermal equiaxed dendrites and faceted equiaxed grains) and the reason for the two morphologies of equiaxed crystal is their low undercooling.Therefore, in this section, we changed the constitutional undercooling by different impurity contents.
The growth mode of the equiaxed Si grains was discovered.They can be faceted or rough and dendritic due to their low contents of impurity.The morphological evolution of the thermal equiaxed dendrites and faceted equiaxed grains during the process of FET are shown in Figure 6.When the contents of carbon C 0 increases from 0.003 wt.% to 0.04 wt.%, the morphology of equiaxed crystal changes from rough equiaxed dendrites to faceted equiaxed grains.This is because that the more contents of carbon impurity, the lower constitutional undercooling.Therefore, the impurity content is an important factor for affecting the equiaxed crystal morphology.
The transition occurs at the same position regardless of the initial carbon concentration.This is due to the fact that nucleation on SiC particles begins at a undercooling of about 0.04 K, regardless of the initial contents of carbon.In the phase diagram [33], when the initial temperature is determined, the undercooling degree of nucleation on SiC particles is determined.

Conclusions
In this paper, we simulate the transition of a planar faceted front to equiaxed growth in directional solidification of a multi-crystalline silicon ingot using a CA-LB coupling model.The modified CA model includes the interfacial kinetics and a strong anisotropy in the surface tension, in order to simulate faceted growth.The simulation results of CA-LB method applied to the faceted growth of silicon are highly consistent with the calculation of phase field model proposed by Lan et al.Our CA method is also competent to the researching of faceted growth.More importantly, the CA method has the ability to calculate large amounts of grain growth, which is fully reflected in the calculation of the transition from columnar faceted front to equiaxed growth.

Figure 1 .
Figure 1.Nucleation at silicon-crucible interface and on SiC particles.

Figure 1 .
Figure 1.Nucleation at silicon-crucible interface and on SiC particles.

2 Solute 1 .
19 × 10 −2 T − 1.21 × 10 −5 T Validation of the Simulation of Faceted Growth by CA In order to verify the correctness of the CA model for faceted growth, this paper simulates the faceted growth for strong anisotropy of kinetic coefficient and compares the simulation results with the literature and experimental observations.Initially, a small circular solid nucleus is placed in an undercooled liquid.The temperature is initially set as K 1683  T in the nucleus and K 1680  T in the surrounding undercooled liquid.Adiabatic boundary conditions are applied to temperature field.As shown in Figure2, a facet is formed having a normal vector in the 45° direction for  forms from an initially circular crystal, and the shape is maintained while growing.
Appl.Sci.2017, 7, 1251 10 of 14In order to observe the morphological transformation from planar to zigzag facets in detail we attempted to use the computational domain size of

Figure 3 .
Figure 3. Evolutions of interface morphologies for s m v e / 241  in the domain of

Figure 4 .
Figure 4. Temperature fields of Si crystal and melt during crystallization for e v = 106 s m 

Figure 3 .
Figure 3. Evolutions of interface morphologies for s m v e / 241  in the domain of

Figure 4 .
Figure 4. Temperature fields of Si crystal and melt during crystallization for e v = 106 s m 

.
Figure5arepresents the position and shape of the solid-liquid interface at every 0.2 s.The growth velocity of the solid-liquid interface is 153

Figure 5 .
Figure 5. (a) Evolution of interface morphologies at ve = 153 s m 

Figure 5 .
Figure 5. (a) Evolution of interface morphologies at v e = 153 µm/s; the time interval is 0.2 s; (b) the corresponding temperature field at 3 s.
Model for the Growth of the Columnar Faceted Interface