Crystal bar manufacturing ending optimization control method

Abstract

The invention provides a crystal bar manufacturing ending optimization control method which comprises the following steps: calculating a temperature field distribution change rule according to a parameter coupling strength matrix, simulating spatial distribution characteristics of a melt temperature field under different parameter combinations by adopting a finite element analysis method, and obtaining growth interface temperature gradient and heat flux density distribution data; a crystal integrity evaluation index is calculated through a process parameter combination scheme, if the crystal integrity index is lower than a set standard, the parameter coupling relation is optimized again, an ending time control strategy is adjusted, and a final multi-parameter cooperative control sequence is determined; real-time parameter adjustment operation is executed according to the multi-parameter cooperative control sequence, the mechanical strength change trend and the process window stability are monitored, and the optimal process control scheme of high-quality cone ending is obtained.

Description

A method for optimizing the control of the finishing stage in crystal rod manufacturing Technical Field

[0001] This invention relates to the field of semiconductor technology, and in particular to a method for optimizing and controlling the final stages of crystal rod manufacturing. Background Technology

[0002] Crystal ingot growth, as a core process in the manufacturing of modern semiconductors and optoelectronic devices, directly determines the quality and industrialization level of key materials such as silicon wafers and gallium arsenide. In this manufacturing process, the crystal ingot is slowly pulled from the melt into shape, and the final stage of growth, namely the finishing operation, has a decisive impact on the material utilization rate of the entire crystal ingot and the efficiency of subsequent processing.

[0003] Current ingot finishing processes commonly suffer from isolated parameter adjustments. Technicians often focus only on optimizing a single parameter, neglecting the synergistic relationship between multiple process parameters such as pulling speed, heating power, and melt level. This fragmented parameter control approach results in a lack of organic coordination among various process stages during the finishing process, making it difficult to establish a stable and controllable process window.

[0004] The dynamic coupling characteristics of process parameters during the final stage pose a significant challenge to traditional control methods. The gradual reduction in pulling speed inevitably leads to a redistribution of the melt temperature field, affecting the stability of the liquid level. Fluctuations in the liquid level, in turn, influence the morphological control of the crystal growth interface. This interrelationship between parameters further exacerbates the non-uniformity of stress distribution. When stress concentration occurs in the final stage of the crystal rod, it not only generates dislocations and microcracks within the crystal but also causes the final taper angle to deviate from design requirements, ultimately affecting the mechanical strength and cutting efficiency of the crystal rod.

[0005] In actual production, operators often encounter the following dilemma: extending the finishing time to obtain the ideal tapered finishing angle results in a decrease in crystal quality; while attempting to shorten the finishing time makes it difficult to guarantee the integrity of the crystal in the finishing section. How to achieve rapid, high-quality tapered finishing through multi-parameter coordinated control while ensuring the integrity and mechanical strength of the crystal in the finishing section has become a key issue in optimizing the crystal rod growth process. Summary of the Invention

[0006] This invention provides a method for optimizing and controlling the finishing stage of crystal rod manufacturing, mainly including: acquiring real-time monitoring data of pulling speed, heating power, and liquid level height during the finishing stage of crystal rod manufacturing; analyzing the coupling relationship between various process parameters through a multi-parameter collaborative optimization algorithm, establishing a mathematical model of the mutual influence between parameters, and obtaining the parameter coupling strength matrix and dynamic correlation coefficient; calculating the temperature field distribution variation law based on the parameter coupling strength matrix, simulating the spatial distribution characteristics of the melt temperature field under different parameter combinations using the finite element analysis method, and obtaining the temperature gradient and heat flux density distribution data of the growth interface; establishing a stress prediction model based on the temperature gradient data of the growth interface, calculating the stress distribution state inside the finishing stage of the crystal rod, and if the stress concentration area exceeds a preset threshold range... The matching relationship between heating power and drawing speed is adjusted to obtain the optimized stress distribution prediction result. Based on the stress distribution prediction result, the control target value of the finishing cone angle is determined. An adaptive control algorithm is used to dynamically adjust the coordinated change trajectory of liquid level height and drawing speed to obtain a process parameter combination scheme that meets the cone angle requirements. The crystal integrity evaluation index is calculated through the process parameter combination scheme. If the crystal integrity index is lower than the set standard, the parameter coupling relationship is re-optimized and the finishing time control strategy is adjusted to determine the final multi-parameter collaborative control sequence. Real-time parameter adjustment operations are performed according to the multi-parameter collaborative control sequence to monitor the mechanical strength change trend and process window stability, thereby obtaining the optimal process control scheme for high-quality cone finishing.

[0007] The technical solution provided by this invention can include the following beneficial effects: This invention discloses an optimized control method for the final stage of crystal rod manufacturing, addressing key technical challenges such as complex coupling of process parameters, difficulty in accurately controlling stress distribution, and unstable taper angle during the final stage of crystal growth. This invention establishes a mathematical model and strength matrix of the coupling relationship between parameters by real-time monitoring of key parameters such as pulling speed, heating power, and liquid level. It uses finite element analysis to simulate the melt temperature field distribution characteristics under different parameter combinations, obtaining temperature gradient and heat flux density data at the growth interface. Based on the temperature gradient data, a stress prediction model is established. When stress concentration exceeds a threshold, the matching relationship of process parameters is automatically adjusted. An adaptive control algorithm dynamically optimizes the coordinated change trajectory of liquid level height and pulling speed, ensuring the accuracy of taper angle control. This invention achieves intelligent coordinated control of the crystal rod finalization process, effectively reducing internal stress concentration, improving crystal integrity and mechanical strength, and obtaining an optimal process control scheme for high-quality taper finalization. Attached Figure Description

[0008] Figure 1 is a flowchart of a crystal rod manufacturing end-of-life optimization control method according to the present invention.

[0009] Figure 2 is a schematic diagram of a crystal rod manufacturing end-of-life optimization control method according to the present invention. Detailed Implementation

[0010] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0011] As shown in Figures 1-2, the crystal rod manufacturing end-of-process optimization control method of this embodiment may specifically include: Step S101, acquiring real-time monitoring data of pulling speed, heating power, and liquid level height during the crystal rod end-of-process, analyzing the coupling relationship between various process parameters through a multi-parameter collaborative optimization algorithm, establishing a mathematical model of mutual influence between parameters, and obtaining the parameter coupling strength matrix and dynamic correlation coefficient; the specific calculation formula is as follows: ;

[0012] This formula represents the method for calculating the dynamic correlation coefficient between parameters, where This represents the dynamic correlation coefficient between parameters i and j at time t. This represents the change of parameter i at time tk. This indicates that parameter j takes time delay into account. The change after that, where n represents the length of the time window. This represents the coupling delay parameter between parameters i and j. Real-time monitoring data from the pulling speed sensor, heating power detector, and liquid level measuring instrument are acquired during the ingot finishing process. The raw monitoring data is filtered and denoised using a data preprocessing module to obtain a standardized monitoring dataset. A sliding time window is used to segment the standardized monitoring dataset, and the variation amplitude and fluctuation frequency of each process parameter are extracted based on time series characteristics to determine the pulling speed change rate, heating power adjustment, and liquid level offset. The temporal correlation between the pulling speed change rate and the heating power adjustment is analyzed using a parameter change rate calculation module. If the pulling speed change rate exceeds a preset threshold, the response delay time of the heating power adjustment at the corresponding moment is calculated to obtain the speed-power coupling delay parameter. A ternary relationship function between the liquid level offset and the first two parameters is established based on the speed-power coupling delay parameter. The least squares method is used to fit the nonlinear relationship curve in the three-dimensional parameter space to obtain a set of mathematical model equations for the interaction between parameters. The specific equations are as follows: This formula represents a ternary function relating the liquid level to the drawing speed and heating power, where h(t) represents the liquid level at time t. Indicates consideration of time delay The pulling speed, Indicates consideration of time delay Heating power, , , , The coefficients are obtained by fitting using the least squares method. The random error term is represented. Partial derivatives between parameter pairs are calculated using a system of mathematical model equations to construct a third-order parameter coupling strength matrix. Matrix elements represent the strength of interaction between parameters, determining dominant and secondary coupling relationships. Dynamic correlation coefficients are calculated using the numerical values ​​in the parameter coupling strength matrix. If the coupling strength of a parameter pair exceeds a critical value, the dynamic correlation coefficient of that parameter pair is set to a strongly correlated state, obtaining complete parameter coupling analysis results. Based on the parameter coupling analysis results, a collaborative control strategy for the process parameters of the crystal rod finishing process is generated, outputting the optimal adjustment sequence and control timing for pulling speed, heating power, and liquid level height.

[0013] Specifically, during the ingot finishing process, a high-precision sensor network was used to collect real-time monitoring data for the pulling speed (1.2 mm / min), heating power (45.8 kW), and liquid level (12.5 mm). The sampling frequency was set to 10 Hz to ensure data continuity. A multi-parameter collaborative optimization framework combining particle swarm optimization and genetic algorithms was used, with a particle swarm size of 50 individuals, 300 iterations, a crossover probability of 0.8, and a mutation probability of 0.1. The coupling relationship of the three key process parameters was analyzed. Pearson correlation coefficient calculations revealed a strong negative correlation between pulling speed and heating power (-0.76), while the correlation coefficient between liquid level and pulling speed was 0.63, showing a moderate positive correlation. A multivariate nonlinear regression mathematical model was established: V = α1P^(-0.45) + α2H^(0.32) + β, where V is the drawing speed, P is the heating power, and H is the liquid level. The least squares fitting yielded α1 = 2.34, α2 = 0.87, and β = 0.15. A 3×3 parameter coupling strength matrix was constructed, with diagonal elements equal to 1 and off-diagonal elements of 0.76, 0.63, and 0.52, reflecting the degree of mutual influence among the parameters. The sliding window method was used to calculate the dynamic correlation coefficients, with a window length of 100 data points. The window was slid every 10 points to obtain a time-varying correlation coefficient sequence. The dynamic correlation coefficient between drawing speed and heating power fluctuated between 0.72 and 0.81, reflecting the dynamic characteristics of parameter coupling relationships during the final stage.

[0014] Step S102: Calculate the temperature field distribution variation law based on the parameter coupling strength matrix, and use the finite element analysis method to simulate the spatial distribution characteristics of the melt temperature field under different parameter combinations to obtain the temperature gradient and heat flux density distribution data at the growth interface; the specific calculation formula is as follows: ;

[0015] G represents the magnitude of the temperature gradient at a certain location on the growth interface, ΔT represents the temperature difference between adjacent nodes, and Δn represents the node spacing. and This represents the temperature values ​​of two adjacent grid nodes. This represents the spatial distance between the two nodes; this formula is used to calculate the spatial distribution of the temperature gradient at the growth interface.

[0016] The radial heat flux density is represented by k, the thermal conductivity coefficient of the material is represented by T, and the radial coordinate is represented by r. and The temperature values ​​at different radial positions, r and This represents the corresponding radial position coordinates; the formula is based on Fourier's law to calculate the radial heat flux density components at each grid node.

[0017] The numerical weights of each element in the parameter coupling strength matrix are obtained, and the boundary condition parameters for temperature field calculation are determined based on the matrix element sizes. A geometric model and material property database for the melt region are then established. A triangular mesh generator is used to spatially discretize the melt region, and the mesh density distribution is adjusted according to the complexity of the geometric model. If the curvature near the growth interface changes drastically, the mesh node density in that region is increased to obtain a high-precision finite element mesh structure. Using the coordinate information of each node in the finite element mesh structure, combined with the boundary condition parameters, a system of partial differential equations for the temperature field is established. The Gauss-Seidel iterative method is used to solve for the temperature values ​​of each mesh node, obtaining the temperature distribution matrix of the entire melt region. The temperature difference between adjacent nodes on the growth interface is calculated based on the temperature distribution matrix. The ratio of the temperature difference to the node spacing is used to determine the magnitude of the temperature gradient at each location, obtaining the spatial distribution data of the temperature gradient at the growth interface. Fourier's law is used to calculate the heat flux density vector at each mesh node. The radial heat flux density value is extracted from the radial component of the heat flux density vector. If the radial heat flux density of a node exceeds the material critical value, the node is marked as a thermal stress concentration region, determining the distribution law of the radial heat flux density. By analyzing the radial heat flux density distribution pattern, the convective heat transfer mode inside the melt is analyzed, the rate of change of heat flux density at different radial positions is calculated, and the direction of heat transfer is determined based on the positive or negative sign of the rate of change of heat flux density, thus obtaining complete data on the temperature field distribution variation.

[0018] Specifically, based on the constructed parameter coupling strength matrix, a three-dimensional melt temperature field calculation model was established using ANSYS finite element analysis software. The mesh was composed of tetrahedral elements, with a total of 156,780 nodes and 892,340 elements. The values ​​0.76, 0.63, and 0.52 in the coupling strength matrix were used as boundary condition weighting coefficients. The crucible bottom temperature was set at 1687 K, the sidewall temperature at 1645 K, and the melt surface temperature at 1683 K. The nonlinear heat conduction equation was solved using the Newton-Raphson iterative algorithm, with a convergence accuracy of 1 × 10⁻⁶ and a maximum iteration step of 500. The calculation results show that the temperature gradient at the growth interface is non-uniformly distributed, with a radial temperature gradient of 23.5 K / cm, an axial temperature gradient of 18.2 K / cm, a central temperature of 1685.3 K, and an edge temperature of 1681.7 K. The heat flux density distribution was calculated using Fourier's law, yielding a heat flux density of 4.8 × 10⁵ W / m² at the interface center and 3.2 × 10⁵ W / m² at the edge. To verify the accuracy of the temperature field distribution, the calculated results were compared with infrared thermal imager measurements, with the temperature deviation controlled within ±2.1 K. Based on the temperature field distribution data, a temperature isosurface was constructed using the Kriging interpolation algorithm, achieving an interpolation accuracy of 99.2%. This generated a three-dimensional temperature field distribution map containing 127 isotherms, providing thermodynamic basis data for subsequent crystal growth quality prediction.

[0019] Step S103: Establish a stress prediction model using the temperature gradient data at the growth interface, calculate the stress distribution within the end section of the crystal rod, and if the stress concentration area exceeds a preset threshold range, adjust the matching relationship between the heating power and the pulling speed to obtain the optimized stress distribution prediction result; the specific calculation method is as follows: ;

[0020] This is the core formula of the stress prediction model based on the theory of elasticity. Let E represent the stress tensor components, E represent the elastic modulus, and ν represent Poisson's ratio. Represents the components of the strain tensor. Indicates volumetric strain. The symbol represents Kronecker, α represents the coefficient of thermal expansion, and ΔT represents the temperature change. This formula is used to calculate the principal stress components of each node inside the end section of the crystal rod.

[0021] The temperature change rate values ​​at each measurement point in the growth interface temperature gradient data are obtained. Based on the spatial distribution characteristics of the temperature change rate, a basic parameter matrix for stress field calculation is established. A stress prediction model framework for the crystal rod's terminal section is constructed using elasticity theory. The principal stress components at each node within the crystal rod's terminal section are calculated using the stress prediction model. The spatial mapping relationship of the stress distribution state is determined based on the magnitude and direction of the principal stress components. If the stress component at a node exceeds the material's yield strength, that node is marked as a stress concentration region. The distribution pattern of the stress concentration region is analyzed using a support vector machine algorithm. A regional risk assessment index is established based on the location coordinates and stress values ​​of the stress concentration region. The safety level of the current stress state is determined by comparing the risk assessment index with a preset threshold. The specific formula is as follows: ;

[0022] This is the formula for calculating the risk assessment index for stress concentration areas. This represents the risk assessment index for the k-th stress concentration region. This indicates the area of ​​the region. This indicates the number of nodes in the region. This represents the equivalent stress at the i-th node. Indicates the yield strength of the material. This represents the node weight coefficient; this formula is used to quantitatively assess the safety level of stress concentration areas.

[0023] Based on the safety level assessment results, the heating power output parameters are adjusted. If the safety level is lower than the preset standard, the heating power value is reduced, and the drawing speed control parameters are adjusted simultaneously to obtain a new matching relationship between heating power and drawing speed. The temperature field distribution data of the ingot's terminal section is recalculated using this new matching relationship. The boundary condition parameters in the stress prediction model are corrected based on the updated temperature field data to obtain the corrected stress distribution calculation results. The corrected stress distribution data is used to verify the changes in the stress concentration region. The quantitative indicators of the optimization effect are evaluated based on the area of ​​the stress concentration region and the reduction in peak stress, determining the final optimized stress distribution prediction results.

[0024] Specifically, based on the distribution data of the radial temperature gradient (23.5 K / cm) and axial temperature gradient (18.2 K / cm) at the growth interface, a stress prediction model for the terminal section of the crystal rod is established using thermoelastic coupling theory. The temperature gradient is input into the finite element solver as a thermal stress source term. The model sets the elastic modulus of silicon crystal to 130 GPa, Poisson's ratio to 0.28, and the coefficient of thermal expansion to 2.6 × 10⁻⁶ / K. The terminal section of the crystal rod is 85 mm long, with its diameter gradually decreasing from 200 mm to zero. Solving the thermoelastic equations using the Lagrange multiplier method, the maximum principal stress inside the terminal section is calculated to be 156 MPa, located radially 12 mm from the growth interface. The peak shear stress is 89 MPa, occurring in the middle region of the terminal section. Since the stress concentration area exceeds the preset threshold of 150 MPa, the system automatically initiates a parameter optimization algorithm, employing a genetic algorithm to jointly optimize the heating power and drawing speed. With an initial heating power of 45kW and a drawing speed of 1.2mm / h, after 200 generations of iteration using a genetic algorithm with a population size of 50, the optimized parameter combination was obtained: heating power of 42.3kW and drawing speed of 0.95mm / h. Substituting these optimized parameters back into the stress prediction model, the maximum principal stress in the final stage was reduced to 138MPa, the peak shear stress was reduced to 76MPa, and the stress distribution uniformity was improved by 23.4%, meeting the requirements of the crystal growth process.

[0025] Step S104: Based on the stress distribution prediction results, determine the target value for the final taper angle. Use an adaptive control algorithm to dynamically adjust the coordinated change trajectory of the liquid level height and drawing speed to obtain a combination of process parameters that meets the taper angle requirements. The specific formula is as follows:

[0026] in, This indicates the target value for cone angle control. This represents the theoretical value of the cone angle calculated based on the stress distribution. From the perspective of process standards, This represents the angle deviation correction factor. This represents the weight coefficient of the i-th node. This represents the stress value at the i-th node. Indicates the average stress value. Let represent the radial position of the i-th node, and n represent the total number of nodes. Obtain the stress values ​​of each node from the stress distribution prediction results. Calculate the theoretical value of the conical angle of the final segment based on the radial distribution pattern of the stress values. Determine the angle control target value by analyzing the deviation between the theoretical conical angle value and the process standard angle. Establish the liquid level adjustment range based on the angle control target value. Use an adaptive control algorithm to analyze the correlation between the current liquid level and the target angle. If the liquid level deviates from the preset range, calculate the correction value for the liquid level. The specific calculation is as follows: ;

[0027] ΔH represents the liquid level correction value. Indicates the liquid level adjustment coefficient. Indicates the current liquid level. This represents the center value of the preset liquid level height range, and α represents the dynamic adjustment parameter. t represents the angular deviation, and t represents time.

[0028] The direction of drawing speed adjustment is determined by the liquid level correction value. Constraints on speed change are established based on the coupling relationship between liquid level and drawing speed, thus obtaining the range of drawing speed variation. A mathematical model of the coordinated change trajectory is constructed using the range of drawing speed variation. Based on the trajectory model, the synchronous change curves of liquid level and drawing speed over time are calculated, resulting in discretized data points for the coordinated change trajectory. The specific model is as follows: ;

[0029] This represents the pulling speed at time t. β represents the base drawing speed, β represents the speed adjustment coupling coefficient, ΔH represents the liquid level height correction amount, and ω represents the angular frequency of the coordinated change. Let λ represent the phase angle, λ represent the attenuation coefficient, and t represent the time variable.

[0030] Candidate combinations of process parameters are extracted from the discretized data points of the coordinated change trajectory. The feasibility of each combination is analyzed by examining the numerical ratios of liquid level height and drawing speed within these combinations, and the process parameter combination that meets the cone angle control requirements is selected. The actual formed value of the cone angle in the final section is calculated using the selected process parameter combinations. The accuracy level of each combination is evaluated based on the error between the actual formed value and the control target value, and the optimal process parameter combination is determined. The liquid level height and drawing speed values ​​from the optimal process parameter combination are used to update the execution instructions of the control system, and a complete control sequence for the forming of the final section is generated based on the parameter settings of the execution instructions.

[0031] Specifically, based on the stress distribution calculation results of the final stage, the system uses a geometric optimization model to determine the optimal cone angle control target. By analyzing the correlation between stress concentration location and cone angle, a cone angle-stress response function relationship is established. The model sets the ideal cone angle range to 15° to 25°. The currently measured cone angle is 28.3°, exceeding the optimal range and causing stress concentration. The system employs a fuzzy PID adaptive control algorithm, using the cone angle deviation as the control input and the liquid level height and drawing speed as coordinated output variables. The controller is set with a proportional coefficient Kp of 0.85, an integral coefficient Ki of 0.12, a derivative coefficient Kd of 0.03, and a sampling period of 30 seconds. The algorithm monitors the change in liquid level height in the crucible in real time. The current liquid level height is 156.8 mm, and the target is adjusted to 162.5 mm. Simultaneously, the drawing speed is adjusted from the current 1.15 mm / h to 0.88 mm / h. The control system employs a Kalman filter to reduce noise in the liquid level sensor signal, with the filter gain matrix K set to 0.75 and the state prediction error covariance set to 0.02. After a 45-minute adaptive adjustment process, the liquid level and pulling speed reach a coordinated equilibrium, the cone angle converges to 22.1°, within the target control range, and the geometric uniformity index of the finishing section improves to 0.92.

[0032] Step S105: Calculate the crystal integrity evaluation index using the process parameter combination scheme. If the crystal integrity index is lower than the set standard, re-optimize the parameter coupling relationship and adjust the end-of-life time control strategy to determine the final multi-parameter collaborative control sequence. The specific calculation method is as follows: ;

[0033] V represents the crystal integrity evaluation index, where V represents the crystal volume. This represents the defect density distribution function within the crystal. Represents the temperature gradient distribution. Indicates reference temperature. The crystallization rate is represented by α, and the rate influence coefficient is represented by α. This formula calculates the crystal integrity evaluation index through the correlation analysis of defect density distribution with temperature gradient and crystallization rate.

[0034] The process involves acquiring temperature gradient and crystallization rate data from the combined process parameters. Based on the temperature gradient data, the internal defect density distribution of the crystal is calculated. A crystal integrity evaluation index is obtained through correlation analysis between the defect density distribution and the crystallization rate data. The current process status is determined based on the difference between the crystal integrity evaluation index and the preset standard. If the evaluation index is lower than the preset standard, the temperature control parameters and time control parameters in the parameter coupling relationship are extracted to obtain a correction benchmark for the parameter coupling relationship. The matching degree between the temperature control parameters and the drawing time is recalculated using the correction benchmark. Matching degree analysis is used to determine the adjustment direction of the finishing time control strategy, obtaining the optimized range of the finishing time. A time series control model is established based on the optimized range of the finishing time. The parameter change requirements at each time node are analyzed using the time series control model to obtain the collaborative change law of multiple parameters in the time dimension. A control sequence generation framework is constructed using the collaborative change law of multiple parameters. The synchronous adjustment amplitude between parameters is calculated based on the temporal relationship of each parameter in the generation framework, determining the discretization nodes of the multi-parameter collaborative control sequence. The specific calculations are as follows:

[0035] This indicates the synchronization adjustment range between parameters, where m represents the total number of parameters. This represents the value of the j-th parameter at time t. This represents the rate of change of the j-th parameter. This represents the coupling coefficient between parameter j and parameter k; this formula is used to calculate the degree of synchronization of changes in each parameter in multi-parameter coordinated control, ensuring coordinated changes between parameters.

[0036] The parameter setpoints for each time period are extracted by discretizing the nodes of the multi-parameter collaborative control sequence. The executability of the control sequence is verified by the matching of the parameter setpoints with the process constraints, resulting in the final multi-parameter collaborative control sequence. A complete control instruction set for the crystal growth process is generated based on the final multi-parameter collaborative control sequence. The automated control process for the final crystal forming stage is completed by setting the execution order and values ​​of each parameter in the control instruction set.

[0037] Specifically, the system employs a multi-dimensional crystal quality assessment model to perform integrity analysis on the combination of process parameters. A comprehensive evaluation system is constructed using three dimensions: X-ray diffraction peak intensity measurement, dislocation density detection, and lattice constant deviation calculation. Currently, the measured full width at half maximum (FWHM) of the crystal diffraction peak is 0.087°, the dislocation density is 2.3 × 10⁴ cm⁻², and the lattice constant deviation rate is 0.15%. The calculated comprehensive integrity index is 0.73, lower than the system's set standard threshold of 0.85. Based on this evaluation result, the system initiates a parameter coupling relationship reconstruction algorithm. A genetic algorithm is used to optimize the nonlinear coupling coefficient between the temperature gradient and the rotation speed, adjusting the temperature gradient from the current 12.5 K / cm to 9.8 K / cm and reducing the crystal rotation speed from 8.5 rpm to 6.2 rpm. Simultaneously, an oscillation frequency compensation mechanism is introduced, setting the oscillation amplitude to ±0.3 rpm and the period to 180 seconds. The finishing time control strategy adopts a segmented decreasing pattern, re-dividing the original 90-minute finishing process into three stages: maintaining a pulling speed of 0.65 mm / h for the initial 30 minutes, linearly decreasing to 0.42 mm / h for the middle 35 minutes, and using an exponential decay function to control the speed to 0.18 mm / h for the final 25 minutes. The system calculates the multi-parameter synergistic effect in real time through a neural network prediction model. The hidden layer of the network has 16 neurons, and the learning rate is set to 0.008. After 1200 iterations of training, the prediction accuracy reaches 96.2%.

[0038] Step S106: Perform real-time parameter adjustment operations according to the multi-parameter collaborative control sequence, monitor the trend of mechanical strength change and process window stability, and obtain the optimal process control scheme for high-quality conical tailing.

[0039] Temperature and pulling speed parameters from a multi-parameter collaborative control sequence are acquired. A real-time monitoring framework is established based on the parameter adjustment time intervals. Mechanical strength data during the crystal pulling process is collected through this framework. The gradient value of strength change is calculated based on the mechanical strength data. The gradient value is used to analyze the fluctuation amplitude of the strength evolution trend. If the fluctuation amplitude exceeds a preset threshold, the current process parameter status is extracted. The specific calculation is as follows:

[0040] The gradient value represents the change in mechanical strength, S(t) represents the mechanical strength data at time t, and Δt represents the time interval for parameter adjustment. This formula is used to calculate the fluctuation range of the strength evolution trend. When the gradient value exceeds the preset threshold, the process parameter status extraction is triggered.

[0041] By analyzing the matching relationship between temperature field distribution and drawing speed through process parameter status analysis, the boundary range data of the process window is obtained, and the stable state of the process window is determined based on the boundary range data. If the stable state of the process window deviates from the standard range, a Kalman filter algorithm is used to correct the parameter adjustment sequence; the specific algorithm is as follows:

[0042] P represents the error covariance matrix. Represents the state transition matrix. Let k represent the process noise covariance matrix and k represent the time step. This Kalman filter algorithm formula is used to correct the parameter adjustment sequence when the steady state of the process window deviates from the standard range.

[0043] The temperature and speed control curves are updated based on the corrected calculation results. The geometric parameters of the conical tail section are determined through intersection analysis of the control curves, obtaining the correlation data between the tail shape and the mechanical strength distribution. A strength prediction model is constructed using this correlation data, and the strength distribution under different tail angles is calculated using the prediction model. Based on the strength distribution, a tail angle range that meets the mechanical performance requirements is selected. The execution sequence of parameter adjustments is recalculated using the angle range constraints, resulting in an optimized control sequence execution scheme. The specific calculations are as follows: ;

[0044] This represents the optimal termination angle, where M represents the number of angle sampling points. Indicates the weighting coefficient. This represents the intensity distribution calculated by the prediction model. This represents the strength value required for mechanical performance; the formula is used to select the range of finishing angles that meet the mechanical performance requirements based on the strength distribution.

[0045] The operating status of the drawing equipment is adjusted using the parameter settings in the execution plan, and the effect of parameter adjustment is verified through real-time feedback of the equipment status. The degree of improvement in process window stability is evaluated based on the parameter adjustment effect. If the degree of stability improvement reaches the expected standard, the current parameter adjustment strategy is solidified to obtain the optimal process control scheme for the conical tail section.

[0046] Specifically, the system initiates a real-time adjustment execution module based on a defined multi-parameter collaborative control sequence. It collects stress distribution data at the bottom of the crucible through a distributed sensor network, measuring a radial stress peak of 145.2 MPa, axial stress of 89.7 MPa, and tangential stress of 112.4 MPa. The comprehensive mechanical strength coefficient is calculated to be 0.847 using a finite element analysis algorithm. The system employs a Kalman filter algorithm to predict and analyze stress change trends, setting the state transition matrix to 3×3 dimensions, the observation noise variance to 0.025, and the diagonal elements of the process noise covariance matrix to 0.012. After 240 seconds of continuous data fusion processing, the predicted stress change slope for the next 600 seconds is -0.032 MPa / min. For process window stability monitoring, a multivariate statistical process control method is used, establishing a T² statistic monitoring model with a control limit of 15.68. The currently calculated T² value is 12.34, which is within the stable control range. The system simultaneously monitors the melt interface shape parameters. The interface curvature radius is measured to be 8.7 mm and the cone angle to be 28.5° using an optical imaging system. Combined with a fluid dynamics simulation model, the interface stability index is calculated to be 0.923. Based on real-time monitoring data, the system employs a fuzzy logic controller to dynamically adjust the heating power distribution ratio, increasing the upper heater power from 2.8 kW to 2.5 kW, the lower heater power from 3.2 kW to 3.6 kW, and keeping the side heater power constant at 1.9 kW. Precise temperature field control is achieved through a PID control algorithm, with a proportional gain of 0.45, an integral time constant of 85 seconds, and a derivative time constant of 12 seconds.

[0047] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application.

Claims

1. A method for optimizing and controlling the final stages of crystal rod manufacturing, characterized in that, The method includes: acquiring real-time monitoring data of pulling speed, heating power, and liquid level during the ingot finishing process; analyzing the coupling relationship between various process parameters through a multi-parameter collaborative optimization algorithm; establishing a mathematical model of the mutual influence between parameters; obtaining the parameter coupling strength matrix and dynamic correlation coefficient; calculating the temperature field distribution variation law based on the parameter coupling strength matrix; simulating the spatial distribution characteristics of the melt temperature field under different parameter combinations using the finite element analysis method; acquiring the temperature gradient and heat flux density distribution data of the growth interface; establishing a stress prediction model based on the temperature gradient data of the growth interface; calculating the stress distribution state inside the ingot finishing section; and adjusting the heating power and pulling speed if the stress concentration area exceeds a preset threshold range. The speed matching relationship is used to obtain the optimized stress distribution prediction result. Based on the stress distribution prediction result, the control target value of the finishing cone angle is determined. An adaptive control algorithm is used to dynamically adjust the coordinated change trajectory of the liquid level height and drawing speed to obtain a process parameter combination scheme that meets the cone angle requirements. The crystal integrity evaluation index is calculated through the process parameter combination scheme. If the crystal integrity index is lower than the set standard, the parameter coupling relationship is re-optimized and the finishing time control strategy is adjusted to determine the final multi-parameter coordinated control sequence. Real-time parameter adjustment operations are performed according to the multi-parameter coordinated control sequence to monitor the mechanical strength change trend and process window stability, thereby obtaining the optimal process control scheme for high-quality cone finishing.

2. The method for optimizing and controlling the finishing stage of crystal rod manufacturing according to claim 1, characterized in that, The process involves acquiring real-time monitoring data on pulling speed, heating power, and liquid level during the crystal rod finishing process. A multi-parameter collaborative optimization algorithm is used to analyze the coupling relationship between various process parameters, establishing a mathematical model of the mutual influence between parameters to obtain the parameter coupling strength matrix and dynamic correlation coefficients. The specific calculation formulas are as follows: This formula represents the method for calculating the dynamic correlation coefficient between parameters, where This represents the dynamic correlation coefficient between parameters i and j at time t. This represents the change of parameter i at time tk. This indicates that parameter j takes time delay into account. The change after that, where n represents the length of the time window. This describes the coupling delay parameter between parameters i and j. It includes acquiring real-time monitoring data from the pulling speed sensor, heating power detector, and liquid level measuring instrument during the ingot finishing process. The raw monitoring data is filtered and denoised using a data preprocessing module to obtain a standardized monitoring dataset. A sliding time window is used to segment the standardized monitoring dataset, and the variation amplitude and fluctuation frequency of each process parameter are extracted based on time series characteristics to determine the pulling speed change rate, heating power adjustment, and liquid level offset. The parameter change rate calculation module analyzes the temporal correlation between the pulling speed change rate and the heating power adjustment. If the pulling speed change rate exceeds a preset threshold, the response delay time of the heating power adjustment at the corresponding moment is calculated to obtain the speed-power coupling delay parameter. Based on the speed-power coupling delay parameter, a ternary relationship function between the liquid level offset and the first two parameters is established. The least squares method is used to fit the nonlinear relationship curve in the three-dimensional parameter space to obtain a set of mathematical model equations for the interaction between parameters. The specific equations are as follows: This formula represents a ternary function relating the liquid level to the drawing speed and heating power, where h(t) represents the liquid level at time t. Indicates consideration of time delay The pulling speed, Indicates consideration of time delay Heating power, 、 、 、 The coefficients are obtained by fitting using the least squares method. The random error term is represented; the partial derivatives between each parameter pair are calculated through the mathematical model equations to construct a third-order parameter coupling strength matrix, where the matrix elements represent the strength of the interaction between parameters, thus determining the dominant and secondary coupling relationships; The dynamic correlation coefficient is calculated numerically using the parameter coupling strength matrix. If the coupling strength of a parameter pair exceeds the critical value, the dynamic correlation coefficient of that parameter pair is set to a strong correlation state, thus obtaining complete parameter coupling analysis results. Based on the parameter coupling analysis results, a collaborative control strategy for the process parameters of the crystal rod finishing process is generated, and the optimal adjustment sequence and control timing for pulling speed, heating power, and liquid level are output.

3. The method for optimizing and controlling the finishing stage of crystal rod manufacturing according to claim 1, characterized in that, The temperature field distribution variation law is calculated based on the parameter coupling strength matrix, and the spatial distribution characteristics of the melt temperature field under different parameter combinations are simulated using the finite element analysis method to obtain the temperature gradient and heat flux density distribution data of the growth interface; the specific calculation formula is as follows: G represents the magnitude of the temperature gradient at a certain location on the growth interface, ΔT represents the temperature difference between adjacent nodes, and Δn represents the node spacing. and This represents the temperature values ​​of two adjacent grid nodes. This represents the spatial distance between the two nodes; the formula is used to calculate the spatial distribution of the temperature gradient at the growth interface. ； The radial heat flux density is represented by k, the thermal conductivity coefficient of the material is represented by T, and the radial coordinate is represented by r. and The temperature values ​​at different radial positions, r and Indicates the corresponding radial position coordinates; This formula calculates the radial heat flux density components at each grid node based on Fourier's law. It includes: obtaining the numerical weights of each element in the parameter coupling strength matrix; determining the boundary condition parameters for temperature field calculation based on the matrix element sizes; establishing a geometric model and material property database for the melt region; using a triangular mesh generator to spatially discretize the melt region; adjusting the mesh density distribution according to the complexity of the geometric model; increasing the mesh node density in the region near the growth interface if the curvature changes drastically to obtain a high-precision finite element mesh structure; establishing a system of partial differential equations for the temperature field using the coordinate information of each node in the finite element mesh structure and the boundary condition parameters; and solving for the temperature values ​​at each grid node using the Gauss-Seidel iteration method to obtain the temperature values ​​for the entire melt region. The temperature distribution matrix is ​​used to calculate the temperature difference between adjacent nodes on the growth interface. The ratio of the temperature difference to the node spacing is used to determine the magnitude of the temperature gradient at each location, thus obtaining the spatial distribution data of the temperature gradient at the growth interface. Fourier's law is used to calculate the heat flux density vector at each grid node. The radial heat flux density value is extracted from the radial component of the heat flux density vector. If the radial heat flux density of a node exceeds the material's critical value, the node is marked as a thermal stress concentration area, thus determining the distribution law of the radial heat flux density. The convective heat transfer mode inside the melt is analyzed through the radial heat flux density distribution law. The rate of change of heat flux density at different radial locations is calculated. The direction of heat transfer is determined based on the sign of the rate of change of heat flux density, thus obtaining complete data on the temperature field distribution variation law.

4. The method for optimizing and controlling the finishing stage of crystal rod manufacturing according to claim 1, characterized in that, The stress prediction model is established using temperature gradient data at the growth interface to calculate the stress distribution within the end section of the crystal rod. If the stress concentration area exceeds a preset threshold range, the matching relationship between heating power and pulling speed is adjusted to obtain an optimized stress distribution prediction result. The specific calculation method is as follows: This is the core formula of the stress prediction model based on the theory of elasticity. Let E represent the stress tensor components, E represent the elastic modulus, and ν represent Poisson's ratio. Represents the components of the strain tensor. Indicates volumetric strain. The symbol for Kronecker is given, α represents the coefficient of thermal expansion, and ΔT represents the change in temperature. This formula is used to calculate the principal stress components at each node within the end section of a crystal ingot. Specifically, it includes: obtaining the temperature change rate at each measurement point in the temperature gradient data of the growth interface; establishing a basic parameter matrix for stress field calculation based on the spatial distribution characteristics of the temperature change rate; constructing a stress prediction model framework for the end section of the crystal ingot using elasticity theory; calculating the principal stress components at each node within the end section of the crystal ingot using the stress prediction model; determining the spatial mapping relationship of the stress distribution state based on the magnitude and direction of the principal stress components; marking the node as a stress concentration region if the stress component at a certain node exceeds the material's yield strength; analyzing the distribution pattern of the stress concentration region using a support vector machine algorithm; establishing a regional risk assessment index based on the location coordinates and stress values ​​of the stress concentration region; and determining the safety level of the current stress state by comparing the risk assessment index with a preset threshold. The specific formula is as follows: This is the formula for calculating the risk assessment index for stress concentration areas. This represents the risk assessment index for the k-th stress concentration region. This indicates the area of ​​the region. This indicates the number of nodes in the region. This represents the equivalent stress at the i-th node. Indicates the yield strength of the material. This represents the node weight coefficient; this formula is used to quantitatively assess the safety level of the stress concentration region; the output parameters of the heating power are adjusted according to the assessment results of the safety level. If the safety level is lower than the preset standard, the heating power value is reduced and the control parameters of the drawing speed are adjusted simultaneously to obtain a new matching relationship combination of heating power and drawing speed; the temperature field distribution data of the crystal rod tail section is recalculated through the new matching relationship combination, and the boundary condition parameters in the stress prediction model are corrected according to the updated temperature field data to obtain the corrected stress distribution calculation results; the changes in the stress concentration region are verified using the corrected stress distribution data, and the quantitative indicators of the optimization effect are evaluated based on the area of ​​the stress concentration region and the degree of reduction of peak stress to determine the final optimized stress distribution prediction results.

5. The method for optimizing and controlling the finishing stage of crystal rod manufacturing according to claim 1, characterized in that, The control target value of the finishing cone angle is determined based on the stress distribution prediction results. An adaptive control algorithm is used to dynamically adjust the coordinated change trajectory of the liquid level height and drawing speed to obtain a combination of process parameters that meets the cone angle requirements. The specific formula is as follows: ;in, This indicates the target value for cone angle control. This represents the theoretical value of the cone angle calculated based on the stress distribution. From the perspective of process standards, This represents the angle deviation correction factor. This represents the weight coefficient of the i-th node. This represents the stress value at the i-th node. Indicates the average stress value. This represents the radial position of the i-th node, and n represents the total number of nodes. The process includes obtaining the stress values ​​of each node from the stress distribution prediction results, calculating the theoretical value of the conical angle of the final segment based on the radial distribution pattern of the stress values, determining the angle control target value through deviation analysis between the theoretical conical angle value and the process standard angle, establishing the liquid level adjustment range based on the angle control target value, using an adaptive control algorithm to analyze the correlation between the current liquid level height and the target angle, and calculating the correction value for the liquid level height if the liquid level height deviates from the preset range. The specific calculation is as follows: ΔH represents the numerical value of the liquid level correction. Indicates the liquid level adjustment coefficient. Indicates the current liquid level. This represents the center value of the preset liquid level height range, and α represents the dynamic adjustment parameter. Let t represent the angle deviation and t represent time. The direction of the drawing speed adjustment is determined by the liquid level correction value. Constraints on speed change are established based on the coupling relationship between liquid level and drawing speed, obtaining the range of drawing speed variation. A mathematical model of the coordinated change trajectory is constructed using the range of drawing speed variation. The synchronous change curves of liquid level and drawing speed over time are calculated based on the trajectory model, yielding discretized data points for the coordinated change trajectory. The specific model is as follows: ； This represents the pulling speed at time t. β represents the base drawing speed, β represents the speed adjustment coupling coefficient, ΔH represents the liquid level height correction amount, and ω represents the angular frequency of the coordinated change. Let λ represent the phase angle, λ represent the attenuation coefficient, and t represent the time variable. Candidate combinations of process parameters are extracted from the discretized data points of the coordinated change trajectory. The feasibility of each combination is analyzed by comparing the numerical ratios of liquid level height and drawing speed within the candidate combinations, and process parameter combinations that meet the cone angle control requirements are selected. The actual formed value of the cone angle of the final section is calculated using the selected process parameter combinations. The accuracy level of each combination scheme is evaluated based on the error between the actual formed value and the control target value, and the optimal process parameter combination scheme is determined. The execution instructions of the control system are updated using the liquid level height and drawing speed values ​​from the optimal process parameter combination scheme, and a complete control sequence for the forming of the final section is generated based on the parameter settings of the execution instructions.

6. The method for optimizing and controlling the finishing stage of crystal rod manufacturing according to claim 1, characterized in that, The crystal integrity evaluation index is calculated through a combination of process parameters. If the crystal integrity index is lower than the set standard, the parameter coupling relationship is re-optimized and the end-of-life control strategy is adjusted to determine the final multi-parameter collaborative control sequence. The specific calculation method is as follows: ； V represents the crystal integrity evaluation index, where V represents the crystal volume. This represents the defect density distribution function within the crystal. Represents the temperature gradient distribution. Indicates reference temperature. Let α represent the crystallization rate and α represent the rate influence coefficient. This formula calculates the crystal integrity evaluation index through correlation analysis between defect density distribution, temperature gradient, and crystallization rate. This includes acquiring temperature gradient and crystallization rate data from the process parameter combination scheme, calculating the internal defect density distribution of the crystal based on the temperature gradient data, obtaining the crystal integrity evaluation index through correlation analysis between the defect density distribution and crystallization rate data, judging the current process status based on the numerical difference between the crystal integrity evaluation index and the preset standard, and extracting the temperature control parameters and time control parameters from the parameter coupling relationship to obtain the correction benchmark for the parameter coupling relationship. The matching degree between temperature control parameters and drawing time is recalculated using a corrected benchmark based on parameter coupling. Matching degree analysis determines the adjustment direction of the finishing time control strategy, yielding the optimized range of the finishing time. A time-series control model is established based on this optimized range, and the parameter change requirements at each time node are analyzed to obtain the coordinated change pattern of multiple parameters over time. A control sequence generation framework is constructed using this coordinated change pattern. The synchronous adjustment amplitude between parameters is calculated based on the temporal relationship of each parameter in the generation framework, determining the discretized nodes of the multi-parameter coordinated control sequence. The specific calculations are as follows: ； This indicates the synchronization adjustment range between parameters, where m represents the total number of parameters. This represents the value of the j-th parameter at time t. This represents the rate of change of the j-th parameter. This represents the coupling coefficient between parameter j and parameter k. This formula is used to calculate the degree of synchronicity of the changes of each parameter in multi-parameter collaborative control, ensuring coordinated changes among parameters. The parameter setpoints for each time period are extracted through the discretized nodes of the multi-parameter collaborative control sequence. The executability of the control sequence is verified based on the matching of the parameter setpoints with the process constraints, resulting in the final multi-parameter collaborative control sequence. A complete control instruction set for the crystal growth process is generated based on the final multi-parameter collaborative control sequence. The automated control process for the final crystal forming stage is completed by setting the execution order and values ​​of each parameter in the control instruction set.

7. The method for optimizing and controlling the finishing stage of crystal rod manufacturing according to claim 1, characterized in that, The process of performing real-time parameter adjustment based on a multi-parameter collaborative control sequence, monitoring the trend of mechanical strength changes and the stability of the process window, and obtaining the optimal process control scheme for high-quality tapered crystal finishing includes: acquiring temperature parameters and pulling speed parameters from the multi-parameter collaborative control sequence; establishing a real-time monitoring framework based on the parameter adjustment time interval; collecting mechanical strength data during the crystal pulling process through the monitoring framework; calculating the gradient value of strength change based on the mechanical strength data; analyzing the fluctuation amplitude of the strength evolution trend using the gradient value; and extracting the process parameter status at the current moment if the fluctuation amplitude exceeds a preset threshold. The specific calculation is as follows: ； The gradient value represents the change in mechanical strength, S(t) represents the mechanical strength data at time t, and Δt represents the time interval for parameter adjustment. This formula is used to calculate the fluctuation amplitude of the strength evolution trend. When the gradient value exceeds a preset threshold, the process parameter state extraction is triggered. By analyzing the matching relationship between the temperature field distribution and the drawing speed through the process parameter state analysis, the boundary range data of the process window is obtained, and the stable state of the process window is determined based on the boundary range data. If the stable state of the process window deviates from the standard range, the Kalman filter algorithm is used to correct the parameter adjustment sequence. The specific algorithm is as follows: P represents the error covariance matrix. Represents the state transition matrix. Let represent the process noise covariance matrix, and k represent the time step. This Kalman filter algorithm formula is used to correct the parameter adjustment sequence when the process window's steady-state deviates from the standard range. The temperature and speed control curves are updated based on the corrected calculation results. The geometric parameters of the conical tail section are determined through intersection analysis of the control curves, obtaining the correlation data between the tail section shape and the mechanical strength distribution. A strength prediction model is constructed using this correlation data, and the strength distribution under different tail angles is calculated using the prediction model. Based on the strength distribution, a tail angle range that meets the mechanical performance requirements is selected. The execution sequence of parameter adjustments is recalculated using the angle range constraint to obtain the optimized control sequence execution scheme. The specific calculations are as follows: ； This represents the optimal termination angle, where M represents the number of angle sampling points. Indicates the weighting coefficient. This represents the intensity distribution calculated by the prediction model. This represents the strength value required for mechanical performance; the formula is used to screen the range of finishing angles that meet the mechanical performance requirements based on the strength distribution; the operating status of the drawing equipment is adjusted using the parameter settings in the execution plan, and the effect of parameter adjustment is verified through real-time feedback of the equipment status; the degree of improvement in the stability of the process window is evaluated based on the effect of parameter adjustment, and if the degree of stability improvement reaches the expected standard, the current parameter adjustment strategy is solidified to obtain the optimal process control plan for the conical finishing section.

Images