Abstract
In order to obtain higher quality silicon single crystal, a hybrid strategy for modeling of Czochralski silicon crystal growth and optimizing of process parameters is presented in the paper. The hybrid strategy includes the computational fluid dynamics (CFD) method, neural network of group method of data handing (GMDH), and improved nondominated sorting genetic algorithm II (NSGA-II). The shape variable of solid-liquid interface and the defect evaluation criteria are set to objective functions according to engineering experience and process requirement. The polynomial of the objective functions is produced by GMDH and CFD. Ultimately, an improved elitist strategy and crowding distance NSGA-II is proposed in the paper to obtain the Pareto optimum solution by the objective functions identified by the GMDH. Compared with other optimization algorithms, the improved NSGA-II can increase the lateral diversity and the uniform distribution of the nondominated solutions. Engineering validation proved that the proposed hybrid strategy can effectively solve the complex uncertain multiobjective optimization problem of system and provides a new method for obtaining high-quality crystal growth process parameters.
1. Introduction
The vigorous development of integrated circuit industry promoted the development of silicon crystal in the direction of larger size and higher quality [1]. The process of silicon crystal growth is the key factor to measure the production level of semiconductor materials. The process parameters are usually set from numerical simulation methods or engineering experience. But the process parameters obtained by those methods are often not optimal. Therefore, it is a new direction to study the intelligent optimization method of crystal growth process parameters to improve the growth process and improve the crystal quality.
At present, there are two main research directions for process parameters in the field of silicon single crystal growth. One is the combination of empirical method and experimental method, the other is numerical simulation method based on computational fluid dynamics (CFD). Empirical and experimental methods rely on specific equipment, which requires experienced technicians and strong experimental funds. The production units often use the experimental method in the final verification stage. The experimental data and experience obtained are also confidential. Because of the advantages of fast, high efficiency, low cost, and flexible design, numerical simulation method has been widely used in the research of crystal growth.
In recent years, the study of crystal growth technology has expanded from two-dimensional global model to three-dimensional local model, from simple fluid-thermal coupling to impurity transport, and from fixed boundary to complex boundary conditions involving functional and chemical reactions [2–4]. At present, the researches on optimization of process parameters are still at the level of numerical simulation and mechanism analysis. Liu et al. studied the effect of temperature and external magnetic field on the deformation of solid-liquid interface [5]. Asadi et al. studied the change of solid-liquid interface at different rotating speeds of crystals and crucibles [6]. Voronkov and Falster studied the effect of growth parameters on crystal quality from the factors of crystal defects and proposed a criterion for evaluating the probability of vacancies and self-gap in crystals [7]. Neural networks and intelligent optimization methods have been applied more and more in the field of crystal growth. Qi et al. established a heat and mass transfer model to simulate the crystal growth process and verified the experimental data with neural network and genetic algorithm to obtain the optimal growth parameters [8]. Dang et al. established the heat transfer model and a thermal stress model in single crystalline silicon growth and studied the mapping relationship between optimization variables and targets [9]. The above researches based on numerical simulation can effectively explain the relationship between single process parameter and crystal quality, but these conclusions are limited to the basis of limited calculation results, which is helpless for multiobjective and multiprocess parameters cooperative optimization.
For parameter optimization problems, researchers always expect a clear mathematical model or objective function, in which CFD cannot achieve. However, the combination of CFD, model identification, and multiobjective optimization is a good strategy. This hybrid strategy has been applied in many industrial fields and has attracted more and more attention [10–12]. Wang et al. simulated the turbine model use CFD and provided training samples for BP neural network. The maximum efficiency and the minimum blocking quality were optimized by using NSGA-II with improved congestion distance [13]. Damavandi et al. established the models of maximum heat transfer efficiency and minimum pressure drop function of corrugated fin elliptical heat exchanger by combining CFD with the group method of data handing (GMDH) neural network. The improved nondominated sorting genetic algorithm II (NSGA-II) is used to optimize the structure design, and the optimal structural design scheme is obtained [14]. Jafari et al. modeling and calculation of DI diesel engine used CFD and response surface fitting Used Kriging interpolation method, and optimize the first efficiency law and the second efficiency law used multiobjective optimization algorithm [15]. Because it is difficult to establish the model of optimization objective and process parameters, the application of this hybrid strategy in the field of crystal growth has not been reported yet.
In the paper, the process of 300 mm silicon crystal growth by the Cz method is studied. The solid-liquid interface deformation and value, which represent the crystal quality, are taken as objective functions. Heating temperature , crystal rotating velocity , melt rotating velocity , and pulling speed are the process parameters to be optimized. Firstly, according to design variables, the three-dimensional local model of the Cz crystal growth was solved by CFD technology, and training samples for GMDH were provided. Then, the polynomial models of two objective functions are obtained by the GMDH algorithm. Finally, an improved NSGA-II is proposed to optimize the process parameters and obtain the Pareto frontier, which satisfies the requirements of two objective functions at the same time. This improved NSGA-II includes controlling elite strategy to limit the number of individuals on nondominated frontier, improving the congestion distance, and using dynamic crowding distance (DCD) to improve the Pareto frontier distribution. The effectiveness of the hybrid algorithm is verified by crystal growth experiments.
2. Crystal Growth Model and Numerical Solution
2.1. Mathematical Models
The Cz silicon crystal growth consists by high-temperature melt and crystal. The structure of the Cz silicon crystal is shown in Figure 1. The melt in crucible is kept liquid state by the high temperature. The crystal rises slowly at a certain pulling speed. In order to make the temperature in the melt evenly distributed, the crucible and the crystal rotate in the opposite direction at different speeds. The continuum equation, momentum conservation equation, and energy conservation equation of melt flow and heat transfer calculation are expressed as follows [16]. where , , , , , , , and are temperature, density, pressure, melt viscosity, gravity acceleration, thermal expansion coefficient, thermal capacity, and thermal conductivity, respectively. The melt velocity is mainly related to the rotational speed of crystals and melts. The last term in Equation (2) denotes thermal buoyancy, which is one of the main factors affecting the conservation of momentum.
At the solid-liquid interface, it is assumed that the temperatures of solid and liquid phases are equal, that is, . The energy balance equation at the solid-liquid interface is as follows: where is the normal vector of interface. and are interface moving speed and crystal lifting speed, respectively. is latent heat. is thermal conductivity, and subscripts and indicate crystals and melts, respectively. The deformation position of solid-liquid interface at each moment can be obtained by calculating in Equation (4).
2.2. Numerical Simulation
The Fluent software is used to simulate the Cz silicon single crystal growth 3D local model in the paper, which is an internationally popular commercial CFD software. Convection and diffusion terms are upwind of second order, and the SIMPLE algorithm [5] is used to solve the coupling between velocity and pressure fields. The two-equation realizable turbulence model is selected, which ensures the continuity of turbulence and has a good performance for rotational flow and boundary layer flow.
In order to ensure the independence of the grid and minimize the computational cost, numerical analysis is carried out under three models with different mesh accuracy. The results of objective functions and on different grids are shown in Table 1, when the computational parameters are , , , and . The errors between the results of 232436 mesh models and those of the denser mesh models can be neglected. So, 232436 mesh models satisfy the requirement of independence precision. The meshing of the three-dimensional model of crystal growth and the result of numerical simulation are shown in Figures 2(a) and 2(b), respectively.
(a)
(b)
3. Design Variables and Objective Functions
Four design variables related to the objective function are selected in the paper, which are the heater temperature , the crystal rotating velocity , the melt rotating velocity , and the lifting speed . These variables and their value ranges are listed in Table 2. Many design points can be obtained by changing the value of the design variables, so the corresponding objective function values of these design points can be calculated by CFD.
The two objective functions are the deformation of solid-liquid interface and the defect evaluation criteria . Because of the change of heat transfer from melt to crystal, the convex deformation of solid-liquid interface is produced under the action of thermal stress. The larger the deformation , the greater the probability of growth defects and the thermoelastic stress in the crystal. is an important index for “self-gap” and “vacancy” in silicon wafers. Literature [17] shows that the range of lifting speed is related to the size of crystal diameter within the required range of . The perfect range of values can be ensured when the lifting speed is set between and for 300 mm diameter silicon crystal. The larger the , the larger lifting speed, and the more conducive to shortening the production cycle, saving energy consumption, and reducing costs. However, the increase of lifting speed will increase and reduce the uniformity of crystal cross-section. The two objective functions are conflicting and constitute a multiobjective optimization problem in the process of silicon crystal growth. Therefore, the purpose of this study is to find the optimal solution which satisfies the two objective functions simultaneously.
4. Object Function Modeling by GMDH
Although the approximate solutions of the two objective functions can be obtained by numerical calculation, the explicit mathematical expressions containing the four design variables mentioned above cannot be obtained. However, the accurate model of objective functions is necessary for system control and process parameter optimization. The artificial neural network cannot be restricted by the nonlinear model because it has the ability to approximate any nonlinear mapping by learning. So, the artificial neural network method is used to determine model has obvious advantages. GMDH is a relatively complete and widely used neural network structure [18]. The GMDH algorithm works by describing the model as the structure of a group of neurons. Each pair of neurons in each layer is connected by a quadratic polynomial and new neurons are created in the next layer. During modeling, the input-output mapping can be represented using this method.
The way to solve the problem of model identification is to establish a function that can approximately replace the actual model , so as to minimize the error between the estimated output and the actual model output when the input vector is specified. Therefore, the -estimated values include inputs, and one output can be expressed as follows:
Estimate of the target value of each input vector by the GMDH neural network is as follows:
The GMDH neural network can minimize the difference square between the actual value and the predicted value, which is
The overall relationship between input and output variables can be expressed as a complex discrete Volterra function:
The above equation can be simplified to the composition form of a partial quadratic polynomial containing only two variables, namely
The parameters in the equation are calculated by the regression method.
Obviously, the greater the amount of input and output data used in the calculation, the more accurate and effective the GMDH method will be. In the paper, each set of input and output data is obtained by the CFD calculation. For each set of numerical simulation calculation, especially for three-dimensional model calculation with large number of grids, there is a certain calculation time cost, and the acquisition of these data cannot be infinite. Therefore, the limited amount of data should be reasonably optimized to obtain the maximum information. In this paper, 245 groups of effective CFD results are used as input and output data, of which 196 groups are training samples and 49 groups are used to test neural networks. The training set trains the neural network according to GMDH algorithm, and the test set is used to analyze the performance of the network. The inputs are four design variables , , , and . The outputs are objective functions and . The structure of GMDH neural network with two hidden layers is shown in Figures 3(a) and 3(b). The GMDH polynomial coefficients and equations of functions can be expressed as follows:
(a)
(b)
Similarly, the GMDH polynomial coefficients and equations of the objective function can be expressed as follows:
The comparison results of objective functions and obtained from CFD and GMDH are shown in Figures 4(a) and 4(b). At the same time, the backpropagation artificial neural network (BP-ANN) algorithm is used to calculate the objective function and compare with CFD results. The result comparisons are shown in Figures 5(a) and 5(b).
(a)
(b)
(a)
(b)
In order to evaluate the performance of the model established by GMDH, two statistical parameters are introduced, that is absolute variance fraction and mean absolute percentage error . The formulas were expressed as follows:
The calculated results of the statistical parameters and of the objective function by GMDH and BP-ANN are shown in Table 3. The results show that the GMDH model has high accuracy in predicting CFD data and the prediction accuracy is slightly higher than that of the BP-ANN algorithm.
5. Optimization Results of NSGA-II and Analysis
5.1. Improved NSGA-II
NSGA-II is a universally recognized multiobjective optimization algorithm with extensive applications and excellent effects [19]. It contains many advanced concepts, such as elite strategy, fast nondominant ranking, and diversity maintenance of the Pareto frontier. However, there are still deficiencies in NSGA-II in maintaining transverse diversity and uniform distribution of nondominant solutions [20]. The main reason for these defects is that they simply judge the individual distribution from the crowded distance but ignore the uneven distribution caused by the individual density, which is easy to fall into the local optimal. To overcome the shortcomings of NSGA-II, an improved scheme is proposed in the paper. This scheme is divided into two aspects. Firstly, to ensure the diversity of nondominant individuals by exert control on elite strategies. Secondly, the congestion distance formula is improved, and dynamic congestion distance (DCD) is used to improve the distribution of noninferior solutions.
Controlling elite strategy is to limit the number of individuals in the current best nondominant frontier and maintain the predistribution of the number of individuals in each frontier. A new population is ordered nondominated, which is formed by combining the parent population and the offspring population . The nondominant frontier number of the merged population is named . According to the geometric distribution, the number of individuals allowed at the th frontier in the new population is the largest and the size of the population is given by where is the decrease rate. So, the maximum number of individuals allowed on the first nondominant front is the maximum. The number of individuals allowed for each subsequent front is exponentially lower.
To eliminate the inhomogeneity of individual distribution and improve lateral diversity in the Pareto frontier, Luo et al. proposed a dynamic congestion distance (DCD) method [21]. This method deletes the individual with the lowest DCD value calculated each time and recalculates the DCD value for the remaining individual. Because the calculation of DCD value still depends on the crowding distance, the shortcomings of crowding distance calculation are not avoided. In the paper, a simple improvement of crowding distance is made, to allocate a characteristic coefficient to the crowding distance to balance the contribution of each target to the crowding distance. where is the improved crowding distance for the th individuals. and are corresponding characteristic coefficients and original congestion distances, respectively. and are the maximum and minimum of the objective function for all individuals on the same front. The improved individual dynamic congestion distance value is where is the variance of crowding distance between neighbors of the th individual, which can be expressed as
The population size is assumed as , and the size of combined population of the -generation nondominant set is . If , individuals are removed from the nondominant concentration by DCD strategy. The process of DCD algorithm is briefly described as follows.
Step1: if , go to Step5; otherwise, go to Step2
Step2: calculating the DCD value of individuals in , use formula (15)
Step3: ranking nondominant set based on DCD strategy
Step4: removing the individual of lowest DCD in
Step5: if , stop population maintenance; otherwise, transfer to Step2 and continue execution
5.2. Implementation of the Improved NSGA-II Strategy
The steps of the improved NSGA-II proposed in the paper for the optimization of the above two objective functions in Cz crystal growth are described as follows.
Step1: set initial parameters, such as population size , probability of crossover and mutation, index of crossover and mutation, iterations, and upper and lower bounds of design variables
Step2: generate random initial population in the range of design variables and iteration numbers
Step3: evaluate the objective functions and with each individual of each generation
Step4: creating offspring population by the tournament selection, simulated binary crossover operator, and polynomial mutation
Step5: nondominant sequencing of merged population
Step6: to limit the number of individuals at the best nondominant frontier at present, use control elite strategy and maintain the predistribution of the number of individuals at each frontier use the formula (13). The decrease rate is 0.5
Step7: if the numbers of nondominant set is larger than the numbers of population, then use the DCD strategy to remove individuals from nondominant sets, and then jump to Step4
Step8: the algorithm stopped when the iteration count reaches the maximum iterations; otherwise, the number of iterations increases 1 and then jumps to Step3
5.3. Result Analysis
The ultimate goal of the paper is to obtain the optimum process parameters with the requirements of minimizing the interface deformation and maximizing . However, the objective functions cannot be optimally satisfied at the same time, so the feasible solution is to seek Pareto optimal sets. The optimization research for Cz crystal growth can be expressed in the following forms
In order to prove the validity and correctness of improved NSGA-II in optimizing the process parameters of crystal growth, four multiobjective optimization algorithms are selected to compared with the results of our algorithm, including NSGA-II, multiobjective particle swarm optimization algorithm (MOPSO) [22], second-generation multiobjective genetic algorithm (MOGA-II) [23], and second-generation strong Pareto evolutionary algorithm (SPEA-II) [24]. The Pareto frontier results of five algorithms are shown in Figure 6, and MNSGA-II is our algorithm. NSGA-II and MOPSO fall into local optimum easily and cannot get the optimum Pareto front. MOGA-II and SPEA-II have similar Pareto frontier trends, but there are obvious shortcomings. Although the MOGA-II algorithm maintains a good individual distance, there are some individuals are missed in the frontier. SPEA-II maintains the continuity of the frontier, but the individual distance is not good. The improved NSGA-II effectively avoids the shortcomings of MOGA-II and SPEA-II and obtains a relatively excellent Pareto frontier.
The Pareto frontier obtained by the improved NSGA-II is shown in Figure 7. Five characteristic optimization points, A, B, C, D, and E are selected. It is clear that each optimization point cannot be called absolute optimization for another point, that is to say, from one point to another, one objective function will become better, and the other one will become worse. The corresponding design variables are shown in Table 4. As can be seen in Figure 7, from point A to point E, the interface deformation increases gradually with the increase of . From point A to point B, the increase of interface deformation is very small (about 10%), while the increase of is very large (about 37.5%). Similarly, from point D to point E, the interface deformation increases significantly (about 67.1%), while has little change (about 16.3%). In order to find the points that fully satisfy these two objective functions, the objective function values corresponding to these optimization points are normalized, and then the norms of these function values are calculated. Among them, the optimum point is the point with the highest norm value. The point C is the optimal point obtained by this method.
6. Experiments
The above research methods are verified by engineering experiments, and 300 mm silicon single crystal is employed as an experimental object. The C point on the Pareto frontier is selected, and the corresponding process parameters are heating temperature , crystal rotation speed , melt rotation speed , and pulling speed . The experimental process of crystal growth is shown in Figure 8(a). From the start stage to the end of body stage of growth, the crystal growth is stable, which shows that the process parameters can meet the needs of large-size silicon single crystal growth. In order to measure the deformation of the solid-liquid interface, the end stage of crystal growth is abandoned. The crystal is pulled away from the silicon melt by rapid pulling, which ensures the visibility of the solid-liquid interface shape and facilitates the measurement of the deformation of the solid-liquid interface. The solid-liquid interface shape is shown in Figure 8(b). The deformation obtained by measuring the solid-liquid interface is about 12.9 mm. The relative error is only 3.1% compared to the optimized results. The actual engineering practice proves that the optimized result of the proposed hybrid strategy has higher precision and is more suitable for the optimization of high-quality silicon single crystal growth process parameters.
(a)
(b)
7. Conclusion
A hybrid strategy including CFD, GMDH, and improved NSGA-II for the process parameter optimization of Cz crystal growth is proposed in the paper. In order to obtain the optimal process parameters, the deformation of solid-liquid interface and the defect evaluation criteria are selected as the optimization functions. The polynomial model of objective function is identified by the CFD and GMDH neural network algorithm. An improved NSGA-II is proposed to obtain Pareto optimal solution of production process parameters. The hybrid strategy proposed in this paper combines the numerical simulation method with an advanced intelligent algorithm, which has the advantage of avoiding the limitation of numerical simulation results, transforming the complex crystal growth process model into a clear mathematical expression by the system identification strategy, and establishing the correspondence relationship between crystal quality and growth process parameters through the intelligent optimization. The experiments prove that the hybrid strategy is a new method to obtain accurate crystal growth process parameters, which can be applied to solve multiobjective optimization problems of other complex systems with uncertain models.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by the National Natural Science Foundation (NNSF) of China (62073258, 62003261, and 62127809), Natural Science Basic Research Program of Shaanxi Province of China (No. 2020JQ-650), and Scientific Research Program Funded of Shaanxi Education Department (No. 20JK0788).