CN102521886B - Three-dimensional simulation method for chemical vapor deposition process - Google Patents

Abstract

The invention relates to a three-dimensional simulation method of a chemical vapor deposition process, which belongs to the field of deposition process simulation in microelectronic processing. Divide into several triangular planes; use the Monte Carlo (MC) method combined with the spatial octree non-uniform subdivision technology to quickly track the transport process of sedimentary particles, calculate the triangular planes that reach the surface of the depositional graphics, and combine According to the corresponding deposition model and the area of the triangular plane, etc., the deposition rate of each triangular plane on the surface of the deposition figure is calculated, which avoids the blindness of the calculation of the deposition rate; the deposition rate is obtained according to the calculation, and the chemical vapor deposition (CVD) is realized by the level set method (LS). ) process for accurate tracking of the motion of the deposited surface. The method can simulate various vapor deposition processes such as low-pressure chemical vapor deposition, plasma-enhanced chemical vapor deposition, and atmospheric pressure chemical vapor deposition, and realizes deposition simulation of deposition patterns composed of complex and various materials.

Description

A 3D Simulation Method of Chemical Vapor Deposition Process

technical field

The invention belongs to the field of deposition process simulation in microelectronic processing, and provides a method of expressing the distribution of Ω materials in the simulation area by using the cell method (CM), realizing the transportation of deposited particles by using the Monte Carlo method (MC), and using the level set method (LS) Simulation method for chemical vapor deposition (CVD) process to achieve surface motion.

Background technique

At present, integrated circuits have been applied to various fields, and its industry has become a basic part of future economic development, and global competition is becoming increasingly fierce. These undoubtedly bring new challenges to the development of large-scale integrated circuits. With the continuous shrinking of the feature size of integrated circuits and the continuous increase of semiconductor wafers, the cornerstone of its production, the integration of integrated circuits is getting higher and higher, which poses new challenges to the processing technology of integrated circuits. In order to increase the integration level of integrated circuits, new requirements are placed on etching, one of which is a higher aspect ratio. However, the increase of the aspect ratio brings greater challenges to the deposition, the factors affecting the deposition quality increase, and the influence of the shadow effect increases.

In order to improve the deposition quality, it is necessary to deeply understand the mechanism of the deposition process and process, and use it to guide the design of deposition process parameters and equipment. The current research on this aspect is mainly realized through trial and error, that is, to find out the appropriate parameters through a large number of deposition experiments. On the one hand, this method is time-consuming, laborious and costly, and on the other hand, the appropriate parameters obtained by experiments are not optimal. Another method is to realize it by simulation, that is, to predict the deposition morphology on the characteristic pattern through certain processing parameters. Compared with the trial and error method, it can reduce the test cost and find out the main optimization parameters. Therefore, the simulation of the deposition process is an effective method to recognize and understand the mechanism of deposition processing. At present, the simulation methods of deposition process are mainly line algorithm, cell method and level set method, etc., and have been developed from two-dimensional to three-dimensional to realize the simulation of deposition process with more complex graphics.

The line algorithm is mainly suitable for two-dimensional simulation. When the surface morphology of the deposition pattern changes smoothly, the interface can be accurately tracked. However, when the time step is large or the surface is undulating, it is easy to appear unstable or inaccurate motion interface. It is difficult to express that the deposition process affects the physical and chemical properties of the deposition and the surface material; the level set method is an effective method for accurately and robustly tracking the complex motion interface, suitable for two-dimensional and three-dimensional simulations, and can accurately track the situation of severe surface changes, However, it cannot represent the physical and chemical properties and surface materials in the deposition process, and is not suitable for the interface evolution of various material simulation graphics; the cell representation method can represent materials and physical chemistry at any position by dividing the simulation area Ω into cells properties, but cannot achieve accurate tracking of motion interfaces. The deposition process is not only affected by processing conditions and processes, but also by the material at the deposition location. Current simulation tools still cannot meet the simulation requirements for complex multi-material deposition patterns.

Contents of the invention

The purpose of the present invention is to overcome the deficiencies of the prior art, and propose a three-dimensional simulation method of chemical vapor deposition process, which has the advantages of simulating low-pressure chemical vapor deposition, plasma enhanced chemical vapor deposition, chemical vapor deposition under atmospheric pressure, etc. A variety of vapor deposition processes can realize the deposition simulation of complex deposition patterns composed of various materials.

A three-dimensional simulation method of a chemical vapor deposition process proposed by the present invention is characterized in that it comprises the following steps:

1) Obtain initial parameters: According to the actual processing equipment and deposition gas, estimate and calculate the boundary parameters of the electric field, the type, quantity and velocity distribution of particles participating in the deposition; input the simulation area and deposition graphics; set the cell side length; set the simulation time T ;Set the maximum forward distance;

2) Divide the simulation area into multiple cells, each cell corresponds to the material at the position of the cell, and obtain the corresponding material deposition model; triangulate the surface of the deposition graph, and divide it into multiple triangular planes , to calculate the deposition velocity of each triangular plane to obtain the motion of the entire deposition pattern;

3) Use the level set function to represent the shortest distance from each point in the simulation area to the deposition pattern, and initialize the level set function so that it satisfies formula (1)

$\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; d</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>& \stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; \&Element;</span>{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>& \stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>\\ <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; d</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>& \stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; \&Element;</span>{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

In formula (1) Indicates the 0 level set at time t=0, that is, the deposition pattern, Ω represents the simulation area, Ω ₁ (t) represents the gas, Ω ₂ (t) represents the material; Ω=Ω ₁ (t)∪Ω ₂ (t)∪Γ (t); Represents a point in the simulation domain, in Cartesian coordinates represents a point in the simulation region Ω The distance to the deposition pattern Γ(t);

4) Use the spatial octree non-uniform subdivision technology to divide the simulation area into grids. Each grid is composed of multiple cells, so as to quickly track the motion trajectory of the deposition particles and obtain the contact between the deposition particles and the surface of the deposition graph. Position; according to the type, quantity and velocity distribution of sedimentary particles, use the Monte Carlo MC method to randomly sample the sedimentary particles to obtain the sedimentary particles to be simulated; according to the type of sedimentary particles to be simulated, use the spatial octree technology to realize fast tracking to be simulated According to the transportation process of the deposited particles, the parameters of the deposited particles arriving at the deposition position are calculated; according to the materials at the deposition position and the parameters of the deposition particles arriving at the deposition position, the corresponding deposition model is selected, and the deposition amount of the triangular plane of the deposition position is calculated and updated; The specified number of deposited particles turns to step 5);

5) Calculate the deposition rate of each triangular plane according to the deposition amount, area and surface unit normal vector of each triangular plane; calculate the stepping time interval ΔT according to the deposition rate of the triangular plane and the maximum forward distance; according to the deposition rate of each triangular plane, use the value The calculation method solves the level set governing equation to realize the movement of the surface of the depositional graph; reinitializes the level set function; re-triangulates the surface of the current depositional graph, and updates the information in the cell according to the current depositional graph; Cells are replaced with freshly deposited material;

6) After updating the simulation time T=T-ΔT, if T=0, go to step 7), otherwise go to step 4);

7) Output the simulation results, generate and display the final 3D depositional graphics according to all the triangular planes of the depositional graphics.

Features and beneficial effects of the present invention:

The invention provides a method of expressing the Ω material distribution in the simulation area by using the cell method (CM), realizing the transport of particles by using the Monte Carlo (MC) method, and combining the deposition models of different materials to calculate the Ω at each position on the surface of the deposition graph The amount of deposition avoids the blindness of the calculation of the deposition rate; according to the obtained deposition rate, the surface evolution of the chemical vapor deposition (CVD) process is realized by the level set method (LS). The invention can simulate various vapor deposition processes such as low-pressure chemical vapor deposition, plasma enhanced chemical vapor deposition, atmospheric pressure chemical vapor deposition, etc., realizes deposition simulation of deposition patterns composed of complex and various materials, and solves the problem of single constituent materials in previous simulations , Inaccurate simulation and other issues.

Description of drawings

Fig. 1 is the program flowchart of the inventive method;

Figure 2 is a two-dimensional representation of the vertical section of the simulation region Ω;

Figure 3 shows the Cartesian coordinate system where the system is located.

Detailed ways

The present invention proposes a three-dimensional simulation method of a chemical vapor deposition process, which is described in detail in conjunction with the accompanying drawings and implementation as follows:

The general flow of the method of the present invention is as shown in Figure 1, and the specific implementation method comprises the following steps:

1) Acquire initial parameters: According to the actual processing equipment and deposition gas, estimate and calculate the electric field boundary parameter E _c , the type S, number N _s and velocity distribution P _v of particles participating in the deposition; input the simulation area Ω and the deposition graph Γ(t); Set the cell side length l; set the simulation time T; set the maximum forward distance dist;

2) Divide the simulated area Ω into cells, divide it into multiple cells, each cell corresponds to the material at the cell position, and obtain the corresponding material deposition model; triangulate the surface of the deposition graph Γ(t), Divide into multiple triangular planes to calculate the deposition velocity of each triangular plane to obtain the movement of the entire depositional figure Γ; specifically include:

2.1) Divide the simulation area Ω into cells: Divide the simulation area Ω into several cells (cubes) with a side length of l, which are used to store the material of each cell position in the simulation area Ω; for example, the current state of a cell is : state=0, indicating that the cell position is the gas phase region; state=1, indicating that the cell position is the first material; state=2, indicating that the cell position is the second material (each material is in the deposition process They have different physical and chemical properties, which are manifested as different deposition models; it is precisely because the reaction gas has different deposition models for different materials, when calculating the deposition amount, it is necessary to select the deposition position (the contact position between the deposition particles and the surface of the deposition pattern) have a material deposition model corresponding to the material);

2.2) Triangulate the surface of the deposition pattern Γ(t): According to the characteristic scale requirement of the deposition pattern Γ(t) to be simulated (this is the minimum line width in the semiconductor processing process), using the surface triangulation algorithm, Divide the surface of the depositional figure Γ(t) to be simulated into multiple triangular planes, and obtain the movement of the entire depositional figure Γ(t) to be simulated by calculating the deposition velocity of each triangular plane;

3) Use the level set function to represent the shortest distance from each point in the simulation area Ω to the deposition pattern Γ(t), and initialize the level set function so that it satisfies formula (1)

$\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; d</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>& \stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; \&Element;</span>{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>& \stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>\\ <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; d</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>& \stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; \&Element;</span>{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span>\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

In formula (1) Indicates the 0 level set at time t=0 (that is, the boundary of two materials, that is, the deposition pattern), the vertical section of the simulation area Ω is shown in Figure 2, the white area 1 in the figure represents the gas Ω ₁ (t), and the oblique line area 2 represents the material Ω ₂ (t), and the dotted line represents the profile; Ω=Ω ₁ (t)∪Ω ₂ (t)∪Γ(t); Represents a point in the simulation domain, in Cartesian coordinates represents a point in the simulation region Ω The distance to the deposition pattern Γ(t);

4) Use the spatial octree non-uniform subdivision technology to divide the simulation area Ω into grids. Each grid is composed of multiple cells, so as to quickly track the trajectory of the deposition particles, and obtain the relationship between the deposition particles and the surface of the deposition graph. Contact position, if there is plasma in the simulation area Ω, according to the boundary conditions of the electric field in the simulation area Ω, use the existing finite difference method to calculate the plasma formation potential in the simulation area Ω; according to the deposition particle type S, number Ns and velocity Distribution Pv, using the Monte Carlo MC method to randomly sample the sediment particles to obtain the sediment particles to be simulated According to the sediment particles to be simulated type, using the spatial octree technology to quickly track the transport process of the deposition particles to be simulated, and calculate the parameters of the deposition particles arriving at the deposition location Depending on the material at the deposition location and the parameters of the deposited particles reaching the deposition location Select the corresponding deposition model, calculate and update the deposition amount of the triangular plane of the deposition position; reach the specified number of deposition particles and turn to step 5); specifically include:

4.1) Use the spatial octree non-uniform subdivision technology to divide the simulation area Ω into a grid, so as to realize the fast tracking of the particle trajectory and obtain the contact position between the particle and the surface of the deposition pattern. The grid division method is to divide the spatial cuboid containing the simulation area Ω into 8 identical sub-cuboid grids according to the three coordinate axes of the Cartesian coordinate system (z-axis direction is downward) as shown in Figure 3, and organize them into an octagonal tree. Tree. If the number of triangular planes contained in a sub-cuboid grid is greater than the given threshold ε _Δ (generally ε _Δ = 4), then the sub-cuboid grid is further subdivided until each grid contains triangles The number of planes is less than or equal to the given threshold ε _Δ ; if there is plasma in the simulation region Ω, according to the electric field boundary condition E _c of the simulation region Ω, using the existing finite difference method, calculate the plasma in the simulation region Ω The electric potential is formed, and the Poisson equation that the electric field satisfies is ▽ ² Φ=-ρ/ε ₀ obtained from the Maxwell equation of the electromagnetic field theory. This ▽ ² represents the Laplacian operator, Φ represents the electric potential (in volts), and ρ is the electric charge Bulk density (in coulomb/m3), and _ε0 is the vacuum permittivity (in farad/m); according to the condition that the sum of positive charge and electron charge in the local simulation region ψ is close to zero in steady state, the electric field can be simplified as Laplace equation ▽ ² Φ = 0 to solve, so that the existing finite difference method can be used in combination with boundary conditions to obtain the potential of each point in the simulation area Ω; in addition to the triangular plane in each cuboid grid, Also preserves the electric field formed by the plasma, which is used to simulate the movement of charged ions;

4.2) According to the sedimentary particle type S, the particle number N _s and the sedimentary particle velocity distribution P _v , use the Monte Carlo MC method to randomly sample the sedimentary particles to obtain the sedimentary particles to be simulated Among them, i represents the serial number of deposited particles; Indicates deposited particles The position of , in the Cartesian coordinate system, Initial z=0, which means the upper boundary of the gas phase in the simulation area Ω, x, y are randomly generated according to the size of the simulation area Ω; _si indicates the type of deposited particles, which is set by the deposited particle components; represent The initial velocity of , given by the isotropic Maxwell velocity distribution $f\left(v_i\right)={\left(\frac{m_i}{2\mathrm{\&pi;k}{T}_i}\right)}^{3/2}{e}^{-\frac{m_i{\left|{\stackrel{\&Right\; Arrow;}{v}}_i\right|}^2}{2k{T}_i}}$ produce, ${\stackrel{\&Right\; Arrow;}{v}}_i=(v_{\mathrm{ix}},v_{\mathrm{iy}},v_{\mathrm{iz}}),v_{\mathrm{ix}},v_{\mathrm{iy}},v_{\mathrm{iz}}$ are the component velocities in the three directions of Cartesian coordinate system x, y, and z, and -∞<v _ix , v _iy <+∞, v _iz >0, M _i is the mass of deposited particles, and k is Boltzmann's constant , T _i is the deposition particle temperature; the deposition particle type includes all neutral active molecules, atoms and ionized ions participating in the action (such as when using H ₂ , N ₂ and CH ₄ plasma-enhanced chemical vapor deposition of Si ₃ N ₄ , The deposition gas contains a large number of N ₂ ⁺ ions, N atoms and CN particles);

4.3) According to the deposited particles type, using the space octree technology to quickly track the transport process of the sediment particles, and calculate the parameters of the sediment particles arriving at the deposition position

If s _i ∈ ion, the deposited particle The movement is affected by the electric field E formed by the plasma, and the movement caused by acceleration (due to the inhomogeneous electric field density, It is not a uniformly accelerated linear motion. For simplicity, this motion can be regarded as consisting of a uniformly accelerated linear motion with a time interval of Δt. Therefore, The motion trajectory is a polyline instead of a straight line); the specific calculation is as follows:

According to the electric field E in the simulation area Ω, Newton's law of motion and the current position of the deposited particles and speed Compute the new position after the Δt time interval and new speed then by and Determine the line segment AB, use the space octree technology to determine whether this line segment intersects with the surface of the depositional graph, if not intersect and this line segment does not leave the simulation area Ω, then set the new position and new speed as current location and speed Repeat this process until it intersects with the surface of the deposition pattern, and find the triangular plane _Δi intersecting with the deposition pattern; then calculate the deposition particles angle θ _i between the incident direction and the normal of the triangular plane; and use the surface reflection model (that is, compare the value of θ _i with the reflection threshold θ ₀ , if θ>θ ₀ , reflect the deposited particles according to the law of specular reflection, otherwise the deposition Particles are deposited at the position of the cell) to calculate reflection or deposition, that is, if θ _i > θ ₀ (θ ₀ > π/3), calculate according to the principle of reflection new location for new speed as current location and speed Go to step 4.3) to continue; otherwise go to step 4.4);

If s _i ∈ neutral active particles, then The motion is constant speed Direction-invariant linear motion, using space octree technology to quickly calculate the triangular plane _Δi where this line intersects with the deposition pattern Γ(t); then according to the adsorption probability model (incident particles with adsorption probability S _c (0 < _S <1) Decide whether the particle is adsorption or scattering, 1- _Sc probability particle scattering, the reflection angle has nothing to do with the incident angle, satisfying the diffuse reflection distribution) Decide whether the particle is scattering or adsorption, that is, generate a random number p (0≤p≤1) , if p>S _c is scattering, randomly generate the scattering angle and velocity and calculate the current position and speed Go to step 4.3) to continue; otherwise go to step 4.4);

4.4) According to the material of the deposition location and the parameters of the deposition particles reaching the deposition location Select a specific deposition model, calculate and update the deposition amount of the triangular plane of the deposition location:

Taking the plasma-enhanced vapor deposition process as an example, the corresponding deposition model can be expressed as formula (2):

$\mathrm{<span\; class="notranslate">\; D.</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}}{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; ion</span>}}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; v</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&theta;</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; state</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; ion</span>}}\mathrm{<span\; class="notranslate">\; ds</span>}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}}{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; state</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}\mathrm{<span\; class="notranslate">\; ds</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

in: Indicates the plasma-enhanced deposition coefficient, which is generally obtained through experiments, and it is related to the particle velocity The angle θ between the plane normal of the particle movement direction is related to the state of the deposition material; K _d (state) represents the deposition coefficient of neutral particles, which is related to the material state, and is also obtained through experiments; F _ion and F _d represent the motion to deposition figure Γ (t) Ion flux and neutral particle flux on the upper triangular plane ds; θ represents the angle between the deposition particles and the normal direction of the triangular plane;

If s _i ∈ ion, then ${\mathrm{D.}}_{{\mathrm{\&Delta;}}_i}={\mathrm{D.}}_{{\mathrm{\&Delta;}}_i}+K_{\mathrm{ion}}({\stackrel{\&Right\; Arrow;}{v}}_i,{\mathrm{\&theta;}}_i,\mathrm{state})*f_{\mathrm{ion}}^i;$

If s _i ∈ neutral active particles, then ${\mathrm{D.}}_{{\mathrm{\&Delta;}}_i}={\mathrm{D.}}_{{\mathrm{\&Delta;}}_i}+K_d\left(\mathrm{state}\right)*f_d^i;$

in: Indicates the deposition amount on the triangular plane _Δi , express Plasma-enhanced deposition coefficient; K _d (state) represents neutral particles deposition coefficient; and respectively represent the current particle Representative ion flux and neutral particle flux;

4.5) Repeat steps 4.2)-4.4) to reach the specified number of particles and turn to step 5);

5) According to the deposition amount of each triangular plane area and the unit normal vector to calculate the deposition velocity of each triangular plane Calculate the step time interval ΔT according to the deposition speed of the triangular plane and the maximum forward distance dist (dist>0); according to the deposition speed of each triangular plane Use the numerical calculation method to solve the level set control equation to realize the movement of the surface of the deposition figure; reinitialize the level set function; re-triangulate the surface of the current deposition figure, and update the information in the cell according to the current deposition figure Γ(t), Cells whose state has changed are replaced with newly deposited materials; specifically include:

5.1) According to the deposition amount of each triangular plane area and the unit normal vector to calculate the deposition velocity of each triangular plane

${\mathrm{<span\; class="notranslate">\; R</span>}}_{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; D.</span>}}_{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; S</span>}}_{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span>{\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; no</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

5.2) Calculate the stepping time interval ΔT according to the deposition speed of the triangular plane and the maximum forward distance dist (dist>0);

$\mathrm{<span\; class="notranslate">\; \&Delta;T</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; dist</span>}}{\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; ls</span>}}{\mathrm{<span\; class="notranslate">\; max</span>}}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; R</span>}}_{{\mathrm{<span\; class="notranslate">\; \&Delta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$

Among them: dist represents the maximum advancing distance, and ls represents the triangular plane set that forms the sedimentary figure Γ(t);

5.3) According to the deposition rate of each triangular plane Using numerical calculation methods to solve the level set governing equation (formula (5)) to realize the movement of the surface;

φ _t + V·▽φ＝0 (5)

Among them: φ _t is the partial derivative of the level set function φ to time t, ▽φ is the gradient of φ at (x, y, z), V represents the velocity at each point (x, y, z) of the domain, approximately The deposition velocity of the triangular plane where it is located When using the numerical calculation method to solve the problem, it is necessary to discretize the domain of the level set function (simulation area Ω); since the entire simulation area Ω and each position velocity (velocity field) are used when solving the level set governing equation, and in During the deposition process, we can only calculate the velocity of each point of the deposition graph Γ(t), therefore, the velocity field of the entire simulation area Ω must be calculated by the surface velocity. In order to improve efficiency, the narrow-band level set method is used here, and only the velocity of the position near the deposition pattern Γ(t) is calculated (the velocity of the position near the deposition pattern Γ(t) is equal to the deposition of the triangular plane on the nearest deposition pattern Γ(t) speed);

5.4) Re-initialize the level set function according to formula (6) so that it satisfies formula (1):

$\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; sign</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

Among them: φ ₀ represents the current level set function, is the level set function to be calculated, and its steady-state solution is a signed distance function satisfying formula (1), sign(φ ₀ ) sign function, for the convenience of solving, it is smoothed into formula (7), ε>0 is very small positive number of

$\mathrm{<span\; class="notranslate">\; sign</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}{\sqrt{{{\mathrm{<span\; class="notranslate">\; \&phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

5.5) Re-triangulate the surface of the current depositional graph, and update the information in the cells according to the current depositional graph Γ(t), and replace the cells whose state has changed with new deposition materials;

6) After updating the simulation time T=T-ΔT, if T=0, go to step 7), otherwise go to step 4);

7) Output the simulation results, generate and display the final 3D depositional graphics according to all the triangular planes of the depositional graphics.

Claims (4)

1. a three-dimensional simulation method of chemical vapor deposition process, is characterized in that, comprises the following steps: 1) Obtain initial parameters: According to the actual processing equipment and deposition gas, estimate the electric field boundary parameters, the type, quantity and velocity distribution of particles participating in the deposition; input the simulation area and deposition graphics; set the cell side length; set the simulation time T; Set the maximum forward distance; 2) Divide the simulation area into multiple cells, each cell corresponds to the material at the cell position, and obtain the corresponding material deposition model; triangulate the surface of the deposition graph, and divide it into multiple triangular planes , to calculate the deposition velocity of each triangular plane to obtain the movement of the entire deposition pattern; 3) Use the level set function to represent the shortest distance from each point in the simulation area to the depositional graph, and initialize the level set function so that it satisfies formula (1) In formula (1) Indicates the 0 level set at time t=0, that is, the deposition pattern, Ω represents the simulation area, Ω ₁ (t) represents the gas, Ω ₂ (t) represents the material; Ω=Ω ₁ (t)∪Ω ₂ (t)∪Γ (t); Represents a point in the simulation domain, in Cartesian coordinates represents a point in the simulation region Ω The distance to the deposition pattern Γ(t); 4) Use the space octree non-uniform subdivision technology to divide the simulation area into grids, each grid is composed of multiple cells, so as to realize the fast tracking of the trajectory of the deposition particles, and obtain the contact between the deposition particles and the surface of the deposition graphics Position; according to the type, quantity and velocity distribution of sedimentary particles, use the Monte Carlo MC method to randomly sample the sedimentary particles to obtain the sedimentary particles to be simulated; according to the type of sedimentary particles to be simulated, use the spatial octree technology to realize fast tracking to be simulated According to the transportation process of the deposited particles, the parameters of the deposited particles arriving at the deposition position are calculated; according to the materials at the deposition position and the parameters of the deposition particles arriving at the deposition position, the corresponding deposition model is selected, and the deposition amount of the triangular plane of the deposition position is calculated and updated; The specified number of deposited particles turns to step 5); 5) Calculate the deposition rate of each triangular plane according to the deposition amount, area and surface unit normal vector of each triangular plane; calculate the stepping time interval ΔT according to the deposition rate of the triangular plane and the maximum advancing distance; according to the deposition rate of each triangular plane, use the value The calculation method solves the level set governing equation to realize the movement of the surface of the depositional graph; reinitializes the level set function; re-triangulates the surface of the current depositional graph, and updates the information in the cell according to the current depositional graph; Cells are replaced with new deposited material; 6) After updating the simulation time T=T-ΔT, if T=0, go to step 7), otherwise go to step 4); 7) Output the simulation results, generate and display the final three-dimensional depositional graphics according to all the triangular planes of the depositional graphics. the 2. method as claimed in claim 1, is characterized in that, described step 2) specifically comprises: 2.1) Divide the simulation area into cells: Divide the simulation area Ω into several cells with a side length of l, which are used to store the materials of each cell position in the simulation area Ω; The material deposition model corresponding to the material, the deposition position is the contact position between the deposition particles and the deposition pattern surface; 2.2) Triangulate the surface of the sedimentary figure: According to the characteristic scale requirements of the sedimentary figure Γ(t), use the surface triangulation algorithm to divide the surface of the sedimentary figure Γ(t) into multiple triangular planes, and calculate the The deposition velocity of , to obtain the movement of the entire deposition pattern Γ(t). the 3. method as claimed in claim 1, is characterized in that, described step 4) specifically comprises: 4.1) Use the spatial octree non-uniform subdivision technology to divide the simulation area into grids to realize fast tracking of the particle trajectory and obtain the contact position between the particles and the surface of the deposition graph; The spatial cuboid in the Cartesian coordinate system is divided into 8 identical sub-cuboid grids according to the three coordinate axes of the Cartesian coordinate system, and is organized into an octree; if the number of triangular planes contained in a sub-cuboid grid is greater than the given threshold ε _Δ , then the sub-cuboid grid is further subdivided until the number of triangular planes contained in each grid is less than or equal to the given threshold ε _Δ ; if there is plasma in the simulation area Ω, according to the simulation area Ω For the electric field boundary condition E _c , use the existing finite difference method to calculate the electric potential formed by the plasma in the simulated region Ω. From the Maxwell equation of the electromagnetic field theory, the Poisson equation that the electric field satisfies is ▽ ² Φ=-ρ/ε ₀ , where ▽ ² represents the Laplacian operator, Φ represents the electric potential, ρ is the charge volume density, and ε ₀ is the vacuum permittivity; according to the condition that the positive charge and electron charge in the local simulation region ψ are close to zero in the steady state, the electric field simplifies To solve the Laplace equation ▽ ² Φ=0, use the existing finite difference method and combine the boundary conditions to obtain the potential of each point in the simulation area Ω; in addition to the triangular plane in each cuboid grid, there are also Preserve the electric field formed by the plasma for simulating the movement of charged ions; 4.2) According to the sedimentary particle type S, the particle number N _s and the sedimentary particle velocity distribution P _v , use the Monte Carlo MC method to randomly sample the sedimentary particles to obtain the sedimentary particles to be simulated Among them, i represents the serial number of deposited particles; Indicates deposited particles The position of , in the Cartesian coordinate system, Initial z=0, which means the upper boundary of the gas phase in the simulation area Ω, x, y are randomly generated according to the size of the simulation area Ω; _si indicates the type of deposited particles, which is set by the deposited particle components; represent The initial velocity of , given by the isotropic Maxwell velocity distribution produce, v _ix , v _iy , v _iz are the component velocities in the three directions of Cartesian coordinate system x, y, and z respectively, and -∞<v _ix , v _iy <+∞, v _iz >0, M _i is the mass of deposited particles , k is the Boltzmann constant, T _i is the deposition particle temperature; 4.3) According to the deposited particles type, use the space octree technology to quickly track the transport process of sediment particles, and calculate the parameters of sediment particles s _i , θ _i , If s _i ∈ ion, the deposited particle The motion is the motion caused by the acceleration caused by the electric field E formed by the plasma; the specific calculation is as follows: According to the electric field E in the simulation area Ω, Newton's law of motion and the current position of the deposited particles and speed Compute the new position after the Δt time interval and new speed then by and Determine the line segment AB, use the space octree technology to determine whether this line segment intersects with the surface of the depositional graph, if not intersect and this line segment does not leave the simulation area Ω, then set the new position and new speed as current location and speed Repeat this process until it intersects with the surface of the deposition pattern, and find the triangular plane _Δi that intersects with the deposition pattern; then calculate the deposition particles The angle θ _i between the incident direction and the normal of the triangular plane; and the surface reflection model is used to calculate the reflection or deposition, that is, if θ _i > θ ₀ , θ ₀ is the reflection threshold, and the deposition particles are calculated according to the reflection principle new location for new speed as current location and speed Go to step 4.3) to continue execution, otherwise go to step 4.4); If s _i ∈ neutral active particles, then the deposition particles The motion is constant speed Straight line motion with constant direction, use the space octree technology to quickly calculate the triangular plane _Δi where this line intersects with the deposition pattern Γ(t); then decide whether to scatter or adsorb according to the adsorption probability model, that is, generate a random number p, if p >S _c is the scattering, S _c is the adsorption probability, then randomly generate the scattering angle and rate and calculate the current position and speed Go to step 4.3) to continue execution; otherwise go to step 4.4); 4.4) According to the material at the deposition location and the deposition particle parameters s _i , θ _i , Select a specific deposition model, calculate and update the deposition amount of the triangular plane of the deposition location: If it is a plasma-enhanced vapor deposition process, the corresponding deposition model can be expressed as formula (2): in: Indicates the plasma-enhanced deposition coefficient, K _d (state) represents the deposition coefficient of neutral particles, and state represents the material; F _ion and F _d represent the ion flow and the neutral particle flow on the triangular plane ds on the deposition graph Γ(t) respectively ; θ represents the angle between the deposition particle and the normal direction of the triangular plane; If s _i ∈ ion, then If s _i ∈ neutral active particles, then in: Indicates the deposition amount on the triangular plane _Δi , express Plasma-enhanced deposition coefficient; K _d (state) represents neutral particles deposition coefficient; and respectively represent the current particle Representative ion flux and neutral particle flux; 4.5) Repeat steps 4.2)-4.4) to reach the specified number of particles and turn step 5). the 4. method as claimed in claim 1, is characterized in that, described step 5) specifically comprises: 5.1) According to the deposition amount of each triangular plane area and the unit normal vector to calculate the deposition velocity of each triangular plane 5.2) Calculate the stepping time interval ΔT according to the deposition speed of the triangular plane and the maximum advancing distance dist; Among them: dist represents the maximum advancing distance, ls represents the set of triangular planes that make up the sedimentary figure Γ(t); 5.3) According to the deposition speed of each triangular plane Using the numerical calculation method to solve the level set governing equation (5) to realize the motion of the surface; φ _t + V·▽φ＝0 (5) Among them: φ _t is the partial derivative of the level set function φ to time t, ▽φ is the gradient of φ at (x, y, z), V represents the velocity at each point (x, y, z) of the domain, approximately The deposition velocity of the triangular plane where it is located ; When using the numerical calculation method to solve, the domain of the level set function is discretized and the narrow-band level set method is used to calculate the velocity of the position near the surface of the deposition figure Γ (t); 5.4) Reinitialize the level set function according to formula (6) to make it satisfy formula (1): Among them: φ ₀ represents the current level set function, is the level set function to be calculated, and its steady-state solution is a signed distance function satisfying formula (1), sign(φ ₀ ) sign function, for the convenience of solving, it is smoothed into formula (7), ε>0 is very small positive number of 5.5) Re-triangulate the surface of the current deposition graph, and update the information in the cells according to the current deposition graph Γ(t), and replace the cells whose state has changed with new deposition materials. the