- Email: [email protected]
Optimization of operating parameters in polysilicon chemical vapor deposition reactor with response surface methodology
Accepted Manuscript Optimization of Operating Parameters in Polysilicon Chemical Vapor Deposition Reactor with Response Surface Methodology Li-sha An, Chun-jiao Liu, Ying-wen Liu PII: DOI: Reference:
S0022-0248(18)30083-6 https://doi.org/10.1016/j.jcrysgro.2018.02.030 CRYS 24499
To appear in:
Journal of Crystal Growth
Received Date: Revised Date: Accepted Date:
23 November 2017 12 February 2018 19 February 2018
Please cite this article as: L-s. An, C-j. Liu, Y-w. Liu, Optimization of Operating Parameters in Polysilicon Chemical Vapor Deposition Reactor with Response Surface Methodology, Journal of Crystal Growth (2018), doi: https:// doi.org/10.1016/j.jcrysgro.2018.02.030
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Optimization of Operating Parameters in Polysilicon Chemical Vapor Deposition Reactor with Response Surface Methodology Li-sha An, Chun-jiao Liu, Ying-wen Liu* Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049 P R China
Abstract: In the polysilicon chemical vapor deposition reactor, the operating parameters are complex to affect the polysilicon’s output. Therefore, it is very important to address the coupling problem of multiple parameters and solve the optimization in a computationally efficient manner. Here, we adopted Response Surface Methodology (RSM) to analyze the complex coupling effects of different operating parameters on silicon deposition rate (R) and further achieve effective optimization of the silicon CVD system. Based on finite numerical experiments, an accurate RSM regression model is obtained and applied to predict the R with different operating parameters, including temperature (T), pressure (P), inlet velocity (V), and inlet mole fraction of H2 (M). The analysis of variance is conducted to describe the rationality of regression model and examine the statistical significance of each factor. Consequently, the optimum combination of operating parameters for the silicon CVD reactor is: T=1400 K, P=3.82 atm, V=3.41 m/s, M=0.91. The validation tests and optimum solution show that the results are in good agreement with those from CFD model and the deviations of the predicted values are less than 4.19 %. This work provides a theoretical guidance to operate the polysilicon CVD process. Keywords: Polysilicon; Optimal design; Response Surface Methodology; Numerical method *Corresponding author. Tel.: (86)13087588436 Email address: [email protected] 1
1. Introduction Minimum greenhouse gas emission, sustainability, economy, and easy storage of solar energy are desired properties to form a potential solution to the pending global energy and environmental crisis [1-3]. The dramatic rise of photovoltaic application opens up prospects for the raw silicon material production [4]. It is reported that the current annual production capacity of polysilicon reaches 350 kiloton, in contrast with the 26 kiloton per year in the past [5]. The cost of terrestrial solar cells can be reduced substantially, if solar cells are fabricated from solar grade silicon, which is mainly produced by chemical vapor deposition (CVD) process. Therefore, large-scale production, high quality, and low cost are acknowledged as the most important indexes in practical silicon (Si) chemical vapor deposition (CVD) system. Furthermore, the crucial technique and thermal economic evaluation of the overall CVD process are inherent of the gaseous precursor reduction to solid silicon, which widely occurs in Siemens process [6]. The formation of solid silicon during Siemens process occurs by three steps. Firstly, volatile chlorosilanes are synthesized from metallurgical grade silicon. Then, they are purified by distillation, mixed with hydrogen, and finally fed into the CVD reactor to deposit solid silicon on high-purity silicon rods at a high temperature about 1398 K [7-9]. A variety of by-products are generated in the deposition process, such as HCl, SiCl2, SiH2Cl2, SiHCl, SiCl4, SiH2, Si2H2Cl4, and SiCl6, which will be returned to the recycled system and used to reproduce Trichlorosilane (TCS). The CVD reaction will be ceased till a certain size of silicon rod and the silicon rods will be removed off the reactor after cooling down [11]. 2
Purities in the range of 99.99 % are within reach, but the polysilicon deposition rate is so slow that the deposition process always takes 60-80 hours when the silicon rod grows to certain diameter (150-200 mm) [10]. Furthermore, the intricate situation results in limited conversions of reactants and huge energy dissipation, especially when the heating is performed by electricity [12]. This attracts more attention on raising the silicon deposition rate, thus increasing the efficiency of silicon production and reducing the energy loss of Si CVD process. Examples include changing the geometric structure and controlling operating parameters to obtain the better production conditions, which will be of great significance for the improvement of Si CVD system. By analyzing and optimizing the operating parameters, such as reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M), many researchers have recently succeeded in enhancing the silicon CVD performance. D Cai et al. [13] investigated the silicon production process in a horizontal tubular CVD reactor through computational and experimental methods. They revealed the correlation between the polysilicon rate and the process flow rate, mass fraction of the silane gas in hydrogen, and the temperature of the substrate silicon tube. G. del Coso et al. [14] discussed the effects of gas velocity, the mixture of gas composition, the reactor pressure, and the surface temperature on the silicon growth rate, deposition efficiency, and energy loss. They also provided the favorable range of these operating parameters for silicon deposition performance. In addition to operating parameters, Fang et al. [15] discussed the influence of deposition time on the silicon growth based on numerical models. Historically, research on polycrystalline silicon was preceded by investigating single factor without consideration of interaction between operating parameters. 3
This implies that using single-factor analysis is very difficult to obtain accurate optimal solution, because the coupling interaction between the parameters is too complex and the huge time has been cost to finish massive amounts of experiments. So, it is very important to solve the coupling problem of many parameters and realize the optimization as quickly and effectively as possible. Then the interest in silicon CVD system was renewed with further optimization by seeking a fast and effective approach, which can consider the coupled impact between different parameters and have the characteristics of statistical analysis and accurate prediction. As one of the most efficient methods, the Response Surface Methodology (RSM) introduced by Box and Wilson is widely adopted to investigate the coupled impact between different parameters, to derive correlations of multiple factors and to optimize the system performance [16-19]. It has been widely applied in numerous manufacturing fields for the design, development, and formulation of new products, as well as improving existing product designs. Based on this method, the regression equation of factors and responses can be obtained. Furthermore, intuitive demonstration of the relationship between multiple factors with tiny amounts of experimental data is another advantage. Therefore, RSM can be applied to solve the complicated problem of coupling interaction of CVD system and further improve its performance. This article establishes a two-dimensional model of heat and mass transfer coupling with complex chemical reaction and further verifies the model by comparing with the results in the literature. Based on RSM, 30 tests are conducted, and the regression formula is obtained to accurately predict the average silicon deposition rate under any combination of different 4
parameters. The variance analysis (ANOVA) is accomplished to identify the goodness of regression formula and analyze influential weight of each parameter quantitatively. The influence of parameters, involving reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M) on the average silicon deposition rate has been presented in detail by using 2D contours and 3D surface plot. Finally, the superior average silicon deposition rate is obtained, which are help for the optimization of actual silicon production.
2. Model 2.1 Geometric model The silicon chemical vapor deposition reactions carry out in a two-dimensional plane plates as illustrated in Fig. 1. It consists of a deposition surface, an inlet, an outlet, and an adiabatic wall in the model. The length of the deposition surface is 1.5 m and the height of the inlet and outlet is 0.1 m. The mixtures of trichlorosilane (SiHCl3) and hydrogen (H2) enters the reactor at the left inlet. Then the gases react in the reactor depositing the desired solid silicon on the surface of the substrate. Finally, the reactant and product gases leave the reactors through the right outlet.
Fig. 1. Computational domain.
2.2 Governing equations and boundary conditions
5
Focusing on the deposition process that the polysilicon deposition rate in the order of microns per minute is very little. Thus, at each instant in the process, the system can be assumed to be in the steady state. The wall (ad) is set to be adiabatic, where neither mass nor heat can pass through. The deposition surface is defined as a constant temperature. At the inlet, a velocity boundary is imposed, in which the uniform velocity, temperature and species concentration are defined. An outflow condition is assigned at the outlet plane. The governing equations for continuity, momentum, energy and species in the computational domain are expressed as follows: Continuity equation:
ui 0 xi
(1)
where
ui is the velocity in the i direction. is the fluid density which is expressed as:
1
(2)
mi' i'
i'
Momentum equation: p ij ui u j gi Fi x j xi x j
(3)
Energy equation:
T (ui ( E p)) (keff h ' J ' u j ( ij )eff Sh xi xi xi j' j j
(4)
where keff is the effective thermal conductivity coefficient (including the molecular and turbulent thermal conductivity). The effect of heat transfer by radiation on the substrate 6
temperature is ignored. J j ' is the diffusion flux of component
j ' expressed as equation (5).
S h represents the consumption and production of heat due to the chemical reactions. J i' Di' ,m i Sct Di' ,m
mi' xi
(5)
1 X i' j , j i '
'
'
(6)
X j' / Di' j '
[T 3 ( Di' , j ' 0.0188
1 1 12 )] M i' M j '
(7)
pop i2' , j ' D
T ( / k )i' j'
TD*
1 2
(8)
i j ( i j ) ' '
Sct
'
(9)
'
Dt t
(10)
where Di' ,m is the diffusion flux of species i and
Sct is the turbulent Schmidt number.
Species transport equation:
( Yi ) J i Ri Si
(11)
N
Ri M i Rir
(12)
i 1
where Ri is the net rate of production of species i by chemical reaction, Mi is the molecular weight of species i, and Rir is the Arrhenius molar rate of reaction/destruction of species i in reaction r.
7
2.3 Gaseous and surface reaction mechanisms The detailed chemical reaction mechanism including Pauline Ho’s detailed [20] gas phase reaction and a set of semi-empirical surface reactions has been taken into account to study the characteristics of silicon deposition process. The equations of reactions and the kinetic parameters are summarized in Table 1. Table 1 The gas phase and surface reaction mechanism of trichlorosilane and hydrogen system. No.
Reaction
A
E （cal•mol-1）
G1 G2 G3 G4 G5 G6 G7 G8 G9 S1 S2 S3 S4 S5 S6 S7 S8 S9 S10
SiHCl3=HCl+SiCl2 SiH2Cl2=SiCl2+H2 SiH2Cl2=SiHCl+HCl H2ClSiSiCl3=SiCl4+SiH2 H2ClSiSiCl3=SiHCl3+SiHCl H2ClSiSiCl3=SiH2Cl2+SiCl2 Si2HCl5=SiCl4+SiHCl Si2HCl5=SiHCl3+SiCl2 Si2Cl6=SiCl4+SiCl2 SiHCl3+4Si(s)→Si(b)+SiH(s)+3SiCl(s) SiH2Cl2+4Si(s)→Si(b)+2SiH(s)+2SiCl(s) SiCl4+4Si(s)→Si(b)+4SiCl(s) SiCl2+2Si(s)→Si(b)+2SiCl(s) 2SiCl(s)+Si(b)→SiCl2+2Si(s) H2+2Si(s)→2SiH(s) 2SiH(s)→H2+2Si(s) HCl+2Si(s)→SiH(s)+SiCl(s) SiH(s)+SiCl(s)→2Si(s)+HCl SiHCl+2Si(s)→Si(b)+SiH(s)+SiCl(s)
3.162e14 3.162e13 6.918e14 1.585e13 3.162e13 6.310e13 5.012e13 7.943e13 1.585e14 4.1e-4 1.2e-3 3.0e-5 2.2 4.758e22 0.1 8.606e23 0.07 4.0e25 0.06
72900 69300 75800 55500 49800 44300 52300 45900 48800 -3800 -3800 -3800 15000 67000 17300 57100 5000 71500 5000
A is pre-exponential factor. E is the reaction activation energy.
2.4 Numerical scheme and grid analysis The preprocessor GAMBIT 2.4.6 is used to generate the computational meshes. When the cell size is reduced further, the average deposition rate of silicon changes less than 0.87 % indicating that the numerical results are independent of the grid size with the grid number of 8
45000. FLUENT 14.0 program is devoted to process the present simulation, combined with CHEMKIN introducing detailed chemical reaction mechanism to FLUENT. Reynold numbers in different CVD reactors are always different. In our plane reactor, the Reynold number range from 2500 to 15000, which means that the turbulence flow exists in two-dimensional model. otherwise, in view of the fact that the complex chemical reaction under high temperature happens in the reactor and the concentration of the substrate changes too largely, the turbulence intensity in the CVD system will be increased. Thus, the realizable k-ε turbulent model was adopted to simulate the turbulent transport [21-25]. The reaction model is the eddy-dissipation concept (EDC) model which is an extension of the eddy-dissipation model including detailed chemical mechanisms in turbulent flows [26, 27].
2.5 Validation of model The heat and mass transfer model coupling with gas and surface reaction mechanisms was validated by comparing with the experimental data from the literature [28]. Referring to Habuka’s research, chemical vapor deposition reactions of polysilicon took place in a horizontal reactor. The gas mixture (TCS, H2) was introduced into the reactor under atmospheric pressure with an inlet velocity of 0.67 m/s. The chemical vapor deposition reactions occur on the surface of a disk substrate held horizontally on a fixed susceptor. The substrate was heated to 1398 K using an infrared furnace through a quartz wall reactor. Fig. 2. shows the comparison of experimental and simulated silicon deposition rate under different inlet mole fractions of hydrogen. It can be seen that the simulated results agree well with the experimental values in literature [28]. So this model can accurately describe the polysilicon CVD process. 9
Silicon deposition rate /μm·min
-1
7 6
Experimental datas Simulation result
5 4 3 2 1 0 0.97
0.98
0.99
1.00
Hydrogen mole fraction
Fig. 2. Deposition rate variation of silicon along with the change of hydrogen mole fraction.
3. Response surface methodology 3.1 Experimental design This article aims at developing fast and effective evaluation procedure for investigating the influence of operating parameters (reaction temperature, pressure, inlet velocity, and inlet mole fraction of H2) on the polysilicon deposition rate, thus achieving the better silicon deposition rate. RSM is employed to build accurate relationship between various parameters and responses through multi-factor nonlinear regression. Generally, a second-order model is adopted and given as follows: k
k 1 k
k
y 0 i xi ij xi x j ii xi2 i 1
i 1 j 2
(13)
i 1
where y is the response of the system (the average deposition rate of Si), xi and xj are coded as independent variables; k is the number of variables;
0 , i , ii and ij
(i=0, 1, 2, … k;
j=0, 1, 2, … k) are the regression coefficients for the intercept, linear, quadratic and interaction terms respectively; and
is the statistical error. The flow chart of the evaluation
procedure of this approach is illustrated in Fig. 3. 10
Start
Objective definition Experimental design BBD Fluent calculation
More design points
ANOVA Quality of fit (R2)
No
Yes Response surface analysis
End
Fig. 3. Flow chart of the evaluation procedure based on RSM.
Firstly, in the optimization procedure, the average silicon deposition rate is selected as the response, and four operating parameters including: reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M) are chosen as the variable factors. According to actual reaction furnace, the effective domains of four factors are determined. Next the experimental matrix shown in Table 2 is designed and generated by the Box-Behnken Design (BBD). Table 2 Ranges and levels of independent variables in BBD. Levels Variables -1
0
1
A, T (K)
1250
1350
1450
B, P (atm)
1
3
5
C, V (m/s)
0.5
2.5
4.5
D, M
0.80
0.88
0.96
According to the design matrix arranged by BBD, the model based on ANSYS 14.0 mentioned above is performed to obtain the corresponding responses. The response results under different combinations of four operating parameters in silicon chemical vapor 11
deposition system are tabulated as shown in Table 3. The first column of the table shows the run number of experiments. The next four columns represent the actual conditions of runs and the last column represents the results of experiments.
Table 3 Design of experimental matrix and its responses. Run
T(K)
P (atm)
V (m/s)
M
Deposition rate (μm/min)
1
1350
3
2.5
0.88
10.43
2
1450
3
4.5
0.88
41.76
3
1250
5
2.5
0.88
2.55
4
1350
5
2.5
0.96
9.87
5
1250
3
2.5
0.80
1.77
6
1350
3
2.5
0.88
10.43
7
1450
1
2.5
0.88
17.47
8
1350
1
2.5
0.80
5.84
9
1250
3
4.5
0.88
2.58
10
1350
3
2.5
0.88
10.43
11
1350
1
4.5
0.88
7.22
12
1350
5
4.5
0.88
16.54
13
1350
1
2.5
0.96
4.81
14
1250
3
0.5
0.88
1.42
15
1350
3
0.5
0.96
4.00
16
1450
3
2.5
0.96
20.38
17
1450
3
0.5
0.88
14.05
18
1250
3
2.5
0.96
1.81
19
1350
3
2.5
0.88
10.43
20
1250
1
2.5
0.88
1.15
21
1450
5
2.5
0.88
40.80
22
1450
3
2.5
0.80
35.75
23
1350
3
0.5
0.80
7.94
24
1350
5
0.5
0.88
7.78
25
1350
3
2.5
0.88
10.43
26
1350
1
0.5
0.88
4.16
27
1350
3
4.5
0.96
10.04
28
1350
3
4.5
0.80
12.79
29
1350
5
2.5
0.80
14.15
12
3.2 Analysis of variance (ANOVA) Analysis of variance is carried out to evaluate the significance of the regression models and influence weight of factors. It discusses the variance of variable and obtains the evaluation index of variance. The ANOVA results of the quadratic model for average silicon deposition rate are shown in Table 4. It suggests that the model was significant, as it is evident from the model’s F-value with a low probability value (P<0.0001). Values of “prob>F” less than 0.05 indicates that model terms are significant. In this case, A, B, C, D, AB, AC, AD and A2 are statistically significant for the average silicon deposition rate. Other interaction term like BC with larger p-values can be considered to have minor influence on the average silicon deposition rate. The goodness of fitting formula is further checked by the correlation coefficient (R2) between the experimental model and predicted values of the response variable. A high R2 value of 0.9730 implies that the regression model is statistically significant and only 2.7 % of the total variations is not explained by the model. The predicted correlation coefficient (pred. R2=0.9536) also shows good agreement with the adjusted correlation coefficient (adj. R2=0.9459). In order to obtain the accurate regression model between factors and responses, the significant and insignificant terms are all considered. The final mathematical formula in terms of actual factors is shown as follows: Average Si deposition rate= 371.76439-1.10227*A-30.39349 *B-44.24245*C+911.04461*D+ 0.027420*A*B+0.03318*A*C0.48168*A*D-0.35643*B*C5.07413*B*D-1.84225*C*D+ 55.5310E-004*A2-0.15365*B20.21824*C2-158.31953*D2
13
(14)
Table 4 ANOVA for average silicon deposition rate. source
Sum of squares
df
Mean square
F value
p-value prob>f
model
3226.51
14
230.46
35.99
< 0.0001
Significant
A-T
2104.96
1
2104.96
328.70
< 0.0001
Significant
B-P
217.21
1
217.21
33.92
< 0.0001
Significant
C-V
221.84
1
221.84
34.64
< 0.0001
Significant
D-M
62.34
1
62.34
9.73
0.0075
Significant
AB
120.30
1
120.30
18.78
0.0007
Significant
AC
176.15
1
176.15
27.51
0.0001
Significant
AD
59.40
1
59.40
9.27
0.0087
Significant
BC
8.13
1
8.13
1.27
0.2788
BD
2.64
1
2.64
0.41
0.5315
CD
0.35
1
0.35
0.054
0.8192
A2
198.44
1
198.44
30.99
< 0.0001
B2
2.45
1
2.45
0.38
0.5461
2
C
4.94
1
4.94
0.77
0.3945
D2
6.66
1
6.66
1.04
0.3251
Residual
89.65
14
6.40
Lack of fit
89.65
10
8.97
Pure error
0.000
4
0.000
Cor total
3316.16
28
R-squared=0.9730
Significant
Adj R-squared=0.9459
3.3 Analysis of regression formula The model of least squares is typically used to estimate the regression coefficients in a multiple linear regression model. The residuals from the least squares fit play an important role in judging model adequacy. The normal probability plot obtained from RSM for the average silicon deposition rate is illustrated in Fig. 4. It can be seen that the residuals plot approximately along a straight line, indicating that the errors are normally distributed. Besides, almost all the internally studentized residuals are in small range near to zero, 14
meaning that the transformation of the response may provide a better analysis.
Fig. 4. Normal probability plot of residuals for average silicon deposition rate.
The relationship between the average silicon deposition rate obtained from our model and the value predicted by regression model is illustrated in Fig. 5. In this figure, the design points are distributed evenly on the diagonal line, and the error between the predicted values from regression formula and the actual CFD results is very little. It signifies that the regression formula can be used to predict the average silicon deposition rate accurately.
Fig. 5. Regression model predictions versus our model’s results for average silicon deposition rate.
4 Results and discussion From the mathematical formula, it is obvious that there exists the coupling effect of operating parameters on the silicon deposition rate. The 2-D contour plots and 3-D surface 15
plots present the interaction effects between any two factors involving reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M). Finally, the relationship between the parameters and silicon deposition rate had been revealed. Four parameters had been optimized to improve the performance of silicon chemical vapor deposition system. 4.1 Influence of temperature and pressure The combined effect of temperature and pressure on the average silicon deposition rate, at inlet velocity of 2.5 m/s and the inlet H2 mole fraction of 0.88, is shown in Fig. 6 (a) and Fig. 6 (b), respectively. From Fig. 6, we can note that compared with the pressure, the reaction temperature plays the major role in affecting the silicon deposition rate as evidenced by the larger lateral gradient in 2D contour plot. It is apparent that the average silicon deposition rate increases as the reaction temperature rises. This is because the enlarging reaction temperature promotes the silicon CVD reaction, which corresponds to the property of endothermic reactions. In addition, higher temperature further results in better crystalline properties of silicon and more metallic silicon surfaces. While, it is important to note that, under the operating condition of high temperature, the silicon will be easily polluted by equipment material and can further result in an increase of phosphorus and boron impurities because of the enhanced chemical activity of silicon. Meanwhile, the high reaction temperature also enhances the silicon corrosion reactions, seriously inhibiting the Si deposition reactions and reducing the quality of silicon. Therefore, the reaction temperature of actual reduction furnace is generally controlled at the range from 1353 to 1423 K. From Fig. 6, it also can be observed that the average silicon deposition rate increases slowly with the increase of pressure. It reveals that the silicon deposition rate is little affected 16
by pressure. The reason is that the raising pressure increases the crash chance for SiHCl3 and H2, and then improves the silicon CVD reactions. however, with low operating reaction temperature, chemical reacting rate is so low that only an amount of feed gas can effectively participate in the complex reaction. though the feed gas becomes more and more with the increase of pressure, the reaction rate is mainly determined by the reaction temperature and chemical mechanism, and the effect of the pressure is less than that of them. Otherwise, considering that the rising pressure will aggravate the equipment load, the pressure from 2 to 4 atm is conductive to the silicon CVD system in actual reduction furnace.
(a) 2-D contour plot (b) 3-D surface plot Fig. 6. Effect of temperature and Pressure on the average silicon deposition rate.
4.2 Influence of inlet velocity and inlet mole fraction of H2 The combined effect of inlet velocity and inlet mole fraction of H 2 on the average silicon deposition rate, at the temperature of 1350 K and the pressure of 3 atm, is shown in Fig. 7 (a) and Fig. 7 (b), respectively. Firstly, it can be observed that the silicon deposition rate increases with inlet velocity rising due to the sufficient reactant induced by the high velocity. But under the mole fraction of 0.8, the relative increment of silicon deposition rate is only 60 % while inlet velocity 17
increases from 1 m/s to 3 m/s. the relative increment of the silicon deposition rate gradually decreases along with the increase of the mole fraction. The reason is that as the inlet velocity continues to increase, the severely reduced reaction time becomes the dominant factor which limits the conversion rate of SiHCl3 to solid silicon, resulting in wasting a lot of reaction gases and electric energy. So, appropriate inlet velocity should be chosen to reduce costs and make rational use of resources. Secondly, the silicon deposition rate reduces as the inlet mole fraction of H 2 increases. The reason is that the rising concentration of H2 will dilute the concentration of SiHCl3 and then reduce the crash probability of SiHCl3 molecule with silicon rod surfaces. In actual reduction furnace, the inlet mole fraction of H2 from 0.89 to 0.95 is selected to reduce cost and save energy.
(a) 2-D contour plot (b) 3-D surface plot Fig. 7. Effect of inlet velocity and inlet mole fraction of H2 on the average silicon rate.
4.3 Influence of reaction temperature and inlet velocity The combined effect of reaction temperature and inlet velocity on the average silicon deposition rate, at the pressure of 3 atm and the inlet H2 mole fraction of 0.88, is shown in Fig. 8 (a) and Fig. 8 (b) in terms of 2-D contour plot and 3-D surface plot, respectively. In this case, the average silicon deposition rate is below 5 μm/min when the temperature is lower 18
than 1300 K, where the reaction temperature is the dominant factor that affect the silicon CVD reaction. High reaction temperature and inlet velocity are conducive to the achieve better silicon deposition rate.
(a) 2-D contour plot (b) 3-D surface plot Fig. 8. Effect of temperature and inlet velocity on the average silicon rate.
4.4 Influence of temperature and inlet mole fraction of H2 The combined effect of temperature and inlet mole fraction of H2 on the average silicon deposition rate, at the pressure of 3 atm and inlet velocity of 2.5 m/s, is shown in Fig. 9 (a) and Fig. 9 (b), respectively. The average silicon deposition rate is less than 5 μm/min in the area where the temperature below 1300 K. Under this condition, the reaction temperature is the main factor corresponding to the analysis above. The silicon deposition rate can be improved by higher temperature and lower inlet mole fraction of H2. It is apparent that the temperature has a more significant effect on the average silicon deposition rate compared with the inlet mole faction of H2.
19
(a) 2-D contour plot (b) 3-D surface plot Fig. 9. Effect of temperature and inlet mole fraction of H2 on the average silicon rate.
4.5 Influence of pressure and inlet velocity The combined effect of pressure and inlet velocity on the average silicon deposition rate, at the temperature of 1350 K and the inlet H2 mole fraction of 0.88, is shown in Fig. 10 (a) and Fig. 10 (b) in terms of 2-D contour plot and 3-D surface plot, respectively. The average silicon deposition rate can be enhanced by increasing the inlet velocity and pressure in the reactor. Furthermore, the different influencing rule of them shows that the inlet velocity has more significant effect on the silicon deposition rate by comparing with the effect of pressure.
(a) 2-D contour plot (b) 3-D surface plot Fig. 10. Effect of pressure and inlet velocity on the average silicon rate.
4.6 Influence of pressure and mole fraction of H2 The combined effect of pressure and mole fraction of H2 on the average silicon 20
deposition rate, at the temperature of 1350 K and inlet velocity of 2.5 m/s, is shown in Fig. 11 (a) and Fig. 11 (b), respectively. The heat and mass transfer characteristics of the silicon chemical vapor deposition process can be enhanced by increasing the pressure and reducing the mole fraction of H2, thus improving silicon deposition rate.
(a) 2-D contour plot (b) 3-D surface plot Fig. 11. Effect of pressure and inlet mole fraction on the average silicon rate.
In conclusion, appropriate design of operating parameters should be selected to control the silicon chemical vapor deposition reactions. Among operating parameters, the influence of reaction temperature on the average silicon deposition rate is most significant. Next is the influence of inlet velocity and the pressure. The influence of H 2 mole fraction is least. In actual operation of silicon chemical vapor deposition system, in order to satisfy the energy saving, high quality of polysilicon and high efficiency, the operating manual should be designed in this sequence. According to the 2-D contour plot and 3-D surface response, higher temperature, higher pressure, larger inlet velocity and little inlet mole fraction of H 2 are beneficial to improve the performance of silicon chemical vapor deposition system.
4.7 Confirmation tests and optimization In this chapter, four experimental combinations including the optimal one are performed 21
to further prove the adequacy of the regression formula obtained by RSM. The average silicon deposition rates and comparison of four tests obtained by regression model and CFD model are shown in Table 5. From this table, the maximum deviation between numerical solutions and predicted values of the silicon deposition rate is 4.19 %, which is acceptable for engineering application. It indicates that the regression model can accurately predict the average silicon deposition rate within the range of given parameters of the silicon CVD process. In addition, the confirmation test 1 is the optimal operating parameter combination, whose average silicon deposition rate increase by 12 percent compared with the maximum value in numerical experiments matrix. Further analyzing and comparing the operating parameters between them, the reaction temperature and inlet velocity of optimal combination are less than these of best combination of numerical experiments matrix. It means that after adopting optimal combination, it is significative to decrease the energy consumption and ensure safe operation of silicon chemical vapor deposition system.
Table 5 Validation results for the average silicon deposition rate. Design parameters
Average Si deposition rate
Confirmation runs
Optimal
RSM
Fluent
/(m/min)
/(m/min)
0.91
45.86
46.78
2.2%
2.5
0.88
18.69
19.39
4.19%
2
3.5
0.9
9.24
9.42
2.00%
1
1.5
0.86
2.93
2.34
-2.00%
T/K
P/atm V/m/s
1
1400
3.82
3.41
2
1400
3
3
1350
4
1300
M
Error *
5. Conclusion In this study, a two-dimensional numerical model was developed to simulate the 22
characteristics of silicon chemical vapor deposition. The numerical model was validated by experimental data. Response surface methodology was selected to establish mathematical models of the average silicon deposition rate as a function of the operating parameters involving temperature, pressure, inlet velocity, and inlet mole fraction of H 2. The developed models were verified by the residual distribution between predicted and actual results. The operating parameters effects and their interactions on the average silicon deposition rate can be explained in detailed by ANOVA. It was found that four operating parameters have different significance on the average silicon deposition rate. Furthermore, confirmation tests were carried out to verify the accuracy of the regression model and the optimized parameters was determined. Results showed a good agreement between numerical and predicted results. This study could provide an efficient and accuracy method to guide and optimize the productive process of CVD.
Acknowledgements: This work was supported by National Key R&D Program of China (Grant No. 2016YFE0204200) and the National Science Fund from the National Natural Science Foundation of China (No. 51576150)
Nomenclature： Abbreviations ANOVA
Analysis of variance
BBD
Box-Behnken Design 23
CVD
Chemical vapor deposition
CFD
Computational Fluid Dynamics
EDC
Eddy-dissipation concept
PV
Photovaltaic
RSM
Response Surface Methodology
TCS
Trichlorosilane
Symbols A
Pre-exponential factor.
Di'm
Diffusion flux of species i
E
Reaction activation energy (cal/mol)
J j'
The diffusion flux of component
kef
Effective thermal conductivity coefficient
k
Number of variable
m
Quality (kg)
M
Molar mass (g/mol)
M
Inlet mole fraction of H2
P
Pressure (atm)
qw
Heat flux (w/g)
R2
Correlation coefficient
Ri
The net rate of production of species i by chemical reaction
Rir
The Arrhenius molar rate of reaction/destruction of species i in reaction r
Sh
The consumption and production of heat due to the chemical reactions
j'
24
Sct
Turbulent Schmidt number
T
Temperature (K)
u
Fluid velocity (m/s)
V
Inlet velocity
xi,xj
Independent variable
y
Response
Greek
The temperature-dependent collision integral
Tangential force (N/m2)
The average collision diameter
Statistical error
Fluid density (kg/m3)
0
Regression coefficient for intercept
i
Regression coefficient for linear
ii
Regression coefficient for quadratic
ij
Regression coefficient for interaction terms
Reference [1]
A. Müller, M. Ghosh, Sonnenschein R, P. Woditsch, Mater. Sci. Eng. B-Adv., 2006, 134(2): 257-262.
[2]
A.F.B. Braga, S.P. Moreira, P.R. Zampieri, J.M.G. Bacchin, P.R. Mei, Sol. Energ. Mat. Sol. C. 2008, 92(4): 418-424.
[3]
P. Woditsch, W. Koch. Sol. Energ. Mat. Sol. C. 2002, 72(1): 11-26.
[4]
W.O. Filtvedt, M. Javidi, A. Holt, Sol. Energ. Mat. Sol. C. 94 (12) (2010) 1980-1995. 25
[5]
S. Ravasio, M. Masi, C. Cavallotti, J. Phys. Chem. A 2013, 117(25): 5221-5231.
[6]
B. Hazeltine, C. Fero, W. Qin. Photovoltaics World, 2010, 3: 16-22.
[7]
G. Del Coso, C. Del Canizo, A. Luque. J. Electrochem. Soc. 2008, 155(6): D485-D491.
[8]
H. Habuka, Y. Aoyama, S. Akiyama, S. Akiyama, T. Otsuka, W.F. Qu, M. Shimada, J. Cryst. Growth 1999, 207(1): 77-86.
[9]
H. Habuka, M. Katayama, M. Shimada M, K. Okuyama, J. Cryst. Growth 1997, 182(3): 352-362.
[10]
H. Ni, S. Lu, C. Chen. J. Cryst. Growth 2014, 404: 89-99.
[11]
A. Müller, M. Ghosh, R. Sonnenschein, Mater. Sci. Eng. B-Adv 2006, 134(2): 257-262.
[12]
G. Del Coso, C. Del Canizo, Luque A, Sol. Energ. Mat. Sol. C. 2011, 95(4): 1042-1049.
[13]
D. Cai, L.L. Zheng, Y. Wan, A.V. Hariharanb, M. Chandrab, J. Cryst. Growth 2003, 250(1): 41-49.
[14]
G. Del Coso, C. Del Canizo, A. Luque, J. Electrochem. Soc. 2008, 155(6): D485-D491.
[15]
M. Fang, Y.Y. Xiong, X.Z. Yuan, Y.W. Liu. Energy Procedia, 2014, 61: 1987-1991.
[16]
D. Zhang, N. Xiao, K.T. Mahbubani, Chem. Eng. Sci. 2015, 130:68-78.
[17]
S. Kasiri, S. Abdulsalam, A. Ulrichb, V. Prasad, Chem. Eng. Sci. 2015, 127:31-39.
[18]
S.K. Rout, B.K. Choudhury, R.K. Sahoo, Cryogenics, 2014, 62: 71-83.
[19]
P. Yang, Y.W. Liu, G.Y. Zhong, Appl. Therm. Eng. 2016, 103: 1004-1013.
[20]
A. Balakrishna, J.M. Chacin, P.B. Comita, US DOE, MA. 1998: 1-6.
[21]
H. Van Santen, C.R. Kleijn, H.E.A. J. Cryst. Growth 2000, 212(1): 299-310.
[22]
S. Kommu, B. Khomami, Ind. Eng. Chem. Res. 2002, 41(4): 732-743.
[23]
P. Zhang, W. Weiwen, G. Cheng, Chinese J. Chem. Eng. 2011, 19(1): 1-9.
[24]
H. Ni, S. Lu, C. Chen. J. Cryst. Growth 2014, 404: 89-99.
[25]
Z. Nie, P.A. Ramachandran, Y. Hou, Int. J. Heat Mass Tran. 2018, 117: 1083-1098.
[26]
I.R. Gran, B.F. Magnussen, Combust. Sci. Technol. 1996, 119(1-6): 191-217.
[27]
B.F. Magnussen, B.W. Hjertager, St. Louis, USA. 1981, 198(1).
[28]
H. Habuka, T. Nagoya, M. Mayusumi, J. Cryst. Growth 1996, 169(1): 61-72.
Highlights
A fast and effective method is performed to optimize the silicon CVD system.
The polynomials for predicting silicon deposition rate is obtained.
Analysis of variance is conducted to identify each factor’s significance.
The influence of operating parameters on silicon deposition rate is analyzed.
S0022-0248(18)30083-6 https://doi.org/10.1016/j.jcrysgro.2018.02.030 CRYS 24499
To appear in:
Journal of Crystal Growth
Received Date: Revised Date: Accepted Date:
23 November 2017 12 February 2018 19 February 2018
Please cite this article as: L-s. An, C-j. Liu, Y-w. Liu, Optimization of Operating Parameters in Polysilicon Chemical Vapor Deposition Reactor with Response Surface Methodology, Journal of Crystal Growth (2018), doi: https:// doi.org/10.1016/j.jcrysgro.2018.02.030
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Optimization of Operating Parameters in Polysilicon Chemical Vapor Deposition Reactor with Response Surface Methodology Li-sha An, Chun-jiao Liu, Ying-wen Liu* Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049 P R China
Abstract: In the polysilicon chemical vapor deposition reactor, the operating parameters are complex to affect the polysilicon’s output. Therefore, it is very important to address the coupling problem of multiple parameters and solve the optimization in a computationally efficient manner. Here, we adopted Response Surface Methodology (RSM) to analyze the complex coupling effects of different operating parameters on silicon deposition rate (R) and further achieve effective optimization of the silicon CVD system. Based on finite numerical experiments, an accurate RSM regression model is obtained and applied to predict the R with different operating parameters, including temperature (T), pressure (P), inlet velocity (V), and inlet mole fraction of H2 (M). The analysis of variance is conducted to describe the rationality of regression model and examine the statistical significance of each factor. Consequently, the optimum combination of operating parameters for the silicon CVD reactor is: T=1400 K, P=3.82 atm, V=3.41 m/s, M=0.91. The validation tests and optimum solution show that the results are in good agreement with those from CFD model and the deviations of the predicted values are less than 4.19 %. This work provides a theoretical guidance to operate the polysilicon CVD process. Keywords: Polysilicon; Optimal design; Response Surface Methodology; Numerical method *Corresponding author. Tel.: (86)13087588436 Email address: [email protected] 1
1. Introduction Minimum greenhouse gas emission, sustainability, economy, and easy storage of solar energy are desired properties to form a potential solution to the pending global energy and environmental crisis [1-3]. The dramatic rise of photovoltaic application opens up prospects for the raw silicon material production [4]. It is reported that the current annual production capacity of polysilicon reaches 350 kiloton, in contrast with the 26 kiloton per year in the past [5]. The cost of terrestrial solar cells can be reduced substantially, if solar cells are fabricated from solar grade silicon, which is mainly produced by chemical vapor deposition (CVD) process. Therefore, large-scale production, high quality, and low cost are acknowledged as the most important indexes in practical silicon (Si) chemical vapor deposition (CVD) system. Furthermore, the crucial technique and thermal economic evaluation of the overall CVD process are inherent of the gaseous precursor reduction to solid silicon, which widely occurs in Siemens process [6]. The formation of solid silicon during Siemens process occurs by three steps. Firstly, volatile chlorosilanes are synthesized from metallurgical grade silicon. Then, they are purified by distillation, mixed with hydrogen, and finally fed into the CVD reactor to deposit solid silicon on high-purity silicon rods at a high temperature about 1398 K [7-9]. A variety of by-products are generated in the deposition process, such as HCl, SiCl2, SiH2Cl2, SiHCl, SiCl4, SiH2, Si2H2Cl4, and SiCl6, which will be returned to the recycled system and used to reproduce Trichlorosilane (TCS). The CVD reaction will be ceased till a certain size of silicon rod and the silicon rods will be removed off the reactor after cooling down [11]. 2
Purities in the range of 99.99 % are within reach, but the polysilicon deposition rate is so slow that the deposition process always takes 60-80 hours when the silicon rod grows to certain diameter (150-200 mm) [10]. Furthermore, the intricate situation results in limited conversions of reactants and huge energy dissipation, especially when the heating is performed by electricity [12]. This attracts more attention on raising the silicon deposition rate, thus increasing the efficiency of silicon production and reducing the energy loss of Si CVD process. Examples include changing the geometric structure and controlling operating parameters to obtain the better production conditions, which will be of great significance for the improvement of Si CVD system. By analyzing and optimizing the operating parameters, such as reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M), many researchers have recently succeeded in enhancing the silicon CVD performance. D Cai et al. [13] investigated the silicon production process in a horizontal tubular CVD reactor through computational and experimental methods. They revealed the correlation between the polysilicon rate and the process flow rate, mass fraction of the silane gas in hydrogen, and the temperature of the substrate silicon tube. G. del Coso et al. [14] discussed the effects of gas velocity, the mixture of gas composition, the reactor pressure, and the surface temperature on the silicon growth rate, deposition efficiency, and energy loss. They also provided the favorable range of these operating parameters for silicon deposition performance. In addition to operating parameters, Fang et al. [15] discussed the influence of deposition time on the silicon growth based on numerical models. Historically, research on polycrystalline silicon was preceded by investigating single factor without consideration of interaction between operating parameters. 3
This implies that using single-factor analysis is very difficult to obtain accurate optimal solution, because the coupling interaction between the parameters is too complex and the huge time has been cost to finish massive amounts of experiments. So, it is very important to solve the coupling problem of many parameters and realize the optimization as quickly and effectively as possible. Then the interest in silicon CVD system was renewed with further optimization by seeking a fast and effective approach, which can consider the coupled impact between different parameters and have the characteristics of statistical analysis and accurate prediction. As one of the most efficient methods, the Response Surface Methodology (RSM) introduced by Box and Wilson is widely adopted to investigate the coupled impact between different parameters, to derive correlations of multiple factors and to optimize the system performance [16-19]. It has been widely applied in numerous manufacturing fields for the design, development, and formulation of new products, as well as improving existing product designs. Based on this method, the regression equation of factors and responses can be obtained. Furthermore, intuitive demonstration of the relationship between multiple factors with tiny amounts of experimental data is another advantage. Therefore, RSM can be applied to solve the complicated problem of coupling interaction of CVD system and further improve its performance. This article establishes a two-dimensional model of heat and mass transfer coupling with complex chemical reaction and further verifies the model by comparing with the results in the literature. Based on RSM, 30 tests are conducted, and the regression formula is obtained to accurately predict the average silicon deposition rate under any combination of different 4
parameters. The variance analysis (ANOVA) is accomplished to identify the goodness of regression formula and analyze influential weight of each parameter quantitatively. The influence of parameters, involving reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M) on the average silicon deposition rate has been presented in detail by using 2D contours and 3D surface plot. Finally, the superior average silicon deposition rate is obtained, which are help for the optimization of actual silicon production.
2. Model 2.1 Geometric model The silicon chemical vapor deposition reactions carry out in a two-dimensional plane plates as illustrated in Fig. 1. It consists of a deposition surface, an inlet, an outlet, and an adiabatic wall in the model. The length of the deposition surface is 1.5 m and the height of the inlet and outlet is 0.1 m. The mixtures of trichlorosilane (SiHCl3) and hydrogen (H2) enters the reactor at the left inlet. Then the gases react in the reactor depositing the desired solid silicon on the surface of the substrate. Finally, the reactant and product gases leave the reactors through the right outlet.
Fig. 1. Computational domain.
2.2 Governing equations and boundary conditions
5
Focusing on the deposition process that the polysilicon deposition rate in the order of microns per minute is very little. Thus, at each instant in the process, the system can be assumed to be in the steady state. The wall (ad) is set to be adiabatic, where neither mass nor heat can pass through. The deposition surface is defined as a constant temperature. At the inlet, a velocity boundary is imposed, in which the uniform velocity, temperature and species concentration are defined. An outflow condition is assigned at the outlet plane. The governing equations for continuity, momentum, energy and species in the computational domain are expressed as follows: Continuity equation:
ui 0 xi
(1)
where
ui is the velocity in the i direction. is the fluid density which is expressed as:
1
(2)
mi' i'
i'
Momentum equation: p ij ui u j gi Fi x j xi x j
(3)
Energy equation:
T (ui ( E p)) (keff h ' J ' u j ( ij )eff Sh xi xi xi j' j j
(4)
where keff is the effective thermal conductivity coefficient (including the molecular and turbulent thermal conductivity). The effect of heat transfer by radiation on the substrate 6
temperature is ignored. J j ' is the diffusion flux of component
j ' expressed as equation (5).
S h represents the consumption and production of heat due to the chemical reactions. J i' Di' ,m i Sct Di' ,m
mi' xi
(5)
1 X i' j , j i '
'
'
(6)
X j' / Di' j '
[T 3 ( Di' , j ' 0.0188
1 1 12 )] M i' M j '
(7)
pop i2' , j ' D
T ( / k )i' j'
TD*
1 2
(8)
i j ( i j ) ' '
Sct
'
(9)
'
Dt t
(10)
where Di' ,m is the diffusion flux of species i and
Sct is the turbulent Schmidt number.
Species transport equation:
( Yi ) J i Ri Si
(11)
N
Ri M i Rir
(12)
i 1
where Ri is the net rate of production of species i by chemical reaction, Mi is the molecular weight of species i, and Rir is the Arrhenius molar rate of reaction/destruction of species i in reaction r.
7
2.3 Gaseous and surface reaction mechanisms The detailed chemical reaction mechanism including Pauline Ho’s detailed [20] gas phase reaction and a set of semi-empirical surface reactions has been taken into account to study the characteristics of silicon deposition process. The equations of reactions and the kinetic parameters are summarized in Table 1. Table 1 The gas phase and surface reaction mechanism of trichlorosilane and hydrogen system. No.
Reaction
A
E （cal•mol-1）
G1 G2 G3 G4 G5 G6 G7 G8 G9 S1 S2 S3 S4 S5 S6 S7 S8 S9 S10
SiHCl3=HCl+SiCl2 SiH2Cl2=SiCl2+H2 SiH2Cl2=SiHCl+HCl H2ClSiSiCl3=SiCl4+SiH2 H2ClSiSiCl3=SiHCl3+SiHCl H2ClSiSiCl3=SiH2Cl2+SiCl2 Si2HCl5=SiCl4+SiHCl Si2HCl5=SiHCl3+SiCl2 Si2Cl6=SiCl4+SiCl2 SiHCl3+4Si(s)→Si(b)+SiH(s)+3SiCl(s) SiH2Cl2+4Si(s)→Si(b)+2SiH(s)+2SiCl(s) SiCl4+4Si(s)→Si(b)+4SiCl(s) SiCl2+2Si(s)→Si(b)+2SiCl(s) 2SiCl(s)+Si(b)→SiCl2+2Si(s) H2+2Si(s)→2SiH(s) 2SiH(s)→H2+2Si(s) HCl+2Si(s)→SiH(s)+SiCl(s) SiH(s)+SiCl(s)→2Si(s)+HCl SiHCl+2Si(s)→Si(b)+SiH(s)+SiCl(s)
3.162e14 3.162e13 6.918e14 1.585e13 3.162e13 6.310e13 5.012e13 7.943e13 1.585e14 4.1e-4 1.2e-3 3.0e-5 2.2 4.758e22 0.1 8.606e23 0.07 4.0e25 0.06
72900 69300 75800 55500 49800 44300 52300 45900 48800 -3800 -3800 -3800 15000 67000 17300 57100 5000 71500 5000
A is pre-exponential factor. E is the reaction activation energy.
2.4 Numerical scheme and grid analysis The preprocessor GAMBIT 2.4.6 is used to generate the computational meshes. When the cell size is reduced further, the average deposition rate of silicon changes less than 0.87 % indicating that the numerical results are independent of the grid size with the grid number of 8
45000. FLUENT 14.0 program is devoted to process the present simulation, combined with CHEMKIN introducing detailed chemical reaction mechanism to FLUENT. Reynold numbers in different CVD reactors are always different. In our plane reactor, the Reynold number range from 2500 to 15000, which means that the turbulence flow exists in two-dimensional model. otherwise, in view of the fact that the complex chemical reaction under high temperature happens in the reactor and the concentration of the substrate changes too largely, the turbulence intensity in the CVD system will be increased. Thus, the realizable k-ε turbulent model was adopted to simulate the turbulent transport [21-25]. The reaction model is the eddy-dissipation concept (EDC) model which is an extension of the eddy-dissipation model including detailed chemical mechanisms in turbulent flows [26, 27].
2.5 Validation of model The heat and mass transfer model coupling with gas and surface reaction mechanisms was validated by comparing with the experimental data from the literature [28]. Referring to Habuka’s research, chemical vapor deposition reactions of polysilicon took place in a horizontal reactor. The gas mixture (TCS, H2) was introduced into the reactor under atmospheric pressure with an inlet velocity of 0.67 m/s. The chemical vapor deposition reactions occur on the surface of a disk substrate held horizontally on a fixed susceptor. The substrate was heated to 1398 K using an infrared furnace through a quartz wall reactor. Fig. 2. shows the comparison of experimental and simulated silicon deposition rate under different inlet mole fractions of hydrogen. It can be seen that the simulated results agree well with the experimental values in literature [28]. So this model can accurately describe the polysilicon CVD process. 9
Silicon deposition rate /μm·min
-1
7 6
Experimental datas Simulation result
5 4 3 2 1 0 0.97
0.98
0.99
1.00
Hydrogen mole fraction
Fig. 2. Deposition rate variation of silicon along with the change of hydrogen mole fraction.
3. Response surface methodology 3.1 Experimental design This article aims at developing fast and effective evaluation procedure for investigating the influence of operating parameters (reaction temperature, pressure, inlet velocity, and inlet mole fraction of H2) on the polysilicon deposition rate, thus achieving the better silicon deposition rate. RSM is employed to build accurate relationship between various parameters and responses through multi-factor nonlinear regression. Generally, a second-order model is adopted and given as follows: k
k 1 k
k
y 0 i xi ij xi x j ii xi2 i 1
i 1 j 2
(13)
i 1
where y is the response of the system (the average deposition rate of Si), xi and xj are coded as independent variables; k is the number of variables;
0 , i , ii and ij
(i=0, 1, 2, … k;
j=0, 1, 2, … k) are the regression coefficients for the intercept, linear, quadratic and interaction terms respectively; and
is the statistical error. The flow chart of the evaluation
procedure of this approach is illustrated in Fig. 3. 10
Start
Objective definition Experimental design BBD Fluent calculation
More design points
ANOVA Quality of fit (R2)
No
Yes Response surface analysis
End
Fig. 3. Flow chart of the evaluation procedure based on RSM.
Firstly, in the optimization procedure, the average silicon deposition rate is selected as the response, and four operating parameters including: reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M) are chosen as the variable factors. According to actual reaction furnace, the effective domains of four factors are determined. Next the experimental matrix shown in Table 2 is designed and generated by the Box-Behnken Design (BBD). Table 2 Ranges and levels of independent variables in BBD. Levels Variables -1
0
1
A, T (K)
1250
1350
1450
B, P (atm)
1
3
5
C, V (m/s)
0.5
2.5
4.5
D, M
0.80
0.88
0.96
According to the design matrix arranged by BBD, the model based on ANSYS 14.0 mentioned above is performed to obtain the corresponding responses. The response results under different combinations of four operating parameters in silicon chemical vapor 11
deposition system are tabulated as shown in Table 3. The first column of the table shows the run number of experiments. The next four columns represent the actual conditions of runs and the last column represents the results of experiments.
Table 3 Design of experimental matrix and its responses. Run
T(K)
P (atm)
V (m/s)
M
Deposition rate (μm/min)
1
1350
3
2.5
0.88
10.43
2
1450
3
4.5
0.88
41.76
3
1250
5
2.5
0.88
2.55
4
1350
5
2.5
0.96
9.87
5
1250
3
2.5
0.80
1.77
6
1350
3
2.5
0.88
10.43
7
1450
1
2.5
0.88
17.47
8
1350
1
2.5
0.80
5.84
9
1250
3
4.5
0.88
2.58
10
1350
3
2.5
0.88
10.43
11
1350
1
4.5
0.88
7.22
12
1350
5
4.5
0.88
16.54
13
1350
1
2.5
0.96
4.81
14
1250
3
0.5
0.88
1.42
15
1350
3
0.5
0.96
4.00
16
1450
3
2.5
0.96
20.38
17
1450
3
0.5
0.88
14.05
18
1250
3
2.5
0.96
1.81
19
1350
3
2.5
0.88
10.43
20
1250
1
2.5
0.88
1.15
21
1450
5
2.5
0.88
40.80
22
1450
3
2.5
0.80
35.75
23
1350
3
0.5
0.80
7.94
24
1350
5
0.5
0.88
7.78
25
1350
3
2.5
0.88
10.43
26
1350
1
0.5
0.88
4.16
27
1350
3
4.5
0.96
10.04
28
1350
3
4.5
0.80
12.79
29
1350
5
2.5
0.80
14.15
12
3.2 Analysis of variance (ANOVA) Analysis of variance is carried out to evaluate the significance of the regression models and influence weight of factors. It discusses the variance of variable and obtains the evaluation index of variance. The ANOVA results of the quadratic model for average silicon deposition rate are shown in Table 4. It suggests that the model was significant, as it is evident from the model’s F-value with a low probability value (P<0.0001). Values of “prob>F” less than 0.05 indicates that model terms are significant. In this case, A, B, C, D, AB, AC, AD and A2 are statistically significant for the average silicon deposition rate. Other interaction term like BC with larger p-values can be considered to have minor influence on the average silicon deposition rate. The goodness of fitting formula is further checked by the correlation coefficient (R2) between the experimental model and predicted values of the response variable. A high R2 value of 0.9730 implies that the regression model is statistically significant and only 2.7 % of the total variations is not explained by the model. The predicted correlation coefficient (pred. R2=0.9536) also shows good agreement with the adjusted correlation coefficient (adj. R2=0.9459). In order to obtain the accurate regression model between factors and responses, the significant and insignificant terms are all considered. The final mathematical formula in terms of actual factors is shown as follows: Average Si deposition rate= 371.76439-1.10227*A-30.39349 *B-44.24245*C+911.04461*D+ 0.027420*A*B+0.03318*A*C0.48168*A*D-0.35643*B*C5.07413*B*D-1.84225*C*D+ 55.5310E-004*A2-0.15365*B20.21824*C2-158.31953*D2
13
(14)
Table 4 ANOVA for average silicon deposition rate. source
Sum of squares
df
Mean square
F value
p-value prob>f
model
3226.51
14
230.46
35.99
< 0.0001
Significant
A-T
2104.96
1
2104.96
328.70
< 0.0001
Significant
B-P
217.21
1
217.21
33.92
< 0.0001
Significant
C-V
221.84
1
221.84
34.64
< 0.0001
Significant
D-M
62.34
1
62.34
9.73
0.0075
Significant
AB
120.30
1
120.30
18.78
0.0007
Significant
AC
176.15
1
176.15
27.51
0.0001
Significant
AD
59.40
1
59.40
9.27
0.0087
Significant
BC
8.13
1
8.13
1.27
0.2788
BD
2.64
1
2.64
0.41
0.5315
CD
0.35
1
0.35
0.054
0.8192
A2
198.44
1
198.44
30.99
< 0.0001
B2
2.45
1
2.45
0.38
0.5461
2
C
4.94
1
4.94
0.77
0.3945
D2
6.66
1
6.66
1.04
0.3251
Residual
89.65
14
6.40
Lack of fit
89.65
10
8.97
Pure error
0.000
4
0.000
Cor total
3316.16
28
R-squared=0.9730
Significant
Adj R-squared=0.9459
3.3 Analysis of regression formula The model of least squares is typically used to estimate the regression coefficients in a multiple linear regression model. The residuals from the least squares fit play an important role in judging model adequacy. The normal probability plot obtained from RSM for the average silicon deposition rate is illustrated in Fig. 4. It can be seen that the residuals plot approximately along a straight line, indicating that the errors are normally distributed. Besides, almost all the internally studentized residuals are in small range near to zero, 14
meaning that the transformation of the response may provide a better analysis.
Fig. 4. Normal probability plot of residuals for average silicon deposition rate.
The relationship between the average silicon deposition rate obtained from our model and the value predicted by regression model is illustrated in Fig. 5. In this figure, the design points are distributed evenly on the diagonal line, and the error between the predicted values from regression formula and the actual CFD results is very little. It signifies that the regression formula can be used to predict the average silicon deposition rate accurately.
Fig. 5. Regression model predictions versus our model’s results for average silicon deposition rate.
4 Results and discussion From the mathematical formula, it is obvious that there exists the coupling effect of operating parameters on the silicon deposition rate. The 2-D contour plots and 3-D surface 15
plots present the interaction effects between any two factors involving reaction temperature (T), pressure (P), inlet velocity (V) and inlet mole fraction of H2 (M). Finally, the relationship between the parameters and silicon deposition rate had been revealed. Four parameters had been optimized to improve the performance of silicon chemical vapor deposition system. 4.1 Influence of temperature and pressure The combined effect of temperature and pressure on the average silicon deposition rate, at inlet velocity of 2.5 m/s and the inlet H2 mole fraction of 0.88, is shown in Fig. 6 (a) and Fig. 6 (b), respectively. From Fig. 6, we can note that compared with the pressure, the reaction temperature plays the major role in affecting the silicon deposition rate as evidenced by the larger lateral gradient in 2D contour plot. It is apparent that the average silicon deposition rate increases as the reaction temperature rises. This is because the enlarging reaction temperature promotes the silicon CVD reaction, which corresponds to the property of endothermic reactions. In addition, higher temperature further results in better crystalline properties of silicon and more metallic silicon surfaces. While, it is important to note that, under the operating condition of high temperature, the silicon will be easily polluted by equipment material and can further result in an increase of phosphorus and boron impurities because of the enhanced chemical activity of silicon. Meanwhile, the high reaction temperature also enhances the silicon corrosion reactions, seriously inhibiting the Si deposition reactions and reducing the quality of silicon. Therefore, the reaction temperature of actual reduction furnace is generally controlled at the range from 1353 to 1423 K. From Fig. 6, it also can be observed that the average silicon deposition rate increases slowly with the increase of pressure. It reveals that the silicon deposition rate is little affected 16
by pressure. The reason is that the raising pressure increases the crash chance for SiHCl3 and H2, and then improves the silicon CVD reactions. however, with low operating reaction temperature, chemical reacting rate is so low that only an amount of feed gas can effectively participate in the complex reaction. though the feed gas becomes more and more with the increase of pressure, the reaction rate is mainly determined by the reaction temperature and chemical mechanism, and the effect of the pressure is less than that of them. Otherwise, considering that the rising pressure will aggravate the equipment load, the pressure from 2 to 4 atm is conductive to the silicon CVD system in actual reduction furnace.
(a) 2-D contour plot (b) 3-D surface plot Fig. 6. Effect of temperature and Pressure on the average silicon deposition rate.
4.2 Influence of inlet velocity and inlet mole fraction of H2 The combined effect of inlet velocity and inlet mole fraction of H 2 on the average silicon deposition rate, at the temperature of 1350 K and the pressure of 3 atm, is shown in Fig. 7 (a) and Fig. 7 (b), respectively. Firstly, it can be observed that the silicon deposition rate increases with inlet velocity rising due to the sufficient reactant induced by the high velocity. But under the mole fraction of 0.8, the relative increment of silicon deposition rate is only 60 % while inlet velocity 17
increases from 1 m/s to 3 m/s. the relative increment of the silicon deposition rate gradually decreases along with the increase of the mole fraction. The reason is that as the inlet velocity continues to increase, the severely reduced reaction time becomes the dominant factor which limits the conversion rate of SiHCl3 to solid silicon, resulting in wasting a lot of reaction gases and electric energy. So, appropriate inlet velocity should be chosen to reduce costs and make rational use of resources. Secondly, the silicon deposition rate reduces as the inlet mole fraction of H 2 increases. The reason is that the rising concentration of H2 will dilute the concentration of SiHCl3 and then reduce the crash probability of SiHCl3 molecule with silicon rod surfaces. In actual reduction furnace, the inlet mole fraction of H2 from 0.89 to 0.95 is selected to reduce cost and save energy.
(a) 2-D contour plot (b) 3-D surface plot Fig. 7. Effect of inlet velocity and inlet mole fraction of H2 on the average silicon rate.
4.3 Influence of reaction temperature and inlet velocity The combined effect of reaction temperature and inlet velocity on the average silicon deposition rate, at the pressure of 3 atm and the inlet H2 mole fraction of 0.88, is shown in Fig. 8 (a) and Fig. 8 (b) in terms of 2-D contour plot and 3-D surface plot, respectively. In this case, the average silicon deposition rate is below 5 μm/min when the temperature is lower 18
than 1300 K, where the reaction temperature is the dominant factor that affect the silicon CVD reaction. High reaction temperature and inlet velocity are conducive to the achieve better silicon deposition rate.
(a) 2-D contour plot (b) 3-D surface plot Fig. 8. Effect of temperature and inlet velocity on the average silicon rate.
4.4 Influence of temperature and inlet mole fraction of H2 The combined effect of temperature and inlet mole fraction of H2 on the average silicon deposition rate, at the pressure of 3 atm and inlet velocity of 2.5 m/s, is shown in Fig. 9 (a) and Fig. 9 (b), respectively. The average silicon deposition rate is less than 5 μm/min in the area where the temperature below 1300 K. Under this condition, the reaction temperature is the main factor corresponding to the analysis above. The silicon deposition rate can be improved by higher temperature and lower inlet mole fraction of H2. It is apparent that the temperature has a more significant effect on the average silicon deposition rate compared with the inlet mole faction of H2.
19
(a) 2-D contour plot (b) 3-D surface plot Fig. 9. Effect of temperature and inlet mole fraction of H2 on the average silicon rate.
4.5 Influence of pressure and inlet velocity The combined effect of pressure and inlet velocity on the average silicon deposition rate, at the temperature of 1350 K and the inlet H2 mole fraction of 0.88, is shown in Fig. 10 (a) and Fig. 10 (b) in terms of 2-D contour plot and 3-D surface plot, respectively. The average silicon deposition rate can be enhanced by increasing the inlet velocity and pressure in the reactor. Furthermore, the different influencing rule of them shows that the inlet velocity has more significant effect on the silicon deposition rate by comparing with the effect of pressure.
(a) 2-D contour plot (b) 3-D surface plot Fig. 10. Effect of pressure and inlet velocity on the average silicon rate.
4.6 Influence of pressure and mole fraction of H2 The combined effect of pressure and mole fraction of H2 on the average silicon 20
deposition rate, at the temperature of 1350 K and inlet velocity of 2.5 m/s, is shown in Fig. 11 (a) and Fig. 11 (b), respectively. The heat and mass transfer characteristics of the silicon chemical vapor deposition process can be enhanced by increasing the pressure and reducing the mole fraction of H2, thus improving silicon deposition rate.
(a) 2-D contour plot (b) 3-D surface plot Fig. 11. Effect of pressure and inlet mole fraction on the average silicon rate.
In conclusion, appropriate design of operating parameters should be selected to control the silicon chemical vapor deposition reactions. Among operating parameters, the influence of reaction temperature on the average silicon deposition rate is most significant. Next is the influence of inlet velocity and the pressure. The influence of H 2 mole fraction is least. In actual operation of silicon chemical vapor deposition system, in order to satisfy the energy saving, high quality of polysilicon and high efficiency, the operating manual should be designed in this sequence. According to the 2-D contour plot and 3-D surface response, higher temperature, higher pressure, larger inlet velocity and little inlet mole fraction of H 2 are beneficial to improve the performance of silicon chemical vapor deposition system.
4.7 Confirmation tests and optimization In this chapter, four experimental combinations including the optimal one are performed 21
to further prove the adequacy of the regression formula obtained by RSM. The average silicon deposition rates and comparison of four tests obtained by regression model and CFD model are shown in Table 5. From this table, the maximum deviation between numerical solutions and predicted values of the silicon deposition rate is 4.19 %, which is acceptable for engineering application. It indicates that the regression model can accurately predict the average silicon deposition rate within the range of given parameters of the silicon CVD process. In addition, the confirmation test 1 is the optimal operating parameter combination, whose average silicon deposition rate increase by 12 percent compared with the maximum value in numerical experiments matrix. Further analyzing and comparing the operating parameters between them, the reaction temperature and inlet velocity of optimal combination are less than these of best combination of numerical experiments matrix. It means that after adopting optimal combination, it is significative to decrease the energy consumption and ensure safe operation of silicon chemical vapor deposition system.
Table 5 Validation results for the average silicon deposition rate. Design parameters
Average Si deposition rate
Confirmation runs
Optimal
RSM
Fluent
/(m/min)
/(m/min)
0.91
45.86
46.78
2.2%
2.5
0.88
18.69
19.39
4.19%
2
3.5
0.9
9.24
9.42
2.00%
1
1.5
0.86
2.93
2.34
-2.00%
T/K
P/atm V/m/s
1
1400
3.82
3.41
2
1400
3
3
1350
4
1300
M
Error *
5. Conclusion In this study, a two-dimensional numerical model was developed to simulate the 22
characteristics of silicon chemical vapor deposition. The numerical model was validated by experimental data. Response surface methodology was selected to establish mathematical models of the average silicon deposition rate as a function of the operating parameters involving temperature, pressure, inlet velocity, and inlet mole fraction of H 2. The developed models were verified by the residual distribution between predicted and actual results. The operating parameters effects and their interactions on the average silicon deposition rate can be explained in detailed by ANOVA. It was found that four operating parameters have different significance on the average silicon deposition rate. Furthermore, confirmation tests were carried out to verify the accuracy of the regression model and the optimized parameters was determined. Results showed a good agreement between numerical and predicted results. This study could provide an efficient and accuracy method to guide and optimize the productive process of CVD.
Acknowledgements: This work was supported by National Key R&D Program of China (Grant No. 2016YFE0204200) and the National Science Fund from the National Natural Science Foundation of China (No. 51576150)
Nomenclature： Abbreviations ANOVA
Analysis of variance
BBD
Box-Behnken Design 23
CVD
Chemical vapor deposition
CFD
Computational Fluid Dynamics
EDC
Eddy-dissipation concept
PV
Photovaltaic
RSM
Response Surface Methodology
TCS
Trichlorosilane
Symbols A
Pre-exponential factor.
Di'm
Diffusion flux of species i
E
Reaction activation energy (cal/mol)
J j'
The diffusion flux of component
kef
Effective thermal conductivity coefficient
k
Number of variable
m
Quality (kg)
M
Molar mass (g/mol)
M
Inlet mole fraction of H2
P
Pressure (atm)
qw
Heat flux (w/g)
R2
Correlation coefficient
Ri
The net rate of production of species i by chemical reaction
Rir
The Arrhenius molar rate of reaction/destruction of species i in reaction r
Sh
The consumption and production of heat due to the chemical reactions
j'
24
Sct
Turbulent Schmidt number
T
Temperature (K)
u
Fluid velocity (m/s)
V
Inlet velocity
xi,xj
Independent variable
y
Response
Greek
The temperature-dependent collision integral
Tangential force (N/m2)
The average collision diameter
Statistical error
Fluid density (kg/m3)
0
Regression coefficient for intercept
i
Regression coefficient for linear
ii
Regression coefficient for quadratic
ij
Regression coefficient for interaction terms
Reference [1]
A. Müller, M. Ghosh, Sonnenschein R, P. Woditsch, Mater. Sci. Eng. B-Adv., 2006, 134(2): 257-262.
[2]
A.F.B. Braga, S.P. Moreira, P.R. Zampieri, J.M.G. Bacchin, P.R. Mei, Sol. Energ. Mat. Sol. C. 2008, 92(4): 418-424.
[3]
P. Woditsch, W. Koch. Sol. Energ. Mat. Sol. C. 2002, 72(1): 11-26.
[4]
W.O. Filtvedt, M. Javidi, A. Holt, Sol. Energ. Mat. Sol. C. 94 (12) (2010) 1980-1995. 25
[5]
S. Ravasio, M. Masi, C. Cavallotti, J. Phys. Chem. A 2013, 117(25): 5221-5231.
[6]
B. Hazeltine, C. Fero, W. Qin. Photovoltaics World, 2010, 3: 16-22.
[7]
G. Del Coso, C. Del Canizo, A. Luque. J. Electrochem. Soc. 2008, 155(6): D485-D491.
[8]
H. Habuka, Y. Aoyama, S. Akiyama, S. Akiyama, T. Otsuka, W.F. Qu, M. Shimada, J. Cryst. Growth 1999, 207(1): 77-86.
[9]
H. Habuka, M. Katayama, M. Shimada M, K. Okuyama, J. Cryst. Growth 1997, 182(3): 352-362.
[10]
H. Ni, S. Lu, C. Chen. J. Cryst. Growth 2014, 404: 89-99.
[11]
A. Müller, M. Ghosh, R. Sonnenschein, Mater. Sci. Eng. B-Adv 2006, 134(2): 257-262.
[12]
G. Del Coso, C. Del Canizo, Luque A, Sol. Energ. Mat. Sol. C. 2011, 95(4): 1042-1049.
[13]
D. Cai, L.L. Zheng, Y. Wan, A.V. Hariharanb, M. Chandrab, J. Cryst. Growth 2003, 250(1): 41-49.
[14]
G. Del Coso, C. Del Canizo, A. Luque, J. Electrochem. Soc. 2008, 155(6): D485-D491.
[15]
M. Fang, Y.Y. Xiong, X.Z. Yuan, Y.W. Liu. Energy Procedia, 2014, 61: 1987-1991.
[16]
D. Zhang, N. Xiao, K.T. Mahbubani, Chem. Eng. Sci. 2015, 130:68-78.
[17]
S. Kasiri, S. Abdulsalam, A. Ulrichb, V. Prasad, Chem. Eng. Sci. 2015, 127:31-39.
[18]
S.K. Rout, B.K. Choudhury, R.K. Sahoo, Cryogenics, 2014, 62: 71-83.
[19]
P. Yang, Y.W. Liu, G.Y. Zhong, Appl. Therm. Eng. 2016, 103: 1004-1013.
[20]
A. Balakrishna, J.M. Chacin, P.B. Comita, US DOE, MA. 1998: 1-6.
[21]
H. Van Santen, C.R. Kleijn, H.E.A. J. Cryst. Growth 2000, 212(1): 299-310.
[22]
S. Kommu, B. Khomami, Ind. Eng. Chem. Res. 2002, 41(4): 732-743.
[23]
P. Zhang, W. Weiwen, G. Cheng, Chinese J. Chem. Eng. 2011, 19(1): 1-9.
[24]
H. Ni, S. Lu, C. Chen. J. Cryst. Growth 2014, 404: 89-99.
[25]
Z. Nie, P.A. Ramachandran, Y. Hou, Int. J. Heat Mass Tran. 2018, 117: 1083-1098.
[26]
I.R. Gran, B.F. Magnussen, Combust. Sci. Technol. 1996, 119(1-6): 191-217.
[27]
B.F. Magnussen, B.W. Hjertager, St. Louis, USA. 1981, 198(1).
[28]
H. Habuka, T. Nagoya, M. Mayusumi, J. Cryst. Growth 1996, 169(1): 61-72.
Highlights
A fast and effective method is performed to optimize the silicon CVD system.
The polynomials for predicting silicon deposition rate is obtained.
Analysis of variance is conducted to identify each factor’s significance.
The influence of operating parameters on silicon deposition rate is analyzed.