Method for solving heat flux of phase transition surface of czochralski silicon single crystal

Abstract

The application discloses a kind of direct pulling silicon single crystal phase transition surface heat flux solving method, first in the position of the observation window of silicon single crystal growth equipment installs temperature sensor, to measure and obtain the temperature of crystal surface measuring point, and establish silicon crystal two-dimensional axis symmetry heat conduction control equation and each boundary condition;Then the silicon crystal two-dimensional axis symmetry heat conduction control equation established is discretized, and iterative solution is solved using finite difference method, and the silicon crystal two-dimensional axis symmetry direct heat transfer model with boundary heat flux as input and measuring point temperature as output is obtained;Finally, an inverse heat transfer model capable of solving the phase transition surface heat flux is established, with the actual measured crystal surface temperature as input, and by continuously optimizing the heat flux information at the current time, the optimal estimation value of the heat flux at the current time is obtained as output. The application solves the problem that the phase transition surface heat flux cannot be directly measured and cannot be obtained in real time in the prior art.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of solving phase transition surface heat flux in the process of semiconductor silicon single crystal growth, and particularly relates to a Czochralski silicon single crystal phase transition surface heat flux solving method. BACKGROUND

[0002] Silicon single crystal is the main semiconductor material for preparing a series of electronic components and integrated circuits. With the development of high-tech fields, high-performance electronic components have higher requirements for the quality of single crystal silicon.

[0003] During the operation of the silicon single crystal production equipment, the internal environment is complex and changeable. Temperature is one of the important factors affecting the quality of the crystal, and temperature plays an important role in the crystal growth process. Due to the constraints of the production environment, a large number of sensors cannot be installed inside the furnace body, and only individual non-contact sensors are used to measure the temperature inside the furnace body. The temperature information obtained is very limited. Because the phase transition surface heat flux is the direct cause of the crystal temperature, the phase transition surface heat flux greatly affects the quality of the crystal, and the change of the heat flux will cause a large number of structural defects in the crystal. At present, there is no direct and effective method to obtain the change trend of the phase transition surface heat flux, and the method for solving the phase transition surface heat flux has the characteristics of complex structure, large model and time-lagging calculation, which leads to great limitations in engineering application. Therefore, it is very important to take an accurate, fast and effective method to solve the change of the phase transition surface heat flux for crystal growth. SUMMARY

[0004] The purpose of the present application is to provide a Czochralski silicon single crystal phase transition surface heat flux solving method, which solves the problem that the phase transition surface heat flux cannot be directly measured and cannot be obtained in real time in the prior art.

[0005] The technical solution adopted by the present application is a Czochralski silicon single crystal phase transition surface heat flux solving method, which is implemented according to the following steps:

[0006] Step 1, a temperature sensor is installed at the position of the observation window of the silicon single crystal growth equipment to measure the temperature of the crystal surface measurement point, and a silicon crystal two-dimensional axisymmetric heat conduction control equation and various boundary conditions are established;

[0007] Step 2, the silicon crystal two-dimensional axisymmetric heat conduction control equation established in step 1 is discretized, and the finite difference method is used for iterative solution to obtain a silicon crystal two-dimensional axisymmetric direct heat transfer model with boundary heat flux as input and measurement point temperature as output;

[0008] Step 3, according to the two-dimensional axisymmetric direct heat transfer model of the silicon crystal obtained in step 2, an inverse heat transfer model capable of solving the heat flux of the phase change surface is established, the actual measured crystal surface temperature is taken as the input, and the optimal estimation value of the heat flux at the current time is obtained by continuously optimizing the heat flux information at the current time.

[0009] The application also has the characteristics that,

[0010] Step 1 is specifically implemented according to the following steps:

[0011] Step 1.1, obtain the thermal physical parameters of the silicon crystal, including the crystal surface convective heat transfer coefficient β of the silicon single crystal growth process, the crystal length l, the crystal radius r, the density ρ, the specific heat capacity c, the heat transfer coefficient λ, the Boltzmann constant σ, and the crystal surface property B;

[0012] Step 1.2, establish a three-dimensional heat conduction control equation of the silicon crystal, the expression is:

[0013]

[0014] In the formula, T represents the temperature of any position of the silicon single crystal, x represents the position of the measuring point in the horizontal coordinate, y represents the position of the measuring point in the vertical coordinate, and z represents the position of the measuring point in the axial coordinate;

[0015] Since the silicon crystal is coaxial with the rotating shaft, the three-dimensional heat conduction control equation of the silicon crystal is converted into a three-dimensional heat conduction control equation in the cylindrical coordinate system, the expression is:

[0016]

[0017] In the formula, θ represents the included angle with the radial direction;

[0018] The three-dimensional heat conduction control equation of the silicon crystal in the cylindrical coordinate system is converted into a two-dimensional axisymmetric heat conduction control equation, the expression is:

[0019]

[0020] Step 1.3, establish a thermal boundary condition according to the crystal growth environment, ignore the water cooling effect of the seed crystal end and the shoulder part, the top of the crystal is considered as a plane and is in contact with the external environment, and the external environment temperature is T env Therefore, the top dissipates heat to the external environment through convection and radiation:

[0021]

[0022] At the side of the crystal, the crystal surface is directly in contact with the external environment, so the crystal side also dissipates heat to the external environment through convection and radiation, the expression is:

[0023]

[0024] On the left side of the crystal, the two-dimensional heat conduction control equation is formed by longitudinal sectioning of the silicon crystal, and the silicon crystal is axisymmetric rotation, so the left side of the crystal is an adiabatic boundary, and the expression is:

[0025]

[0026] Below the crystal, the phase transition surface continuously transmits heat to the silicon crystal, so the lower part of the crystal is a heat flux input boundary, and the expression is:

[0027]

[0028] Thus, the two-dimensional axisymmetric heat conduction control equation of the silicon crystal is established.

[0029] Step 2 is implemented according to the following steps:

[0030] Step 2.1, the two-dimensional space calculation domain is discretized by using the finite difference method, and the spatial step in the radial direction is Δr, then the number of grids in the radial direction is Then in the two-dimensional axisymmetric heat conduction control equation of the silicon crystal, is converted to:

[0031]

[0032] In the formula, i is the position in the radial direction after discretization, j is the position in the axial direction after discretization, and k is the position on the time axis after discretization;

[0033] is converted to:

[0034]

[0035] Step 2.2, in the axial direction, the spatial step is Δz, then the number of grids in the axial direction is Then in the two-dimensional heat conduction control equation, is converted to:

[0036]

[0037] Step 2.3, in the time axis, the time step is Δt, then the number of grids in the time axis is Then in the two-dimensional axisymmetric heat conduction control equation of the silicon crystal, is converted to:

[0038]

[0039] Step 2.4, according to step 2.1, 2.2, 2.3, the two-dimensional axisymmetric heat conduction equation of the silicon crystal can be converted into the two-dimensional axisymmetric finite difference heat conduction equation of the silicon crystal, the expression is as follows:

[0040]

[0041] In the formula, the temperature information at the current time can be derived from the temperature information at the last time, and the two-dimensional axisymmetric finite difference heat conduction equation of the silicon crystal is solved by iteration, that is, the temperature change of the silicon crystal at any position under the change of heat flux is obtained. Thus, a two-dimensional axisymmetric direct heat transfer model of the silicon crystal is established, with the boundary heat flux as the input and the temperature information of the silicon single crystal as the output.

[0042] Step 3 is implemented according to the following steps:

[0043] Step 3.1, when the boundary heat flux q n at time n is inverted, the length of the future time step is selected as p, and it is assumed that the heat flow form of the initial future p time steps is the form of constant heat flow, as follows:

[0044] q n =q n+1 =q n+2 =…=q n+p-1 (13)

[0045] q n , q n+1 , q n+2 , and q n+r-1 are the initial boundary heat flow estimation values;

[0046] Step 3.2, select different positions on the surface of the silicon crystal as temperature measurement points, and measure the temperature change at different positions through temperature sensors, and assume that the number of temperature measurement points is M;

[0047] Step 3.3, select the initial heat flow q n at time n as the heat flow estimation result at the last time, if n=1, q n is a random value, then the initial predicted heat flow of the future time step is:

[0048]

[0049] In the formula, rand(n) represents the initial heat flow random value at time n, and the two-dimensional axisymmetric finite difference heat conduction equation of the silicon crystal given in step two is solved by iteration to obtain the temperature field at each time within the future p time steps, and the measurement point temperature estimation value at each time within the future p time steps under the current estimated heat flow is obtained;

[0050] Step 3.4, calculate the sensitivity coefficient H of the mth temperature measuring point in the future p time steps m (i,j,k), the sensitivity coefficient H m (i,j,k) represents the degree of influence of the boundary heat flux on the temperature T(i,j,k) of a certain point of the silicon crystal, the sensitivity coefficient H m (i,j,k) in the time period [t n t n+p-1 ] satisfies the following equation:

[0051]

[0052] H m (i,j,k) = 0, 0≤r≤R = 0, 0≤z≤L, t = t n-1 (16)

[0053]

[0054]

[0055]

[0056]

[0057] wherein H m (i,j,k) represents the sensitivity coefficient with the horizontal coordinate i, the vertical coordinate j and the time step k, H env is the environmental sensitivity, and the initial value is 0;

[0058] The above formula is also discretized by the method described in Step 2, and is solved iteratively to obtain the sensitivity coefficient of the temperature measuring point in the time period [t n t n+r-1 ];

[0059] Step 3.5, optimize the heat flux value at the future p time, expressed as:

[0060]

[0061]

[0062] b represents the iteration number, and Δq is the estimated increment value of the boundary heat flux. Through the above formula, the estimated value of the heat flux in the future p time steps is optimized;

[0063] Step 3.6, establish a convergence criterion,

[0064] Step 3.7, optimize the boundary heat flux in the future p time steps After that, it is judged whether the convergence criterion is reached through step 3.6, if the convergence criterion is reached, the current heat flux estimation value is considered as the accurate value, and the initial temperature of the silicon crystal two-dimensional axisymmetric direct heat transfer model is updated, and the heat flux q of the next time is inverted n+1 If the convergence criterion is not reached, the boundary heat flux in the future p time steps is optimized through step 3.5 Steps 3.5 to 3.6 are repeated until the convergence criterion is reached.

[0065] Step 3.6 is specifically implemented according to the following steps:

[0066] The temperature Y of the measuring point in the future time step is measured i,j The estimated temperature T of the measuring point in the future time step i,j (q) between the deviation between the heat flux of the current time and the accurate value, that is, the convergence is reached, and the expression is:

[0067]

[0068] Wherein, M represents the number of measuring points, T i,j represents the theoretical calculation value of the temperature of the measuring point, Y i,j represents the actual measured value of the temperature of the measuring point, when J(q)≤10 -1 , it is considered that the estimation value of the heat flux of the current time is accurate.

[0069] The beneficial effects of the present application are that the Czochralski silicon single crystal phase transition surface heat flux solving method can solve the phase transition surface heat flux change in real time. First, the temperature of the crystal surface is measured by a temperature sensor, the silicon single crystal forward heat transfer model is solved by using finite difference, and the initial temperature field distribution of the crystal is determined. And by estimating the heat flux of the future finite time step, the estimated value of the temperature change of the measuring point is obtained. Then, the sensitivity coefficient equation is established and solved, and the sensitivity coefficient which can represent the influence degree of the boundary heat flow on the temperature field is obtained. According to the sensitivity coefficient, the heat flux of the future time step is optimized, so that the deviation between the estimated value of the temperature change of the measuring point and the actual temperature measurement value can meet the convergence criterion. Finally, it is verified through experiments that the method proposed is simple and efficient, and can meet the requirement of real-time estimation of the phase transition surface heat flux. In the process of crystal growth, according to the temperature of the measuring point on the crystal surface, the heat flux change process of the phase transition surface is estimated in real time by the method. According to the heat flux change, it is judged whether the crystal growth environment is suitable, so as to improve the control ability of the crystal quality and reduce the generation of structural defects. The method proposed has the characteristics of simplicity, efficiency and real-time, and has high engineering application value. BRIEF DESCRIPTION OF DRAWINGS

[0070] Figure 1 is a temperature measuring point schematic diagram;

[0071] Figure 2 is the inverse result when the boundary heat flux is in the form of a triangular wave;

[0072] Figure 3 is the inverse result when the boundary heat flux is in the form of a triangular wave;

[0073] Figure 4 is the inverse result when the boundary heat flux is in the form of a triangular wave;

[0074] Figure 5 is the inverse result when the boundary heat flux is in the form of a triangular wave;

[0075] Figure 6 is the actual temperature change of measuring point A and the temperature change reconstructed by the estimated heat flux;

[0076] Figure 7 is the actual temperature change of measuring point A and the temperature change reconstructed by the estimated heat flux;

[0077] Figure 8 is the actual temperature change of measuring point A and the temperature change reconstructed by the estimated heat flux. DETAILED DESCRIPTION

[0078] The present application will be described in detail below in combination with the drawings and specific embodiments.

[0079] The present application is a method for solving the inverse problem of heat conduction of a direct silicon single crystal based on a sequence function specification method. In the process of preparing a silicon single crystal by a direct method, the heat transferred to the silicon single crystal cannot be measured, and the heat transferred to the silicon single crystal directly affects the thermal stress inside the silicon single crystal, thereby affecting the crystal quality and the growth process. Therefore, an inverse heat transfer model is established to obtain the heat flux at the phase change surface from the temperature change of the measuring points on the crystal surface. The input of the inverse heat transfer model is the temperature change history of the measuring points on the crystal surface, and the output is the heat flux change process at the phase change surface. First, a two-dimensional axisymmetric direct heat transfer model of a silicon single crystal is established, and the finite difference method is used to solve it. The temperature estimation value of the measuring points on the crystal surface under the estimated heat flux can be obtained. According to the sensitivity of the measuring points and the deviation between the temperature estimation value of the measuring points and the actual temperature value of the measuring points, the heat flux estimation value at the current time is continuously optimized, so as to obtain the optimal estimation value of the heat flux. The heat flux estimation value at the next time is calculated in the same way. Thus, the effect of real-time inversion of the heat flux at the phase change surface is achieved.

[0080] The method for solving the heat flux at the phase change surface of a direct silicon single crystal is implemented according to the following steps:

[0081] Step 1, a temperature sensor is installed at the position of the observation window of the silicon single crystal growth equipment to measure the temperature of the measuring points on the crystal surface, as shown in the accompanying drawings. And the two-dimensional axisymmetric heat conduction control equation of the silicon crystal and the boundary conditions are established; Figure 1

[0082] ​Step 1 is specifically implemented according to the following steps:

[0083] Step 1.1, obtaining the thermal physical parameters of the silicon crystal, including the heat transfer coefficient β of the crystal surface of the silicon single crystal growth process, the crystal length l, the crystal radius r, the density ρ, the specific heat capacity c, the heat transfer coefficient λ, the Boltzmann constant σ, and the crystal surface property B;

[0084] Step 1.2, establishing a three-dimensional heat conduction control equation for the silicon crystal, the expression is:

[0085]

[0086] In the formula, T represents the temperature of any position of the silicon single crystal, x represents the position of the measuring point in the horizontal coordinate, y represents the position of the measuring point in the vertical coordinate, and z represents the position of the measuring point in the axial coordinate;

[0087] Since the silicon crystal is coaxial with the rotation axis, the three-dimensional heat conduction control equation of the silicon crystal is converted into a three-dimensional heat conduction control equation in cylindrical coordinates, the expression is:

[0088]

[0089] In the formula, θ represents the included angle with the radial direction;

[0090] Since the silicon crystal is a three-dimensional cylinder and coaxial with the rotation axis, it is irrelevant to the included angle θ in the heat conduction process, and the three-dimensional heat conduction control equation of the silicon crystal in the cylindrical coordinate system is converted into a two-dimensional axisymmetric heat conduction control equation, the expression is:

[0091]

[0092] Step 1.3, establishing a thermal boundary condition according to the crystal growth environment, ignoring the water cooling effect of the seed crystal end and the shoulder part, the top of the crystal is considered as a plane and is in contact with the external environment, and the temperature of the external environment is T env Therefore, the top dissipates heat to the external environment through convection and radiation:

[0093]

[0094] At the side of the crystal, the crystal surface is directly in contact with the external environment, so the side of the crystal also dissipates heat to the external environment through convection and radiation, the expression is:

[0095]

[0096] On the left side of the crystal, it is a two-dimensional heat conduction control equation formed by longitudinally sectioning the silicon crystal, and the silicon crystal is axisymmetric rotation, so the left side of the crystal is an adiabatic boundary, the expression is:

[0097]

[0098] Under the crystal, the phase transition surface constantly transmits heat to the silicon crystal, so the under the crystal is a heat flux input boundary, and the expression is:

[0099]

[0100] Thus, the two-dimensional axisymmetric heat conduction control equation of the silicon crystal is established.

[0101] Step 2, discretize the two-dimensional axisymmetric heat conduction control equation of the silicon crystal established in step 1, and solve iteratively by using the finite difference method to obtain a two-dimensional axisymmetric direct heat transfer model of the silicon crystal with boundary heat flux as input and measurement point temperature as output;

[0102] Step 2 is implemented according to the following steps:

[0103] Step 2.1, discretize the two-dimensional space calculation domain by using the finite difference method, and set the space step in the radial direction as Δr, so the number of grids in the radial direction is Then in the two-dimensional axisymmetric heat conduction control equation of the silicon crystal, is converted to:

[0104]

[0105] In the formula, i is the position in the radial direction after discretization, j is the position in the axial direction after discretization, and k is the position on the time axis after discretization;

[0106] is converted to:

[0107]

[0108] Step 2.2, set the space step in the axial direction as Δz, so the number of grids in the axial direction is Then in the two-dimensional heat conduction control equation, is converted to:

[0109]

[0110] Step 2.3, set the time step on the time axis as Δt, so the number of grids on the time axis is Then in the two-dimensional axisymmetric heat conduction control equation of the silicon crystal, is converted to:

[0111]

[0112] Step 2.4, according to steps 2.1, 2.2, 2.3, the two-dimensional axisymmetric heat conduction equation of the silicon crystal can be converted into a two-dimensional axisymmetric finite difference heat conduction equation of the silicon crystal, expressed as follows:

[0113]

[0114] In the formula, the temperature information at the current time can be derived from the temperature information at the last time, and the two-dimensional axisymmetric finite difference heat conduction equation of the silicon crystal is solved by iteration, that is, the temperature change of the silicon crystal at any position under any heat flux change is obtained. Thus, a two-dimensional axisymmetric direct heat transfer model of the silicon crystal is established, with the boundary heat flux as the input and the temperature information of the silicon single crystal as the output.

[0115] Step 3, according to the two-dimensional axisymmetric direct heat transfer model of the silicon crystal obtained in step 2, an inverse heat transfer model capable of solving the heat flux of the phase change surface is established. The inverse heat transfer problem is an ill-posed problem, the temperature at the current time is related to not only the heat flux at the last time but also the heat flux before, so the initial temperature field of the silicon crystal at the current time is obtained by solving step two according to the known heat flow information. To solve the ill-posed problem, it is assumed that the future heat flow change is linear, parabolic or constant. The actual measured crystal surface temperature is input, and the current heat flux information is continuously optimized to obtain the optimal estimation value of the current heat flux as output.

[0116] Step 3 is implemented according to the following steps:

[0117] Step 3.1, when the inverse n-time boundary heat flux q n is selected, the length of the future time step is p, and it is assumed that the heat flow form of the initial future p time steps is a constant heat flow form, as follows:

[0118] q n =q n+1 =q n+2 =…=q n+p-1 (36)

[0119] q n , q n+1 , q n+2 , and q n+r-1 are initial boundary heat flow estimation values;

[0120] Step 3.2, select different positions on the surface of the silicon crystal as temperature measurement points, and measure the temperature change process at different positions through temperature sensors, and set the number of temperature measurement points as M;

[0121] Step 3.3, select the initial heat flow q n at time n as the heat flow estimation result at the last time, if n=1, then q nIf the value is random, then the initial predicted heat flux for the future time step is:

[0122]

[0123] In the formula, rand(n) represents the initial random value of heat flux at time n. By iteratively solving the two-dimensional axisymmetric finite difference heat conduction equation of silicon crystal given in step two, the temperature field at each time step in the next p time steps is obtained, and the estimated temperature of the measuring point at each time step in the next p time steps is obtained under the current estimated heat flux.

[0124] Step 3.4: Calculate the sensitivity coefficient H of the m-th temperature measurement point for p future time steps. m (i,j,k), sensitivity coefficient H m (i,j,k) represents the degree of influence of the boundary heat flux on the temperature T(i,j,k) at a certain point in the silicon crystal. Sensitivity coefficient H m (i,j,k) in the time interval [t] n t n+p-1 The following equations are satisfied:

[0125]

[0126] H m (i,j,k)=0,0≤r≤R=0,0≤z≤L,t=t n-1 (39)

[0127]

[0128]

[0129]

[0130]

[0131] In the formula, H m (i,j,k) represents the sensitivity coefficient with x-axis i, y-axis j, and time step k. env Environmental sensitivity, initial value is 0;

[0132] The above equation is also discretized using the method described in step two, and then iteratively solved to obtain the result in [t]. n t n+r-1 Sensitivity coefficient of temperature measurement points within a time period;

[0133] Step 3.5: Optimize the heat flux values ​​at future p time points, expressed as:

[0134]

[0135]

[0136] b represents the iteration number, Δq represents the estimated increment of the boundary heat flux, and the above formula is used to optimize the heat flux estimate over the next p time steps.

[0137] Step 3.6: Establish convergence criteria by measuring temperature Y at measurement points within future time steps. i,j Estimated temperature T at measurement points within future time steps i,j The deviation between (q) determines whether the heat flux at the current moment is an accurate value, i.e., whether convergence has been achieved. The expression is:

[0138]

[0139] Where M represents the number of measuring points, T i,j Y represents the theoretically calculated temperature at the measuring point. i,j This represents the actual measured temperature at the measuring point, when J(q) ≤ 10. -1 If the estimated value of the heat flux at the current moment is considered to be the accurate value, then the estimated value of the heat flux at the current moment is considered to be the accurate value.

[0140] Step 3.7: Optimize the boundary heat flux in each of the next p time steps. Next, step 3.6 determines whether the convergence criterion has been met. If the convergence criterion is met, the current estimated heat flux is considered to be the accurate value, and the initial temperature of the silicon crystal two-dimensional axisymmetric direct heat transfer model is updated to invert the heat flux q at the next moment. n+1 If the convergence criterion is not met, then step 3.5 is used to optimize the boundary heat flux in the next p time steps. Repeat steps 3.5 to 3.6 until the convergence criterion is met; the flowchart for solving the heat flux change process of the phase transition surface is attached. Figure 2 As shown.

[0141] Thus, steps 1 to 3 have established a method for calculating the heat flux of the Czochralski silicon single crystal phase transformation surface based on the temperature change at the measuring point.

[0142] To verify the accuracy and practicality of the method proposed in this invention, triangular and rectangular wave heat fluxes were used as boundary heat fluxes to obtain the temperature changes at the measuring points under two different waveforms. The obtained temperature changes at the measuring points were then used as experimental data to invert the heat flux at the phase transition surface. (See attached diagram.) Figure 3 Appendix Figure 4As shown in the figure, the horizontal axis represents time, the vertical axis represents the heat flux value of the phase transition surface, the dashed line represents the actual boundary heat flux, and the marked line represents the estimated boundary heat flux value. The figure shows that the heat flux of the phase transition surface can be accurately retrieved, and the estimated heat flux variation trend is consistent with the actual heat flux variation trend, demonstrating the accuracy of the algorithm. To further verify the practicality of the method, temperature changes on the crystal surface during the heating stage were collected using a temperature sensor as experimental data to retrieve the heat flux of the phase transition surface. The results are shown in the attached figure. Figure 5 As shown, the heat flux of the phase change surface exhibits a slow increasing trend during the heating stage, and this increasing trend gradually levels off. To verify the accuracy of the inversion results, the estimated heat flux value was used as input to a two-dimensional axisymmetric direct heat transfer model, yielding temperature estimates at three different measuring points, which were then compared with the actual values. (See attached figure.) Figure 6 Appendix Figure 7 Appendix Figure 8 As shown, the estimated temperature values ​​at the measuring points and the actual values ​​at the measuring points exhibit the same trend, indicating the accuracy of the inverted phase change surface heat flux results.

Claims

1. A method for solving the heat flux of the phase transformation surface of a Czochralski-grown silicon single crystal, characterized in that, The specific steps are as follows: Step 1: Install a temperature sensor at the observation window of the silicon single crystal growth equipment to measure the temperature at the measuring point on the crystal surface, and establish the two-dimensional axisymmetric heat conduction control equation of the silicon crystal and various boundary conditions. Step 2: Discretize the two-dimensional axisymmetric heat conduction control equation of silicon crystal established in Step 1, and solve iteratively using the finite difference method to obtain a two-dimensional axisymmetric direct heat transfer model of silicon crystal with boundary heat flux as input and measuring point temperature as output. Step 3: Based on the two-dimensional axisymmetric direct heat transfer model of silicon crystal obtained in Step 2, establish an inverse heat transfer model that can solve the heat flux of the phase change surface. Take the actual measured crystal surface temperature as input, and continuously optimize the heat flux information at the current moment to obtain the output as the optimal estimate of the heat flux at the current moment. Step 3 is implemented in the following steps: Step 3.1, in the inversion Time boundary heat flux When choosing the length of the future time step, it is... Assuming the initial future The heat flow at each time step is a constant heat flow, as shown in the following equation: (13) , , Until This is the initial boundary heat flux estimate; Step 3.2: Select different locations on the silicon crystal surface as temperature measurement points, and measure the temperature change process at different locations using a temperature sensor. Let the number of temperature measurement points be... ; Step 3.3, Select Initial heat flow at time If the heat flow estimation result from the previous moment is true, then... ,but If the value is random, then the initial predicted heat flux for the future time step is: (14) In the formula, express The initial random value of heat flux is determined, and the future value is obtained by iteratively solving the two-dimensional axisymmetric finite difference heat conduction equation of silicon crystal given in step two. The temperature field at each moment within a time step is obtained, and the future temperature field is calculated under the current estimated heat flux. The estimated temperature of the measuring point at each moment within a time step; Step 3.4, Calculate the future The time step, the first Sensitivity coefficient of each temperature measuring point Sensitivity coefficient This represents the boundary heat flux with respect to the temperature at a certain point in the silicon crystal. The extent of the impact Sensitivity coefficient In time period Satisfy the following equation: (15) (16) (17) (18) (19) (20) In the formula, The convective heat transfer coefficient at the crystal surface during the silicon single crystal growth process. For the crystal surface density during silicon single crystal growth process, This refers to the specific heat capacity of the crystal surface during the silicon single crystal growth process. Where is the crystal radius, The heat transfer coefficient is... For Boltzmann constant, For crystal length, As a crystal surface property, The x-coordinate is The vertical axis is The time step is The sensitivity coefficient, Environmental sensitivity, initial value is 0; The above formula is also discretized using the method described in step 2, and then iteratively solved to obtain the result. Sensitivity coefficient of temperature measurement points within a time period; Step 3.5, Optimize the Future The heat flux value at time is expressed as: (22) (23) Indicates the number of iterations. The estimated increment of the boundary heat flux, obtained through the above formula, optimizes the future... Estimated heat flux over a time step; Step 3.6: Establish convergence criteria. Step 3.7, in each optimization of the future Boundary heat flux within a time step Next, step 3.6 determines whether the convergence criterion has been met. If the convergence criterion is met, the current estimated heat flux is considered to be the accurate value, and the initial temperature of the silicon crystal two-dimensional axisymmetric direct heat transfer model is updated to invert the heat flux at the next moment. If the convergence criterion is not met, then optimize the future through step 3.

5. Boundary heat flux within a time step Repeat steps 3.5 to 3.6 until the convergence criterion is met.

2. The method for solving the heat flux of the phase transformation surface of a Czochralski-grown silicon single crystal according to claim 1, characterized in that, Step 1 is implemented in the following steps: Step 1.1: Obtain the thermal properties of silicon crystals, including the convective heat transfer coefficient at the crystal surface during the silicon single crystal growth process. Crystal length Crystal radius ,density Specific heat capacity Heat transfer coefficient Boltzmann constant Crystal surface properties ; Step 1.2: Establish the three-dimensional heat conduction control equation for silicon crystal, the expression of which is: (1) In the formula, This represents the temperature at any location in a silicon single crystal. Indicates the position of the measuring point on the horizontal axis. Indicates the position of the measuring point on the vertical axis. Indicates the position of the measuring point in axial coordinates; Since the silicon crystal is coaxial with the rotation axis, the three-dimensional heat conduction control equation of the silicon crystal is transformed into the three-dimensional heat conduction control equation in cylindrical coordinates, expressed as: (2) In the formula, Indicates the angle with respect to the radial direction; The three-dimensional heat conduction control equations for silicon crystals in cylindrical coordinates are transformed into two-dimensional axisymmetric heat conduction control equations, expressed as follows: (3) Step 1.3: Establish thermal boundary conditions based on the crystal growth environment. Ignoring the water-cooling effect and shoulder formation at the seed crystal end, the top of the crystal is considered a plane in contact with the external environment, where the temperature is... Therefore, the top dissipates heat to the external environment through convection and radiation: (4) At the side of the crystal, the crystal surface is in direct contact with the external environment, so the side of the crystal also dissipates heat to the external environment through convection and radiation, as expressed in the following expression: (5) On the left side of the crystal, since the two-dimensional heat conduction governing equation is formed by longitudinally sectioning the silicon crystal, and the silicon crystal is axisymmetrically rotated, the left side of the crystal is an adiabatic boundary, expressed as: (6) Beneath the crystal, the phase transition surface continuously transfers heat to the silicon crystal, therefore the area beneath the crystal is the heat flux input boundary, expressed as: (7) Thus, the two-dimensional axisymmetric thermal conduction control equation for silicon crystal was established.

3. The method for solving the heat flux of the phase transition surface of a Czochralski silicon single crystal according to claim 2, characterized in that, Step 2 is implemented in the following steps: Step 2.1: Discretize the two-dimensional computational domain using the finite difference method, setting the spatial step size in the radial direction as . Then the number of grids in the radial direction In the two-dimensional axisymmetric heat conduction governing equation of silicon crystal, Transform into: (8) In the formula, This represents the position along the radial direction after discretization. This represents the axial position after discretization. This represents the position on the time axis after discretization. Transform into: (9) Step 2.2: In the axial direction, set the spatial step size to... The number of grids in the axial direction In the two-dimensional heat conduction governing equations, Transform into: (10) Step 2.3: On the timeline, set the time step to... The number of grids on the time axis In the two-dimensional axisymmetric thermal conduction governing equation of silicon crystal, Transform into: (11) Step 2.4: Based on steps 2.1, 2.2, and 2.3, the two-dimensional axisymmetric heat conduction governing equation of silicon crystal can be transformed into the two-dimensional axisymmetric finite difference heat conduction equation of silicon crystal, as shown in the following expression: (12) In the formula, the temperature information at the current moment can be derived from the temperature information at the previous moment. By iteratively solving the two-dimensional axisymmetric finite difference heat conduction equation of silicon crystal, the temperature change at any position of silicon crystal under any heat flux change is obtained. Thus, a two-dimensional axisymmetric direct heat transfer model of silicon crystal with boundary heat flux as input and silicon single crystal temperature information as output is established.

4. The method for solving the heat flux of the phase transition surface of a Czochralski silicon single crystal according to claim 3, characterized in that, Step 3.6 is implemented in the following steps: Temperature measured at points within a future time step Temperature estimated by measuring points within future time steps The deviation between the two values ​​determines whether the heat flux at the current moment is an accurate value, i.e., convergence has been achieved. The expression is: (21) in, Indicates the number of measuring points. This represents the theoretically calculated temperature value at the measuring point. This indicates the actual measured temperature value at the measuring point. If the estimated value of the heat flux at the current moment is considered to be the accurate value, then the estimated value of the heat flux at the current moment is considered to be accurate.

Images