Numerical Simulation of Multi-crystalline Silicon Crystal Growth Using a Macro–micro Coupled Method during the Directional Solidification Process

In this work, the crystal growth of multi-crystalline silicon (mc-Si) during the directional solidification process was studied using the cellular automaton method. The boundary heat transfer coefficient was adjusted to get a suitable temperature field and a high-quality mc-Si ingot. Under the conditions of top adiabatic and bottom constant heat flux, the shape of the crystal-melt interface changes from concave to convex with the decrease of the heat transfer coefficient on the side boundaries. In addition, the nuclei form at the bottom boundary while columnar crystals develop into silicon melt with amzigzag-faceted interface. The higher-energy silicon grains were merged into lower energy ones. In the end, the number of silicon grains decreases with the increase of crystal length.


Introduction
Multi-crystalline silicon (mc-Si) is a very useful material owing to its excellent integrated properties [1].Compared with other solar energy materials, mc-Si grown via directional solidification is of higher quality and cheaper than others in the large market of solar cell products [2,3].In past years, numerical simulation for silicon crystal growth has been carried out by many researchers.It is well known that the grain structure of mc-Si ingot has a strong effect on the conversion efficiency of solar cells because the grain boundaries are the locations of lattice defects which interact with minority carriers [4,5].Ma et al. designed an insulation partition in a seeded directional solidification (DS) furnace and the result shows that it can reduce the total heating power consumption as well as improve the efficiency of solar cells [6].Liu et al. designed different moving partition blocks to investigate the solidification front shape and thermal stress of quasi-single-crystal silicon ingots [7].Jun et al. studied the heat and mass transfer and solute segregation during single-crystal growth by a numerical method [8].Nadri et al. have investigated the grain structure in an mc-Si ingot to predict the evolution of grain structures [9].Lin et al. have investigated the facet formation mechanism and grain competition during the silicon crystal growth via adaptive phase field modeling, which presented the morphological evolution and thermal field near the interface as well as the comparison with experimental results [10].Their work focuses on the improvement of macro-processes or the study of micro-mechanisms.Anisotropy, faceting, and twining are important features in the crystal growth of many materials, and they cannot be overlooked when studying the microscopic growth mechanism of silicon.The phase-field method is a good method for the study of crystal micro-growth, but it is important to calculate the efficiency of computation and cannot simulate large-scale crystal growth.Macroscopic simulation of crystal growth is an important method for improving the growth process, but this method ignores some crystal growth mechanisms, such as anisotropy and faceting, etc. which result in the physical mechanism of crystal growth not being accurately studied.Mesoscale simulation combining the cellular automation and finite difference or finite element methods, which is widely used in metal casting to simulate crystal growth and crystal grain structures, efficiently bridges the gap between micro-and macroscopic scale method.
Rappaz and Gandin were the first ones to apply the cellular automaton (CA) method to predict solidification grain structures [11].The material was considered to be a group of cells, and the motion of each cell was calculated based on the behavior of their neighboring environment.The CA method is not only simple; it is also a powerful tool to achieve high level of complexity and diversity.It has been used to simulate the dendritic growth in the solidification process as well as the grain boundary motions in polycrystalline materials [12,13].Yukinobu et al. developed a 3-D cellular automaton finite difference (3-D CAFD) model to predict the solidification grain structure of Al-3wt%Si binary alloy [14].Their research results can be of considerable help for controlling grain structure and segregation in complex casting products.Bei et al. used the Cellular Automata Finite Element (CAFE) modeling technique to simulate silicon grain growth.They ignored the effect of anisotropy on silicon crystal growth in the calculation model, and compared the simulation results with real ingot grain structures [15].The simulation results were in good agreement with the reality.The CA is a mesoscale method, which neglects some details in microscopic simulation method and has important guiding significance for macroscopic simulation.
Many researchers have studied the growth process of polycrystalline silicon at the macroscopic or microscopic scale [6,10,16,17].However, little research has been done on mc-Si during the directional solidification process using the CA method in the mesoscale.In this work, we therefore employed the CA method to study the process of mc-Si grain growth, including the effect of temperature field, nucleation, columnar zone growth and grain competition.

Mathematical 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 microcosmic time step.In this work, the model needs to firstly calculate temperature distribution of macrocosmic node in order to simulate the mc-Si microstructure in the computational area.To speed up the calculation, two nested grids were employed for heat transfer and silicon growth analysis: macrocosmic grids and microcosmic grids.As shown in Figure 1, the computational area was chosen as 50 mm × 50 mm, which was meshed with quadrilateral elements.The macrocosmic scale includes 25 grids and the dimension of each grid is 10 mm × 10 mm, and each macrocosmic grid includes 1 × 10 4 microcosmic grids.Compared to other mesh plug-in attempts, the design of that grid is efficient and feasible to simulate multi-crystalline silicon crystal growth [12].

Model of Temperature Distribution
The temperature field calculation includes two parts: the macrocosmic temperature field and the microcosmic temperature field.The calculation of the macrocosmic temperature field is mainly based on the finite difference method and heat diffusion equation, while the calculation of microcosmic temperature is based on the temperature of macrocosmic node and the position of microcosmic node in macrocosmic scale.Then, the temperature of microcosmic node is updated through the temperature of its neighboring nodes in computational area.The temperature field calculation includes two parts: the macrocosmic temperature field and the microcosmic temperature field.The calculation of the macrocosmic temperature field is mainly based on the finite difference method and heat diffusion equation, while the calculation of microcosmic temperature is based on the temperature of macrocosmic node and the position of microcosmic node in macrocosmic scale.Then, the temperature of microcosmic node is updated through the temperature of its neighboring nodes in computational area.

Macrocosmic Temperature Field
To simplify the simulations, this model ignores convective and radiative heat transfer, and the heat transfer is considered only by the conduction in neighboring grids.The heat diffusion equation that controls macrocosmic temperature field during solidification can be described as and the macrocosmic time step meets the requirements where x and y are the coordinates, ∆x mac and ∆y mac are macrocosmic grid spacing, ∆t is macrocosmic time step, T is temperature, κ is thermal conductivity, ρ is density, C is heat capacity, and Q is latent heat source in crystallization area.

Microcosmic Temperature Field
As shown in Figure 2, the weighted value method was used to calculate microcosmic temperature.It is assumed that the centers of macrocosmic nodes are named A, B, C, and D, whose temperature are represented respectively by T A , T B , T C , and T D , while the areas are represented by M 1 , M 2 , M 3 , and M 4 , respectively.The temperature of microcosmic node P is expressed by The coordinates of A, B, C, and D are (X i , Y i ), P is (X P , Y P ), and the equation for calculating area is

Model of Nucleation
It is difficult to form nuclei in the bulk of the pure melt at a low undercooling, so we only consider the nearby wall nucleation and ignored internal nucleation.In order to describe the process of heterogeneous nucleation, the continuous nucleation model is used in our model, and the density

Model of Nucleation
It is difficult to form nuclei in the bulk of the pure melt at a low undercooling, so we only consider the nearby wall nucleation and ignored internal nucleation.In order to describe the process of heterogeneous nucleation, the continuous nucleation model is used in our model, and the density of nucleation n(∆T) can be expressed by where ∆T is undercooling, n max is the maximum density nuclei, ∆T σ is the standard deviation, and ∆T mea is the mean undercooling.
Then, the density nuclei of any undercooling can be given by The probability of nucleation is written as where V CA is a cellular volume.Compared with probability of nucleation p and random number r, the cellular nucleate if p > r; otherwise, it does not nucleate [11].In this work, there is a constant heat flux at the bottom that results in the temperature being gradually reduced.When the temperature is below the melting point, there will be a lot of potential nuclei randomly generated, and with the loop calculation, only some will continue to grow to be a part of nuclei, and then the cell reaches a solid state [18].A randomized nucleation model was used, and each potential crystal nucleus growth may have different growth orientations.Considering the symmetry, the angle between the growth direction of the cyrstal grain and the x-axis changes only between −π/4 and +π/4, and the state value is assigned an integer among 48 classes.The change of the angle between different grain growth directions is π/96, which can meet the needs of grain growth and transformation.

Model of mc-Si Ingot Growth
The mc-Si ingot growth is governed by where dz dt represents velocity of silicon growth, ρ S is the density of solid silicon, Q is latent heat, κ s is thermal conductivity of solid silicon and κ m is thermal conductivity of melt silicon, dT dz s is temperature gradient in solid silicon and dT dz m is temperature gradient in melt silicon.
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 the Equation ( 8), and the horizontal projection length of crystal growth, l(t), can be given by where N is the iteration number, δt is the microcosmic time step, and θ is the angle of the solid cell growth direction to the linking line between this solid cell and its liquid neighbor cell.In this work, micro and macro time are equal.The solid fraction, f s (t), can be obtained with where L is the distance between two adjacent cell.In this work, the probability of growth of cell capture for the nearest neighbor liquid cell is 59.8%, and the L is equal to ∆x mic , while the corner neighbor liquid cell is 40.2%, and the L is √ 2∆x mic when f s (t) ≥ 1, which indicates that the cell has been completely transformed into a solid state and captured the adjacent liquid cell; the solid cell stops growing and its adjacent liquid cell starts to transform liquid to solid grown in the same way [19].

Temperature Distribution
The temperature distribution has a strong influence on crystal growth velocity, interface shape and grain boundaries during the crystal growth process.A suitable temperature field may decrease the defects in silicon ingot and increase the quality of silicon wafers.Therefore, it is necessary to discuss the temperature field during the process of silicon crystal growth.
According to the Procast software, the range of the heat transfer coefficient is about 0-2000 W•m −2 •K −1 between inorganic nonmetallic materials and air.In this part, we analyzed the effect of the side boundary conditions on temperature field within the computational area.The heat transfer coefficients were 0 and 1000 W•m −2 •K −1 at the top and bottom, respectively.The temperature distribution is shown in Figures 3 and 4 with different heat transfer coefficients on the side boundaries, where the temperature isotherms of the melting point represent the crystal-melt interface at different times.The color of the isotherms reflects the temperature value, with red representing the maximum temperature and deep blue representing the minimum temperature.
As clearly shown in Figures 3a and 4a, the computational area is full of silicon melt at an early stage of crystal growth.The temperature isotherms are straight and the maximum temperature gradient appears at the bottom boundary, which is attributed to the fact that heat transfer between environment and silicon melt is more intense.It is beneficial for crystal growth in that it releases the latent heat of crystallization with a high temperature gradient.In the columnar crystal zone, heat transfer is the dominant factor for crystal growth and the isothermal surface is likely to reflect the shape of the crystal-melt interface.
As shown in Figure 3, when the heat transfer coefficient is 100 W•m −2 •K −1 on the side boundaries, the shape of the crystal-melt interface changes from straight to slightly curved because the heat transport near the side boundaries is faster than that in the bottom boundary with the increase of crystal length.It can be seen that the temperature distribution at the center point of silicon melt is more uniform when the process of solidification is finished.In this case, the quality of the mc-Si ingot in the center position is much better than that near the side and bottom boundaries where the silicon ingot having higher temperature gradient.
at different times.The color of the isotherms reflects the temperature value, with red representing the maximum temperature and deep blue representing the minimum temperature.
As clearly shown in Figures 3a and 4a, the computational area is full of silicon melt at an early stage of crystal growth.The temperature isotherms are straight and the maximum temperature gradient appears at the bottom boundary, which is attributed to the fact that heat transfer between environment and silicon melt is more intense.It is beneficial for crystal growth in that it releases the latent heat of crystallization with a high temperature gradient.In the columnar crystal zone, heat transfer is the dominant factor for crystal growth and the isothermal surface is likely to reflect the shape of the crystal-melt interface.
As shown in Figure 3, when the heat transfer coefficient is 100 W•m −2 •K −1 on the side boundaries, the shape of the crystal-melt interface changes from straight to slightly curved because the heat transport near the side boundaries is faster than that in the bottom boundary with the increase of crystal length.It can be seen that the temperature distribution at the center point of silicon melt is more uniform when the process of solidification is finished.In this case, the quality of the mc-Si ingot in the center position is much better than that near the side and bottom boundaries where the silicon ingot having higher temperature gradient.In order to get more straight crystal-melt interface, The temperature field was optimized by means of adjusting the heat transfer coefficient of the side boundaries.In shown in Figure 4, the temperature field is more suitable when the lower heat transfer coefficient is dropped to 50 W•m −2 •K −1 .In this case, not only does the isothermal surface become more uniform, but the larger temperature gradient also appears at the bottom boundary.Compared with the bottom boundary as presented in Figure 3, the heat transport is also balanced on the side boundaries.Moreover, the undercooling in the horizontal direction is slower than that at the bottom boundary.The isotherm is flatter than that with a heat transfer coefficient of 100 W•m −2 •K −1 as shown in Figure 3.The evolution of the crystal-melt interface from the central axis to the right sidewall with a different heat transfer coefficient on the side boundaries during mc-Si ingot growth can be seen from Figure 5. Based on this, we predict that a high-quality mc-Si ingot will be obtained with lower heat transfer coefficients on the side boundaries.
Figure 3, the heat transport is also balanced on the side boundaries.Moreover, the undercooling in the horizontal direction is slower than that at the bottom boundary.The isotherm is flatter than that with a heat transfer coefficient of 100 W•m −2 •K −1 as shown in Figure 3.The evolution of the crystalmelt interface from the central axis to the right sidewall with a different heat transfer coefficient on the side boundaries during mc-Si ingot growth can be seen from Figure 5. Based on this, we predict that a high-quality mc-Si ingot will be obtained with lower heat transfer coefficients on the side boundaries.The ideal mc-Si ingot consists of uniform vertical columnar crystal.The temperature gradient is the driving force for crystal growth and it is desirable to have vertical growth instead of horizontal.Hence, the vertical direction needs to maintain an appropriate temperature gradient and the temperature gradient in a horizontal direction needs to constantly decrease.Figure 6 shows the temperature gradient in a horizontal direction at different stages during the mc-Si ingot solidification.The vertical axis represents the temperature gradient and the horizontal one is the distance from central axis to the side boundaries along the crystal-melt interface (the pink line in Figure 6a).The ideal mc-Si ingot consists of uniform vertical columnar crystal.The temperature gradient is the driving force for crystal growth and it is desirable to have vertical growth instead of horizontal.Hence, the vertical direction needs to maintain an appropriate temperature gradient and the temperature gradient in a horizontal direction needs to constantly decrease.Figure 6 shows the temperature gradient in a horizontal direction at different stages during the mc-Si ingot solidification.The vertical axis represents the temperature gradient and the horizontal one is the distance from central axis to the side boundaries along the crystal-melt interface (the pink line in Figure 6a).It can be seen that the maximum negative temperature gradient in the horizontal direction usually appears at the side boundaries with the higher heat transfer coefficient of 100 W•m −2 •K −1 .The growth orientation is perpendicular to the crystal-melt interface as usual, and the interface is concave with a higher heat transfer coefficient because the ability of heat transport on the side boundaries is strong.The temperature gradient in the vertical direction is positive and the crystal grows along this direction.However, the temperature gradient in the horizontal direction is negative because the heat flux transports effectively to the side boundaries with higher heat transfer coefficients and the crystal grows faster than that in the central region.The temperature of silicon melt in the central region is higher than that near side boundaries and the heat transports from the central to the side region.Based on this condition, the negative temperature gradient is formed and its value increases constantly with the growth of the mc-Si ingot.Moreover, the crystal-melt interface becomes more and more concave to the silicon melt with higher heat transfer coefficients on the side boundaries as shown in Figure 3.
In order to get a suitable interface shape that helps to obtain a perfect silicon crystal, the heat transfer coefficient is decreased.The heat transfer becomes more balanced between central melt and side boundaries with an adjusted heat transfer coefficient where the temperature gradient in the horizontal direction is around 0.5 K/mm and the crystal-melt interface is close to flat.Hence, it is It can be seen that the maximum negative temperature gradient in the horizontal direction usually appears at the side boundaries with the higher heat transfer coefficient of 100 W•m −2 •K −1 .The growth orientation is perpendicular to the crystal-melt interface as usual, and the interface is concave with a higher heat transfer coefficient because the ability of heat transport on the side boundaries is strong.The temperature gradient in the vertical direction is positive and the crystal grows along this direction.However, the temperature gradient in the horizontal direction is negative because the heat flux transports effectively to the side boundaries with higher heat transfer coefficients and the crystal grows faster than that in the central region.The temperature of silicon melt in the central region is higher than that near side boundaries and the heat transports from the central to the side region.Based on this condition, the negative temperature gradient is formed and its value increases constantly with the growth of the mc-Si ingot.Moreover, the crystal-melt interface becomes more and more concave to the silicon melt with higher heat transfer coefficients on the side boundaries as shown in Figure 3.
In order to get a suitable interface shape that helps to obtain a perfect silicon crystal, the heat transfer coefficient is decreased.The heat transfer becomes more balanced between central melt and side boundaries with an adjusted heat transfer coefficient where the temperature gradient in the horizontal direction is around 0.5 K/mm and the crystal-melt interface is close to flat.Hence, it is desirable to grow silicon with a suitable temperature field and flat crystal-melt interface with a lower transfer coefficient of 50 W•m −2 •K −1 .

Characteristics of Silicon Crystal Growth
Based on above analysis, we know that the shape of the crystal-melt interface mainly depends on the boundary conditions.In this part, we keep adiabatic conditions on the top and side boundaries while the heat transfer coefficient of 1000 W•m −2 •K −1 is used for the bottom boundary.
Figure 7 shows the process of microstructure of silicon grain growth.Blue represents the melt, and other different colors are used to represent different crystal grains.Thus, the interface between blue and other colors is the crystal-melt interface.The nuclei usually appear at the bottom boundary.At the beginning, most nuclei grow gradually with different sizes.The shape of the crystal-melt interface is zigzag faceted in the preliminary stage and the heat of crystallization is taken away by the melt and mc-Si ingot.The heat flux transports slowly and the zigzag faceted interface propels into silicon melt with lower velocity during the crystal growth, which is in line with references [24,25].The undercooling is large at the bottom of the computational area because heat transfer is intense between the melt and bottom boundary.Therefore, many fine grains are formed at the bottom of silicon melts.The silicon grains grow quickly due to the higher transfer coefficient.Being influenced by the thermal field, the nuclei appear firstly at the center of the bottom and then the nucleation is expanded to the side boundaries until the process of nucleation is finished at the bottom of the silicon During crystal growth, the sharp corners in front of interfaces disappear gradually and the shape of the crystal interface becomes smoother.The crystal-melt interface is convex to silicon melt and the convexity decreases with the increase of silicon crystal length.Some grain boundaries disappeared because the array of atoms is similar among neighboring columnar crystal zones with similar phase, and the silicon grains with the higher energy may be merged with the lower energy ones.In addition, structure and energy fluctuations exist in the crystal-melt interface during crystal growth.It will release the latent heat in front of the crystal-melt interface during crystallization, resulting in some silicon grains with high energy which will be fused and stop the growth.
The undercooling is large at the bottom of the computational area because heat transfer is intense between the melt and bottom boundary.Therefore, many fine grains are formed at the bottom of silicon melts.The silicon grains grow quickly due to the higher transfer coefficient.Being influenced by the thermal field, the nuclei appear firstly at the center of the bottom and then the nucleation is expanded to the side boundaries until the process of nucleation is finished at the bottom of the silicon melt.
As shown in Figures 3, 4 and 7, in conditions of top adiabatic and bottom constant heat flux, the shape of the crystal-melt interface changes from concave to convex with the decrease of the heat transfer coefficient on the side boundaries.In the industrial production of polycrystalline silicon, the concave crystal-melt interface easily leads to many stresses and dislocations, which reduce the crystal quality.Conversely, the flat or a slightly convex crystal-melt interface can reduce stresses and dislocations [26], so reducing the heat transfer coefficient on side boundaries is necessary to produce high-quality polycrystalline silicon crystal.

Improvement of Silicon Crystal Growth
As mentioned above, the top and side boundaries are adiabatic.Figure 8 shows the silicon crystal growth process with a different heat transfer coefficient at the bottom boundary.Controlling the morphology of the crystal-melt interface during silicon growth is crucial for obtaining high-quality crystals, which affects the macrostructure and microstructure of the products.In the end, the mechanical, optical and electrical properties of materials are also affected.Therefore, the heat transfer coefficient at the bottom boundary is adjusted in the simulations.
As shown in Figure 8a, at 10 s, the number of nuclei with a lower heat transfer coefficient of 750 W•m −2 •K −1 at the bottom boundaries is less than those with a higher heat transfer coefficient.Because the boundary condition is adiabatic on the top and side, the heat exchange via the bottom boundary is the main source for heat transfer.It is not desirable for nucleation at the bottom boundary due to poor heat transfer ability.Therefore, the energy for nucleation is inadequate and a small number of nuclei are formed.Based on this, the heat transfer coefficient at the bottom boundary is increased in order to obtain more nuclei.It can be seen that the number of nuclei increases significantly when the heat transfer coefficient is adjusted to 1000 At 100 s, the silicon grains grow into columnar zone and they propel stably to silicon melt as evidenced in Figure 8a-c.The convexity of the crystal-melt interface and the length of the small columnar crystal zone increase with the increase of heat transfer coefficient at the bottom boundary because the growth energy is mainly provided by larger undercooling through crucible bottom.Furthermore, it increases the crystal growth velocity of the silicon melt.In practice, we not only want to improve the velocity of crystal growth, but also to get a high-quality mc-Si ingot.To do so, we can adjust the heat transfer coefficient at the bottom boundary.
A small number of nuclei were obtained and the velocity of crystal growth was slow with lower heat transfer coefficient of 750 W•m −2 •K −1 .Also, the length of small columnar crystal is short and crystal-melt interface convex to silicon melt is small during the silicon grains growth.When increasing the bottom heat transfer coefficient to 1250 W•m −2 •K −1 , although the crystal growth velocity increases with large columnar crystal formed in the top zone, most part of the mc-Si ingot is undesirable.This is because the fine-grain crystals are formed easily due to a great number of nuclei that are generated at the bottom boundary with large undercooling.In addition, the crystal-melt interface is more convex to silicon melt and large temperature gradient in horizontal direction is generated with the increase of crystal length.Also, fast growth leads to inhomogeneous crystal.Therefore, it is suitable for silicon grain growth with the heat transfer coefficient of 1000 W•m −2 •K −1 at the bottom boundary.This guarantees that reasonable number of silicon nuclei are formed at early stage and it may provide enough energy for crystal growth and relatively improvement in uniformity.As a result, a relatively perfect mc-Si ingot will be achieved where the impurities and dislocation are released with the moderate convexity of crystal-melt interface.
concave crystal-melt interface easily leads to many stresses and dislocations, which reduce the crystal quality.Conversely, the flat or a slightly convex crystal-melt interface can reduce stresses and dislocations [26], so reducing the heat transfer coefficient on side boundaries is necessary to produce high-quality polycrystalline silicon crystal.

Improvement of Silicon Crystal Growth
As mentioned above, the top and side boundaries are adiabatic.Figure 8 shows the silicon crystal growth process with a different heat transfer coefficient at the bottom boundary.Controlling the morphology of the crystal-melt interface during silicon growth is crucial for obtaining high-quality crystals, which affects the macrostructure and microstructure of the products.In the end, the mechanical, optical and electrical properties of materials are also affected.Therefore, the heat transfer coefficient at the bottom boundary is adjusted in the simulations.
As shown in Figure 8a, at 10 s, the number of nuclei with a lower heat transfer coefficient of 750 W•m −2 •K −1 at the bottom boundaries is less than those with a higher heat transfer coefficient.Because the boundary condition is adiabatic on the top and side, the heat exchange via the bottom boundary is the main source for heat transfer.It is not desirable for nucleation at the bottom boundary due to poor heat transfer ability.Therefore, the energy for nucleation is inadequate and a small number of nuclei are formed.Based on this, the heat transfer coefficient at the bottom boundary is increased in According to Figure 9, the number of silicon grains decreases with increasing crystal length.It also affects the process of grains competition and the change of grain number with different heat transfer coefficient at bottom boundary.At the final stage, only three grains exist with heat exchange coefficient of 1000 W•m −2 •K −1 .The size of the grains is in the range of 1-20 mm.Therefore, the results are reasonable, and it is similar with the theories of crystal growth.
dislocation are released with the moderate convexity of crystal-melt interface.
According to Figure 9, the number of silicon grains decreases with increasing crystal length.It also affects the process of grains competition and the change of grain number with different heat transfer coefficient at bottom boundary.At the final stage, only three grains exist with heat exchange coefficient of 1000 W•m −2 •K −1 .The size of the grains is in the range of 1-20 mm.Therefore, the results are reasonable, and it is similar with the theories of crystal growth.

Conclusions
In the present work, the process of multi-crystalline silicon ingot growth in the microstructure was investigated using the cellular automaton method, which includes temperature distribution, nucleation, columnar zone growth and grain competition.Our results are identical to the existing experimental ones.Firstly, temperature field has a significant effect on the shape of the crystal-melt

Conclusions
In the present work, the process of multi-crystalline silicon ingot growth in the microstructure was investigated using the cellular automaton method, which includes temperature distribution, nucleation, columnar zone growth and grain competition.Our results are identical to the existing experimental ones.Firstly, temperature field has a significant effect on the shape of the crystal-melt interface and the quality of the mc-Si ingot.In conditions of top adiabatic and bottom constant heat flux, the shape of the crystal-melt interface changes from concave to convex with the decrease of the heat transfer coefficient on the side boundaries.Secondly, a majority of nuclei usually appear at the bottom boundary at an early stage and decrease in the following stage, which is due to the different energy of silicon grains.The columnar zone grows stably to silicon melt.The higher-energy silicon grains will be merged by the lower-energy ones and the number of silicon grains decreases with an increase in crystal length.Finally, compared with different heat transfer coefficients at the bottom boundary, an appropriate value of 1000 W•m −2 •K −1 is obtained for pursuing a high quality mc-Si ingot.

Figure 2 .
Figure 2. Weighted value method to calculate microscopic temperature.

Figure 5 .
Figure 5.The evolution of the crystal-melt interface with different heat transfer coefficients on the side boundaries during mc-Si ingot growth.

Figure 5 .
Figure 5.The evolution of the crystal-melt interface with different heat transfer coefficients on the side boundaries during mc-Si ingot growth.

Figure 6 .
Figure 6.Distribution of temperature gradient in horizontal direction from central crystal-melt interface in different time with different heat transfer coefficient on the side boundaries: (a) marked x-direction of crystal-melt interface; (b) t = 100 s; (c) t = 200 s.

Figure 6 .
Figure 6.Distribution of temperature gradient in horizontal direction from central crystal-melt interface in different time with different heat transfer coefficient on the side boundaries: (a) marked x-direction of crystal-melt interface; (b) t = 100 s; (c) t = 200 s.

Figure 9 .
Figure 9.The grain number at different times with different heat transfer coefficients (Hc) at the bottom boundary.

Figure 9 .
Figure 9.The grain number at different times with different heat transfer coefficients (Hc) at the bottom boundary.

Table 1 .
Physical properties and calculation parameters used in the present model.