Next Article in Journal
The Study of the Mechanism of Protein Crystallization in Space by Using Microchannel to Simulate Microgravity Environment
Previous Article in Journal
High-Performance Thermal Management Nanocomposites: Silver Functionalized Graphene Nanosheets and Multiwalled Carbon Nanotube
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Cellular Automaton Modeling of Silicon Facet Formation during Directional Solidification

Hebei Engineering Laboratory of Photoelectronic Functional Crystals, School of Materials Science and Engineering, Hebei University of Technology, 8 DingZiGu 1st Road, Tianjin 300130, China
*
Authors to whom correspondence should be addressed.
Crystals 2018, 8(11), 399; https://doi.org/10.3390/cryst8110399
Submission received: 6 September 2018 / Revised: 11 October 2018 / Accepted: 16 October 2018 / Published: 23 October 2018

Abstract

:
Silicon facet formation during directional solidification is simulated by cellular automaton (CA) modeling in which anisotropic interfacial energy and kinetics are considered. Numerical simulations were performed with different anisotropy strengths of interfacial energy and they show good agreement with analytical equilibrium shapes obtained by the Gibbs-Thomson equation. We also compare our results of anisotropic kinetics with in situ observation experiments and the results of the phase model to verify the accuracy of our model. Simulation results of facet formation show that perturbation is promoted to the corner by the negative temperature gradient of the interface and the heat accumulation location leads to the disappearance of small corners.

1. Introduction

Directional solidification of silicon has a striking development in photovoltaics for its low-cost. The morphology of solid-liquid interfaces during directional solidification is the key factor of the formation of defects, which has an important role in the photon-to-electron conversion efficiency. For silicon, the morphology that tends to be faced with the Jackson factor is greater than two [1], and facet formation has attracted a widespread concern in recent years. Tokairin et al. [2] observed the zigzag facet morphology of silicon with an in situ observation system, and they considered that the negative temperature gradient in front of the crystal-melt interface is the controlling factor for the growth velocity and the wavelength of facets. Fujiwara et al. [3] found the critical growth velocities of morphological transition for silicon from planar to zigzag facets. Lin et al. [4] conducted an adaptive phase field modeling with anisotropic interfacial energy and kinetics, and they examined the necessary condition for morphological instability based on the classic theory for the reported experiments.
The Cellular automaton (CA) model is a powerful tool used to study solidification problems. Rappaz and Gandin [5] developed a cellular automata-finite element (CA-FE) model to study solidification problems. Wang et al. had developed a CA model for describing the primary dendrite spacing selection of superalloys during directional solidification process, in which a modified decentered square growth algorithm was used [6]. Dong and Lee extended Wang’s model to study the stray grain formation in cast turbine blades [7]. Sun et al. developed a lattice-Boltzmann CA model to study dendritic growth process with a forced melt convection [8].
As mentioned above, all the CA models are developed for alloy solidification. So far, none of the CA models evaluated the facet formation of silicon. To simulate it, the key factors used are the anisotropic interfacial free energy and kinetics, which has been extensively studied in the phase field model. Eggleston et al. [9] developed a modified phase field method to resolve the problem of missing orientations in the four-fold symmetry with high anisotropy of the interfacial free energy. Uehara and Sekerka [10] drove two functions, having sharp minima for special directions, to describe the anisotropic kinetic coefficient, and the result showed excellent agreement with the kinetic Wulff shapes. As for kinetic coefficient, Weinstein and Brandon [11] presented a new approach for two–dimensional computations.
For the growth of silicon, our team has already performed multiple attempts. Lian et al. [12] studied the directional solidification process of multi-crystalline silicon by the CA model with different temperature boundary conditions. Zhang et al. [13] then added velocity and solute content field to the CA model to study the facet front to equiaxed structure transition (FET) during directional solidification of multi-crystalline silicon. In this paper, we present a cellular automaton model including both anisotropies to simulate the silicon facet formation. Numerical simulations with different anisotropic interfacial free energy and kinetics were performed to compare with analytical equilibrium shapes and the results of in situ observation experiment, which show good consistency.

2. Materials and Methods

A cellular automata-lattice-Boltzmann (CA-LB) coupled modeling, having anisotropic interfacial free energy and kinetics, was presented for studying the facet formation of silicon during directional solidification process with a small Peclet number and a low Reynols number. The LB method was used in temperature field simulation and the CA method was used to simulate the phase change of crystal and melt. Different from phase field method, CA modeling is an open model, which can be modified to accommodate the evolutionary approach required. At the same time, high flexibility is a challenge to the theoretical foundation. In order to verify the usability of the CA model, we will compare the results of the experiment and simulation in the subsequent content.

2.1. LBM Model

A single relaxation time lattice Bhatnagar-Gross-Krook (LBGK) method, the D2Q9 model, was employed [14]. In consideration of the latent heat source term, the particles’ space-time evolution equation used to calculate the temperature field’s distribution function is given as follows [15].
f i ( r + e i Δ t , t + Δ t ) = f i ( r , t ) + 1 τ f ( f i e q ( r , t ) f i ( r , t ) + F i ) ,
In the above equation, r and t denote position vector and time, respectively. Δ t denotes time step. e i denotes the speed of i direction. τ f denotes the relaxation time of temperature, then the relaxation time shows the speed level of distribution functions for particles ( f i ) close to equilibrium distribution functions of particles ( f i e q ). f i denotes the distribution functions of particles of temperature. f i e q denotes the equilibrium distribution functions of temperature. F i denotes the source terms of heat, which denote the latent heat release during the growth of multi-crystalline silicon.

2.2. CA Model

To simulate crystal growth, a two-dimensional (2-D) CA modeling is employed. In order to observe the facet formation, the computational domain size, 500 μ m × 500 μ m , was divided into 500 × 500 uniform square cells with a length of 1 μ m . The governing equations and numerical algorithms for the models of growth and capture to simulate silicon facet formation are introduced in detail below.

2.2.1. Model for the Growth

The undercooling of the front of the crystal-melt interface includes thermodynamic undercooling and curvature undercooling.
Δ T = T L T + Γ K f ( θ , φ ) ,
where T L denotes the liquidus temperature, T denotes the node temperature of the cell, Γ denotes Gibbs-Thomson coefficient, K denotes the interface curvature given by the following Equation (3) [16].
K = [ 2 f s x f s y 2 f s x y ( f s x ) 2 2 f s x y ( f s y ) 2 2 f s x y ] [ ( f s x ) 2 + ( f s y ) 2 ] 3 2 ,
In order to accurately calculate the partial derivative of the solid fraction, a bilinear interpolation method is used to reduce the anisotropy of tetragonal grid [17]. As shown in Figure 1, f i ( i = 1 ~ 25 ) denotes the solid fraction of each cell, and F i ( i = 1 ~ 12 ) denotes the linearly interpolated variable obtained by averaging the solid fractions of four cells surrounding a mesh node, such as F 1 = ( f 2 + f 3 + f 7 + f 8 ) / 4 . Thus, the partial derivatives of the solid fraction can be calculated by the finite-difference method based on the discretization of bilinear interpolation. For example, the first and second partial derivatives of cell f 13 in Figure 1 can be solved as
f ¯ 8 = ( F 1 + F 2 + F 4 + F 5 ) / 4 ,
f ¯ 12 = ( F 3 + F 4 + F 7 + F 8 ) / 4 ,
f ¯ 14 = ( F 5 + F 6 + F 9 + F 10 ) / 4 ,
f ¯ 18 = ( F 8 + F 9 + F 11 + F 12 ) / 4 ,
f s x = ( f ¯ 14 f ¯ 12 ) / ( 2 Δ x ) ,
f s y = ( f ¯ 8 f ¯ 18 ) / ( 2 Δ x ) ,
2 f s x y = ( F 4 + F 9 F 5 F 8 ) / ( Δ x ) 2 ,
2 f s x 2 = ( f ¯ 14 + f ¯ 12 2 f ¯ 13 ) / ( Δ x ) 2 ,
2 f s y 2 = ( f ¯ 8 + f ¯ 18 2 f ¯ 13 ) / ( Δ x ) 2 ,
f ( θ , ϕ ) is the interfacial anisotropy function given by the following Equation (13).
f ( θ , φ ) = 1 + ε cos ( 4 ( φ θ ) ) ,
In the above equation, ε denotes the anisotropy strength, ϕ denotes the angle between the local growth direction and the X axis direction, θ denotes the crystallographic orientation.
The relationship between the growth rate and the undercooling can be written as Fujiwara’s form [18,19]
υ e = μ ( φ ) Δ T ,
where υ e denotes the growth rate, Δ T denotes the kinetic undercooling. μ ( ϕ ) denotes the growth rate coefficient obtained by Equation (15)
μ ( φ ) = μ ¯ ( a n 2 ( φ ) + W 0 a n ( φ ) t n ( φ ) ) 1 ,
where μ ¯ denotes the angular average of the kinetic coefficient, a n ( ϕ ) denotes the anisotropic function for interfacial free energy and t n ( ϕ ) denotes the anisotropic function for kinetic coefficient, W 0 denotes a parameter which is related to mean diffusivity and kinetic [20].
a n ( φ ) = 1 + ε cos ( 4 φ ) ,
t n ( φ ) = { 2 [ 1 ( cos ( 4 φ ) 1 2 ) δ ] } λ ,
Here, δ and λ are parameters. For silicon, the kinetic coefficient of (100) is three-order of magnitude larger than (111) for two-dimensional crystal growth [21]. We set λ = 10 , which means that the maximum of kinetics is 1024 times as large as the minimum in the anisotropy function.

2.2.2. Model for the Capture

The cells are divided into three states for different solid fractions in cellular automaton modeling. The cell is pure solid state or liquid state when solid fraction f s ( t ) = 1 or f s ( t ) = 0 , respectively, and it is a mixed liquid-solid state when 0 < f s ( t ) < 1 .
The growth rate is calculated by Equation (17), and the horizontal projection length of crystal growth l ( t ) can be given by Equation (18) [12]
l ( t ) = t t + δ t υ e ( Δ T ( t ) ) d t cos θ + sin θ ,
Where δ t denotes time step. The solid fraction f s ( t ) can be obtained with
f s ( t ) = l ( t ) δ x ,
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. Limited neighbor solid fraction (LNSF) method [22] is used to determine whether the state transforms. If the average solid fraction of liquid cell and its neighboring cells is f ¯ s > f s L N S F , the state of the neighboring liquid cell transforms from liquid to solid. In this paper, f s L N S F is 0.225.

2.3. The Physical Parameters of Silicon Melt

The physical properties used in the present computations are listed in Table 1 [23,24].

3. Results and Discussion

3.1. Anisotropy of Interfacial Free Energy

The analytical equilibrium shape of two-dimensional domains of a phase in a matrix can be described by the Gibbs-Thomson equation in parametric form [25]:
x = υ m μ e μ 0 ( γ cos φ γ φ sin φ ) ,
y = υ m μ e μ 0 ( γ cos φ + γ φ sin φ ) ,
where υ m denotes the molar volume, γ denotes the interfacial energy and the subscripts ϕ denotes differentiation with respect to ϕ .
γ ( φ ) = γ 0 ( 1 + ε cos ( 4 φ ) ) ,
Where γ 0 is the interfacial energy of a flat interface.
Under high anisotropy conditions ( ε > 1 / 15 ), the equilibrium shape has corners and missing orientations. Figure 2a shows the analytical equilibrium shape for five degrees of anisotropy without removal of missing orientations, the part shown in red line. For ε = 0 , 0.05 , the shapes are smooth without missing orientations.
For the missing orientations, we adopt a new interfacial anisotropy function, which was proposed by Chen and Lan to adopt negative stiffness for facet simulation [26]. a is a parameter related to negative stiffness and n is a parameter related to ε ( n = 9.5 for ε = 0.25 ). As shown in Figure 2b. The result ( a = 1 ) obtained by this function are consistent with Eggleston’s function. When a = 1.5 , it makes the anisotropy stronger.
α s = ( n x n + n y n + n z n ) a / n ,
To examine the effect of the anisotropy function for interfacial free energy, numerical simulations for four degrees of anisotropy ( ε = 0 , 0.05 , 0.15 , 0.5 ) were performed. Figure 3 shows the simulation results for ε = 0.5 . Figure 3a shows the numerical equilibrium shapes, with a circular initial shape and the comparison of numerical and analytical interface positions for the quardrant 0 ϕ π / 2 . The analytical interface positions are obtained by Equations (20) and (21). In Figure 3b the isochorones of interface are constructed from a contour of constant solid fraction, f s ( t ) = 0.5 .
Due to the anisotropy interfacial energy, the temperature distribution is different in each normal direction. Figure 4 shows the result of temperature with ε = 0.5 . Figure 4a is the temperature distribution in the entire simulation area, which shows the different temperature gradients. Figure 4b is the temperature distribution in the normal direction of 0° and 45°, in which the position of star is the solid-liquid interface. The release of latent heat of crystallization leads to the negative temperature gradient in front of the solid–liquid interface, which increases the instability of the interface. The heat from 0° direction inhibits the formation of the negative temperature gradient in 45° direction leading to the formation of facet growth.
Similarly Figure 5 shows the numerical results for ε = 0 , 0.05 , 0.15 respectively. As shown in Figure 5, there is a high consistency between numerical and analytical interface positions. When the value of ε is small, the result is nearly circular. With ε increases, a nearly flat face form becomes larger in the normal direction of 45 ° . For ε = 0 , 0.05 , the shapes are smooth without corners and missing orientations which is consistent with the analytical results showed in Figure 2 and the results obtained by Eggleston’s phase field model [9]. Figure 5c shows the effect of highly anisotropic interfacial energy. The equilibrium shape is almost a square.

3.2. Anisotropy of Kinetics

To examine the effect of anisotropy function for kinetics, anisotropy function for interfacial free energy is assumed to be a constant ( ε = 0 ).

3.2.1. Effect of λ

Figure 6a–d show the results for λ = 0 , 1 , 4 , 10 with δ = 4 , respectively, with a circular initial shape. With λ = 0 , the function for kinetics is isotropic and the result is nearly circular. When λ > 0 , the effect of anisotropic is obvious. With λ increases, a nearly flat face form becomes larger in the normal direction of 45 ° . As shown in Table 2, the tip angles of the corner are 127 ° , 110 ° and 90 ° for λ = 1 , 4 , 10 , respectively, which is in good agreement with the experiment results obtained by Tokairin and Fujiwara’s in situ observation system [25].

3.2.2. Effect of δ

Figure 7a shows the effect of δ . Figure 7b–d show the simulation results for δ = 4 , 16 , 64 , respectively, with a circular initial shape, where λ = 10 in every case. As shown in Table 3, the tip angles of the corner are 90 ° , 118 ° and 124 ° for δ = 4 , 16 , 64 , respectively, which is in good agreement with the results obtained by phase field model [4,26].
The differences between CA and phase field are caused by two points. First of all, they have different types of interfaces. Phase field is diffuse interface, whose boundary layer has a thickness containing several nodes (solid fraction between 0 and 1). Additionally, CA is a sharp interface, whose solid fraction has a mutation between 0 and 1. This means that CA has lower precision in the process of derivatives. Secondly, phase field is a highly centralized solution model, which applies environmental parameters directly to obtain results. Additionally, CA is an open model, which is a dynamic growth process involving transformation and capture.

3.3. Simulations for Facet Formation

According to the previous comparison, a computation was performed using a sinusoidal interface as an initial condition with a = 1.5 , n = 9.5 , δ = 64 and λ = 10 , see Figure 8. The anisotropic interfacial energy and kinetics were chosen to make sure that corners could form during the directional solidification process. At the initial stage of growth, the growth rates of the top and bottom of the crystal-melt interface are very quick which meets the local equilibrium condition of regularization of the gradient energy. Coarsening and annexation occur in the growth process as time passes by. Figure 8 also shows that corners will be introduced by the method over time and only stable orientations are contained in the final interface.

4. Conclusions

Facet formation of silicon was simulated by the CA model by introducing anisotropy function for interfacial free energy and kinetic coefficient. A comparison between the results of anisotropic interfacial energy with Gibbs-Thomson equation shows a good agreement. For the anisotropic kinetic coefficient, the accuracy of our model compared with in situ observation experiment and the results of phase model is acceptable, whose differences should be caused by the anisotropy of the mesh and the capture rules even if we use the bilinear interpolation discretization mesh and Limited neighbor solid fraction (LNSF) method to reduce it.
The next steps of our research are three aspects: (1) Use the adaptive-meshing technique to arise the accuracy of our model. (2) Enlarge the model to three-dimensional model by parallel computing technique. (3) Add velocity and solute content field to our model.

Author Contributions

R.L. and H.C. conceived and designed the experiments; J.W. and W.M. performed the experiments and analyzed the data; N.L. and W.Y. supervised the work. All the authors have contributed to manuscript revision.

Funding

This research was funded by the National Natural Science Foundation of China, grant number No. 51475138, and the “Blue Fire Plan” Program of the Ministry of Education, grant number No. 2014-LHJH-HSZX-011.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Jackson, K.A. Crystal growth kinetics. Mater. Sci. Eng. 1984, 65, 7–13. [Google Scholar] [CrossRef]
  2. Tokairin, M.; Fujiwara, K.; Kutsukake, K.; Kodama, H.; Usami, N.; Nakajima, K. Pattern formation mechanism of a periodically faceted interface during crystallization ofsi. J. Cryst. Growth 2010, 312, 3670–3674. [Google Scholar] [CrossRef]
  3. Fujiwara, K.; Gotoh, R.; Yang, X.B.; Koizumi, H.; Nozawa, J.; Uda, S. Morphological transformation of a crystal–melt interface during unidirectional growth of silicon. Acta Mater. 2011, 59, 4700–4708. [Google Scholar] [CrossRef]
  4. Lin, H.K.; Chen, H.Y.; Lan, C.W. Phase field modeling of facet formation during directional solidification of silicon film. J. Cryst. Growth 2014, 385, 134–139. [Google Scholar] [CrossRef]
  5. Rappaz, M. Modelling of microstructure formation in solidification processes. Int. Mater. Rev. 1989, 34, 93–124. [Google Scholar] [CrossRef]
  6. Wang, W.; Lee, P.D.; Mclean, M. A model of solidification microstructures in nickel-based superalloys: Predicting primary dendrite spacing selection. Acta Mater. 2003, 51, 2971–2987. [Google Scholar] [CrossRef]
  7. Yang, X.L.; Dong, H.B.; Wang, W.; Lee, P.D. Microscale simulation of stray grain formation in investment cast turbine blades. Mater. Sci. Eng. A Struct. Mater. 2004, 386, 129–139. [Google Scholar] [CrossRef]
  8. Sun, D.; Zhu, M.; Pan, S.; Raabe, D. Lattice boltzmann modeling of dendritic growth in a forced melt convection. Acta Mater. 2009, 57, 1755–1767. [Google Scholar] [CrossRef]
  9. Eggleston, J.J.; Mcfadden, G.B.; Voorhees, P.W. A phase-field model for highly anisotropic interfacial energy. Physica D 2001, 150, 91–103. [Google Scholar] [CrossRef]
  10. Uehara, T.; Sekerka, R.F. Phase field simulations of faceted growth for strong anisotropy of kinetic coefficient. J. Cryst. Growth 2003, 254, 251–261. [Google Scholar] [CrossRef]
  11. Weinstein, O.; Brandon, S. Dynamics of partially faceted melt/crystal interfaces i: Computational approach and single step–source calculations. J. Cryst. Growth 2004, 268, 299–319. [Google Scholar] [CrossRef]
  12. Lian, Q.; Liu, W.; Li, R.; Yan, W.; Liu, C.; Zhang, Y.; Wang, L.; Chen, H. Numerical simulation of multi-crystalline silicon crystal growth using a macro–micro coupled method during the directional solidification process. Appl. Sci. 2016, 7, 21. [Google Scholar] [CrossRef]
  13. Zhang, Y.X.; Li, R.; Wang, J.; Wang, L.X.; Yan, W.B.; Liu, C.C.; Chen, H.J. Cellular automaton modeling of the transition of multi-crystalline silicon from a planar faceted front to equiaxed growth. Appl. Sci. 2017, 7, 13. [Google Scholar]
  14. Eshraghi, M.; Jelinek, B.; Felicelli, S.D. Large-scale three-dimensional simulation of dendritic solidification using lattice boltzmann method. JOM 2015, 67, 1786–1792. [Google Scholar] [CrossRef]
  15. Sun, D.K.; Zhu, M.F.; Pan, S.Y.; Yang, C.R.; Raabe, D. Lattice boltzmann modeling of dendritic growth in forced and natural convection. Comput. Math. Appl. 2011, 61, 3585–3592. [Google Scholar] [CrossRef]
  16. Pan, S.; Zhu, M. A three-dimensional sharp interface model for the quantitative simulation of solutal dendritic growth. Acta Mater. 2010, 58, 340–352. [Google Scholar] [CrossRef]
  17. Wei, L.; Lin, X.; Wang, M.; Huang, W. A cellular automaton model for the solidification of a pure substance. Appl. Phys. A 2010, 103, 123–133. [Google Scholar] [CrossRef]
  18. Fujiwara, K.; Obinata, Y.; Ujihara, T.; Usami, N.; Sazaki, G.; Nakajima, K. Grain growth behaviors of polycrystalline silicon during melt growth processes. J. Cryst. Growth 2004, 266, 441–448. [Google Scholar] [CrossRef]
  19. Fujiwara, K.; Obinata, Y.; Ujihara, T.; Usami, N.; Sazaki, G.; Nakajima, K. In-situ observations of melt growth behavior of polycrystalline silicon. J. Cryst. Growth 2004, 262, 124–129. [Google Scholar] [CrossRef]
  20. Karma, A.; Rappel, W.J. Quantitative phase-field modeling of dendritic growth in two and three dimensions. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 1998, 57, 4323–4349. [Google Scholar] [CrossRef]
  21. Buta, D.; Asta, M.; Hoyt, J.J. Kinetic coefficient of steps at the si(111) crystal-melt interface from molecular dynamics simulations. J. Chem. Phys. 2007, 127, 5262. [Google Scholar] [CrossRef] [PubMed]
  22. Lei, W.; Xin, L.; Meng, W.; Huang, W. Orientation selection of equiaxed dendritic growth by three-dimensional cellular automaton model. Physica B 2012, 407, 2471–2475. [Google Scholar] [Green Version]
  23. Qiu, S.; Wen, S.; Fang, M.; Zhang, L.; Gan, C.; Jiang, D.; Tan, Y.; Li, J.; Luo, X. Process parameters influence on the growth rate during silicon purification by vacuum directional solidification. Vacuum 2016, 125, 40–47. [Google Scholar] [CrossRef]
  24. Desai, P.D. Thermodynamic properties of iron and silicon. J. Phys. Chem. Ref. Date 1986, 15, 967–983. [Google Scholar] [CrossRef]
  25. Burton, W.K.; Cabrera, N.; Frank, F.C. The growth of crystals and the equilibrium structure of their surfaces. Philos. Trans. R. Soc. Lond. A 1951, 243, 299–358. [Google Scholar] [CrossRef]
  26. Chen, G.Y.; Lan, C.W. Understanding the facet formation mechanisms of si thin-film solidification through three-dimensional phase-field modeling. J. Cryst. Growth 2017, 474, 166–170. [Google Scholar] [CrossRef]
Figure 1. The scheme of the bilinear interpolation discretization mesh for calculating partial derivatives of the solid fraction.
Figure 1. The scheme of the bilinear interpolation discretization mesh for calculating partial derivatives of the solid fraction.
Crystals 08 00399 g001
Figure 2. (a) Analytical equilibrium shape for different values of ε ; (b) the effect of different interfacial anisotropy functions.
Figure 2. (a) Analytical equilibrium shape for different values of ε ; (b) the effect of different interfacial anisotropy functions.
Crystals 08 00399 g002
Figure 3. Simulation results for ε = 0.5 . (a) numerical equilibrium shapes and comparison of numerical and analytical interface positions; (b) isochorones of interface at every 2000 time steps.
Figure 3. Simulation results for ε = 0.5 . (a) numerical equilibrium shapes and comparison of numerical and analytical interface positions; (b) isochorones of interface at every 2000 time steps.
Crystals 08 00399 g003
Figure 4. Simulation results of temperature: (a) entire simulation area; (b) normal direction of 0° and 45°.
Figure 4. Simulation results of temperature: (a) entire simulation area; (b) normal direction of 0° and 45°.
Crystals 08 00399 g004
Figure 5. Simulation results for equilibrium shapes with different values of ε and the comparison of numerical and analytical interface position: (a) ε = 0 ; (b) ε = 0.05 ; (c) ε = 0.15 .
Figure 5. Simulation results for equilibrium shapes with different values of ε and the comparison of numerical and analytical interface position: (a) ε = 0 ; (b) ε = 0.05 ; (c) ε = 0.15 .
Crystals 08 00399 g005
Figure 6. Simulation results for different values of λ : (a) λ = 0 ; (b) λ = 1 ; (c) λ = 4 ; (d) λ = 10 ; δ = 4 .
Figure 6. Simulation results for different values of λ : (a) λ = 0 ; (b) λ = 1 ; (c) λ = 4 ; (d) λ = 10 ; δ = 4 .
Crystals 08 00399 g006
Figure 7. Simulation results for δ : (a) the kinetic anisotropic functions with different values of δ ; (b) the simulation shape with δ = 4 ; (c) δ = 16 ; (d) δ = 64 ; λ = 10 .
Figure 7. Simulation results for δ : (a) the kinetic anisotropic functions with different values of δ ; (b) the simulation shape with δ = 4 ; (c) δ = 16 ; (d) δ = 64 ; λ = 10 .
Crystals 08 00399 g007
Figure 8. Isochorones of interface at every 200 time steps.
Figure 8. Isochorones of interface at every 200 time steps.
Crystals 08 00399 g008
Table 1. Physical properties and calculation parameters used in the present model.
Table 1. Physical properties and calculation parameters used in the present model.
PropertyValue
Melting temperature, T m e l ( K ) 1683
Density of silicon solid, ρ s ( k g m 3 ) 2330 − 2.19 × 10−2T
Density of silicon melt, ρ m ( k g m 3 ) 2330 − 2.19 × 10−2T − 1.21 × 10−5T2
Thermal diffusivity of silicon solid, a s ( m 2 / s ) 9.6 × 10−6
Thermal diffusivity of melt, a l ( m 2 / s ) 2.134 × 10−5
Molar liquid heat capacity, c p ( J m o l 1 K 1 ) 26.6
Molar melting entropy, Δ S f ( J m o l 1 K 1 ) 29.8
Free interfacial energy, γ S L ( J / m 2 ) 0.438
Kinetic coefficient, μ ( m K 1 s 1 ) 2.4 × 10−4 (Estimated)
Molar latent heat, Δ H f ( J / m o l ) 4.444 × 104
Liquidus slope, m ( K / m a s s % ) −0.7637
Table 2. Comparison of calculated tip angles with Fujiwara’s experiment results.
Table 2. Comparison of calculated tip angles with Fujiwara’s experiment results.
Tip Angle Comparison
λ 1410
Experiment-110°90°
CA simulation127°110°90°
Table 3. Comparison of calculated tip angles with Lin’s phase field model.
Table 3. Comparison of calculated tip angles with Lin’s phase field model.
Tip Angle Comparison
δ 41664
Phase field126°108°90°
CA simulation124°118°90°

Share and Cite

MDPI and ACS Style

Wang, J.; Li, R.; Li, N.; Yan, W.; Ma, W.; Chen, H. Cellular Automaton Modeling of Silicon Facet Formation during Directional Solidification. Crystals 2018, 8, 399. https://doi.org/10.3390/cryst8110399

AMA Style

Wang J, Li R, Li N, Yan W, Ma W, Chen H. Cellular Automaton Modeling of Silicon Facet Formation during Directional Solidification. Crystals. 2018; 8(11):399. https://doi.org/10.3390/cryst8110399

Chicago/Turabian Style

Wang, Jia, Ri Li, Ning Li, Wenbo Yan, Wang Ma, and Hongjian Chen. 2018. "Cellular Automaton Modeling of Silicon Facet Formation during Directional Solidification" Crystals 8, no. 11: 399. https://doi.org/10.3390/cryst8110399

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop