A method for accelerating 3D deposition simulation based on a hybrid algorithm of ray tracing and level set

Abstract

This invention relates to the field of semiconductor process simulation technology, and discloses a method and system for accelerating 3D deposition simulation based on a hybrid algorithm of ray tracing and level set. The method maps the substrate structure to a 3D mesh using a signed distance function; it traverses the mesh points to detect changes in the signed distance function value, extracting surface patches traversed by the zero isosurface; it emits probe rays from the gas source region and surface patches, and absorbs and scatters the rays according to the Lambertian scattering model and adhesion coefficient when they intersect the surface, statistically analyzing the gas flux at each patch; it maps the gas flux to a normal velocity field, solves the level set convection equation to update the signed distance function value, and iterates until the target time. This invention replaces the solution of 3D full-space partial differential equations, reducing computational complexity, while simultaneously simulating the shading effect and multiple reflections in the molecular flow region, making it suitable for deposition and etching simulations of high aspect ratio structures.

Description

Technical Field

[0001] This invention relates to the field of semiconductor process simulation technology, and discloses a method for accelerating 3D deposition simulation based on a hybrid algorithm of ray tracing and level set. Background Technology

[0002] In semiconductor process simulation, such as 3D topography simulation for LPCVD processes, current mainstream techniques typically employ a pure level set method combined with a continuum transport model to simulate the evolution and growth of thin films. In traditional methods, the entire gas space and solid structure are divided into a dense 3D uniform voxel mesh, with the physical surface implicitly represented by a signed distance function defined on this entire mesh. When solving for the gas flux reaching the surface, current techniques require establishing diffusion-reaction equations on the entire 3D gas phase mesh. This involves discretizing the massive 3D mesh containing the surface and the entire space using the finite difference method or finite volume method, transforming it into solving a giant sparse linear equation system AC=b. For 3D scenes with tens of millions of voxels, the time complexity of solving this linear equation system approaches O(N³), resulting in extremely long computation times and massive memory consumption. Actual test data from mainstream TCAD software shows that for a 3D trench structure with 10 million voxels, a single simulation using a pure level set method combined with a continuum transport model takes no less than 48 hours and consumes no less than 64GB of memory. Furthermore, LPCVD is a low-pressure process, resulting in a relatively long mean free path for gas molecules, which have entered the molecular flow region. Traditional global diffusion equations based on continuous media are applicable only when the mean free path of gas molecules is much smaller than the structural feature size (viscous flow region). This premise no longer holds under the low-pressure LPCVD process, making it difficult to accurately describe the shielding effect and multiple reflection mechanism of gas molecules within high aspect ratio trenches. The conformal error between simulation results and experimental SEM measurements is typically no less than 20%.

[0003] To address the shortcomings of the aforementioned continuous medium models, the industry has proposed a technical solution combining the level set method with Monte Carlo particle transport. This approach uses the Monte Carlo method to simulate the transport process of gas molecules to calculate surface flux, while combining it with the level set method to track the evolution of surface morphology. However, the Monte Carlo method relies on statistical sampling of a large number of random particles, resulting in inherent statistical noise in its flux calculations. Reducing this noise requires increasing the number of sampled particles, leading to a significant increase in computational load and slow convergence speed. Furthermore, while this approach can simulate the spatial motion of particles, it still lacks an accurate physical modeling method based on molecular kinematics for the shading effect and multiple reflections of gas molecules within the trenches in LPCVD molecular flow regions. Summary of the Invention

[0004] The main objective of this invention is to solve the technical problem of accurately simulating the shading effect and multiple reflections inside high aspect ratio structures in LPCVD molecular flow region deposition simulation, while avoiding the statistical noise and convergence speed problems of traditional Monte Carlo methods, and reducing computational complexity.

[0005] This invention provides a method for accelerating 3D deposition simulation based on a hybrid algorithm of ray tracing and level set, comprising: Construct a 3D simulation scene and configure the gas source area; map the initial base structure to a 3D mesh using a signed distance function, and set the level set evolution narrowband and narrowband width; Traverse the grid points of the 3D grid and check whether the sign distance function value of adjacent grid points changes sign; if so, determine that the zero isosurface passes through the grid edge between adjacent grid points, and use linear interpolation to calculate the coordinates of the intersection point of the zero isosurface and the grid edge; extract the surface patch formed by all intersection points, and normalize the gradient of the sign distance function as the normal vector of the surface patch. In a non-sequential physical simulation environment, probe rays are emitted from the gas source region and surface patches respectively. The rays propagate in a straight line in space, and the intersection detection between the rays and the surface is accelerated using a hierarchical bounding box tree. When the ray intersects the surface, the reflection direction is determined according to the scattering model, and the deposition absorption ratio and scattering ratio are determined according to the surface reaction coefficient. The corresponding weight of the deposition absorption ratio is added to the flux accumulator of the corresponding surface patch, and the flux accumulation result is normalized according to the area of ​​the surface patch. The ray weight of the scattered part is decayed and the tracing continues until the termination condition is met. The gas flux of each surface patch is calculated based on the final weight of each ray received by each surface patch and the incident angle. The gas flux is multiplied by the deposition rate coefficient and mapped to the normal velocity field of the level set evolution. The normal velocity field is extrapolated to the grid points in the narrow band along the gradient direction of the sign distance function. The level set convection equation is solved in the narrow band, the sign distance function value is updated, and the operation of traversing the grid points and updating the sign distance function value is repeated until the target deposition time is reached.

[0006] Furthermore, the gas source region is a surface parallel to the top of the substrate. The emission point of the rays is uniformly sampled on the gas source region, and the emission direction of the rays is cosine distributed along the normal hemisphere space of the gas source region. The scattering model adopts the Lambert scattering model, and the reflection direction of the Lambert scattering model is cosine distributed along the normal hemisphere space of the irradiated surface.

[0007] Furthermore, the sign change represents the sign distance function value changing from positive to negative or from negative to positive; Adjacent sign-changing grid points are x1 and x2, with corresponding sign distance function values ​​Φ1>0 and Φ2<0, respectively. The intersection point x0 of the zero isosurface and the grid edge is calculated using the following formula: In the formula, Indicates the sign change grid point The sign distance function value; Indicates the sign change grid point The sign distance function value; Represents the absolute value function; The surface patches are triangulated, and adjacent surface patches share vertices.

[0008] Furthermore, the surface reaction coefficient is the adhesion coefficient, and the value range of the adhesion coefficient is [0,1]; the deposition absorption ratio is the adhesion coefficient, and the scattering ratio is 1 minus the adhesion coefficient.

[0009] Furthermore, the ray weights are updated according to the following rules: In the formula, Indicates the weight of the scattered rays; Indicates the weight of the incident ray; This represents the adhesion coefficient.

[0010] Furthermore, the termination conditions include: The ray position exceeds the simulation boundary; the ray weight is lower than the preset threshold; the number of ray scatterings reaches the maximum scattering limit.

[0011] Furthermore, methods for obtaining gas flux include: The gas flux is obtained using the gas flux calculation formula, which is shown below: In the formula, Indicates the first Gas flux per surface patch; Indicates the first The area of ​​each surface patch; Indicates hitting the first The number of rays on each surface patch; Indicates hitting the first The first surface patch The weight carried by a ray at the moment of impact; Indicates the first The cosine of the angle of incidence of the ray; This represents the cosine function.

[0012] Furthermore, methods for obtaining the normal velocity field include: The normal velocity field is obtained using the formula for calculating the normal velocity field, which is shown below: In the formula, Indicates surface position Normal velocity field at the location; Represents the deposition rate constant; Indicates surface position Gas flux at the location; Indicates the molar mass of the deposited material; Indicates the density of the deposited material; This represents the adhesion coefficient.

[0013] Furthermore, the hierarchical bounding box tree adopts incremental updates for the dynamic surface updated after each evolution step; the level set convection equation is spatially discretized using the fifth-order WENO scheme and temporally discretized using the third-order TVD Runge-Kutta scheme. The normal velocity field is extrapolated to the grid points within the narrow band along the gradient direction of the sign distance function using a fast sweep method. Every preset number of steps, the symbolic distance function within the narrow band is reinitialized using a fast sweep method.

[0014] A 3D deposition simulation acceleration system based on a hybrid algorithm of ray tracing and level set includes: The scene construction and initialization module is used to build a 3D simulation scene and configure the gas source area; it maps the initial base structure to a 3D mesh through a signed distance function and sets the level set evolution narrowband and narrowband width. The surface patch extraction module is used to traverse the grid points of the three-dimensional grid, detect whether the sign distance function value of adjacent grid points has changed sign; if so, it determines that the zero isosurface passes through the grid edge between the adjacent grid points, calculates the coordinates of the intersection point of the zero isosurface and the grid edge using linear interpolation; extracts the surface patch formed by all intersection points, and normalizes the gradient of the sign distance function as the normal vector of the surface patch. The flux calculation module is used to emit probe rays from the gas source region and the surface patches respectively in a non-sequential physical simulation environment. The rays propagate in a straight line in space, and the intersection detection between the rays and the surface is accelerated by using a hierarchical bounding box tree. When the ray intersects the surface, the reflection direction is determined according to the scattering model, and the deposition absorption ratio and scattering ratio are determined according to the surface reaction coefficient. The corresponding weight of the deposition absorption ratio is added to the flux accumulator of the corresponding surface patch, and the flux accumulation result is normalized according to the area of ​​the surface patch. The ray weight of the scattered part is decayed and the tracing continues until the termination condition is met. The gas flux of each surface patch is calculated according to the final weight of each ray received by each surface patch and the incident angle. The morphology evolution driving module is used to multiply the gas flux by the deposition rate coefficient and map it into a normal velocity field of the level set evolution; the normal velocity field is extrapolated to the grid points in the narrow band along the gradient direction of the sign distance function; the level set convection equation is solved in the narrow band; the sign distance function value is updated; and the operation of traversing the grid points to update the sign distance function value is repeated until the target deposition time is reached.

[0015] In the technical solution provided by this invention, during the surface extraction stage, the zero isosurface is located by traversing three-dimensional grid points and detecting the sign change of the sign distance function value of adjacent points. Linear interpolation is then used to accurately calculate the intersection coordinates, enabling rapid location of the material's physical surface. This provides precise geometric boundaries for subsequent ray tracing, avoiding flux calculation errors caused by surface gaps. During the flux calculation stage, non-sequential ray tracing replaces the traditional three-dimensional full-space partial differential equation solution, reducing the calculation scope from the entire volume grid to two-dimensional surface patches. This fundamentally avoids the construction and solution of giant sparse linear equation systems, reducing computational complexity from O(N³) to the O(N²) order of magnitude, thus reducing computation time and memory usage. Simultaneously, through the Lambertian scattering model and the adhesion coefficient weight attenuation mechanism, the shading effect and multiple reflections of gas molecules in the LPCVD molecular flow region can be accurately simulated, overcoming the physical distortion problem of continuous medium models under low-pressure conditions. In the velocity field mapping stage, the local gas fluxes of each patch obtained by ray tracing are directly mapped to the normal velocity field of the level set evolution, establishing a quantitative correlation between physical flux and geometric evolution, ensuring the comparability of simulation results with actual processes. In the evolution update stage, the level set convection equations are solved only in a narrow band region near the zero isosurface, and the velocity field extrapolation and periodic reinitialization of the sign distance function are combined with the fast sweep method, further reducing the computational load while ensuring numerical accuracy; coupled with the incremental update of the hierarchical bounding box tree, the intersection detection efficiency of the dynamic surface after each evolution step is maintained, shortening the overall simulation time. Attached Figure Description

[0016] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings: Figure 1 This is a flowchart of a 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set in an embodiment of the present invention; Figure 2 This is a block diagram of a 3D deposition simulation acceleration system based on a hybrid algorithm of ray tracing and level set in an embodiment of the present invention. Detailed Implementation

[0017] Exemplary embodiments of the invention will now be described more fully with reference to the accompanying drawings. However, these exemplary embodiments can be implemented in many forms and should not be construed as limiting the invention to the embodiments set forth herein. Rather, these exemplary embodiments are provided to make the invention more comprehensive and complete, and to facilitate a full communication of the inventive concept to those skilled in the art. The same reference numerals in the drawings denote the same or similar elements, components, or parts, and therefore repeated descriptions of them will be omitted.

[0018] Subject to the technical concept of this invention, the features, structures, characteristics or other details described in a particular embodiment may be combined in one or more other embodiments in a suitable manner.

[0019] In the description of specific embodiments, the features, structures, characteristics, or other details described in this invention are intended to enable those skilled in the art to fully understand the embodiments. However, it is not excluded that those skilled in the art can practice the technical solutions of this invention without one or more of the specific features, structures, characteristics, or other details.

[0020] The flowcharts shown in the accompanying drawings are merely illustrative and do not necessarily include all content and operations / steps, nor do they necessarily have to be performed in the described order. For example, some operations / steps can be broken down, while others can be combined or partially combined; therefore, the actual execution order may change depending on the specific circumstances.

[0021] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically independent entities. That is, these functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices.

[0022] The terms “and / or” or “and / or” include all combinations of any one or more of the listed items.

[0023] See Figure 1 The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set described in this embodiment has the following specific steps: Step S1: Construct a three-dimensional simulation scene and configure the gas source area; map the initial base structure to the three-dimensional mesh through the symbolic distance function, and set the level set evolution narrowband and narrowband width.

[0024] This invention is used to simulate a multi-material hybrid deposition and etching process for high aspect ratio structures in semiconductor manufacturing. To adapt to this scenario, this application first constructs a three-dimensional simulation scenario and configures the gas source area.

[0025] In one embodiment of this invention, the gas source region adopts a Lambertian diffuse emission model, with the ray emission direction following a cosine distribution, consistent with the actual molecular emission law of LPCVD gas sources. The gas source region is a user-configured process input parameter, including its geometry, spatial location, normal direction, and total emission intensity. Specifically, the gas source region is a surface parallel to the top of the substrate. The user configures the shape, size, spatial coordinates, and total emission intensity of this surface according to the actual reaction chamber configuration. It should be noted that this configuration process can be set by the implementer according to the specific implementation scenario, and will not be elaborated here. The total emission intensity is measured in particles per second, used to convert ray statistics into absolute physical flux, ensuring that the simulated deposition rate can be directly correlated with experimental data.

[0026] Furthermore, the initial substrate structure is mapped onto a 3D mesh using a signed distance function. The signed distance function is a scalar function defined on a spatial mesh, where its absolute value is the distance from a spatial point to the nearest material surface. Positive values ​​represent the gas phase region (outside the material surface), negative values ​​represent the solid region (inside the material surface), and zero iso-faces correspond to the physical surface of the material. Through this mapping, the geometric information of the substrate structure is transformed into a numerical field on the mesh, providing a data foundation for subsequent surface extraction.

[0027] Secondly, the narrowband width for level set evolution is set. The narrowband refers to a strip-shaped region near the zero isosurface. Level set calculations and signed distance function updates are performed only at grid points within the narrowband; grid points outside the narrowband are not involved in the calculation. The narrowband method reduces the computational scope from the entire three-dimensional space to a thin layer near the surface, significantly reducing the computational cost of level set evolution. Since the selection of the narrowband width needs to balance computational efficiency and numerical stability—too small a width may cause the surface to move out of the narrowband during evolution, while too large a width increases unnecessary computation—3 to 5 grid layers exhibit optimal overall performance in typical applications. Therefore, in this embodiment of the invention, the narrowband width is set to 3 to 5 grid layers.

[0028] Finally, the simulation process parameters, surface material properties, target deposition time, and time step are defined. The simulation process parameters include the adhesion coefficient. The value range is [0,1], representing the probability of a deposition reaction occurring after gas molecules collide with the surface, determined by specific process conditions (temperature, pressure, type of reactant gas). Surface material properties include the molar mass and density of the deposited material. Simultaneously, a simulation boundary is defined; rays exceeding this boundary are considered escapes and tracking is terminated. The specific shape and size of the simulation boundary are determined by those skilled in the art based on the actual simulation scenario.

[0029] Step S2: Traverse the grid points of the 3D grid and check whether the sign distance function value of adjacent grid points changes sign; if so, determine that the zero isosurface passes through the grid edge between adjacent grid points, and use linear interpolation to calculate the coordinates of the intersection point of the zero isosurface and the grid edge; extract the surface patch formed by all intersection points, and normalize the gradient of the sign distance function as the normal vector of the surface patch.

[0030] Traverse all grid points of the 3D mesh and check whether the sign of the distance function value between adjacent grid points changes from positive to negative or from negative to positive. If a change of sign occurs, it indicates that a transition has occurred from the gas phase to the solid phase. At this point, determine the zero isosurface, that is, the physical surface of the material, which passes through the grid edge between the adjacent grid points.

[0031] Since the signed distance function changes approximately linearly at the scale of the mesh edges, the location of the zero isosurface can be precisely determined by the proportional relationship between the signed distance function values ​​on both sides. Compared to higher-order interpolation, linear interpolation minimizes computation while maintaining sufficient accuracy, making it suitable for frequent use in each evolution step.

[0032] Let adjacent sign-changing grid points be x1 and x2, with corresponding sign distance function values ​​Φ1>0 and Φ2<0, respectively. The intersection point x0 of the zero isosurface and the grid edge is calculated using the following formula: In the formula, Indicates the sign change grid point The sign distance function value; Indicates the sign change grid point The sign distance function value; This represents the absolute value function.

[0033] In the formula, the position of the zero isosurface is determined by the ratio of the absolute values ​​of the sign distance function of the grid points on both sides. That is, the side with the smaller absolute value of the sign distance function is closer to the material surface, and the intersection of the zero isosurface and the grid edge is closer to that side.

[0034] After extracting all intersection points, these intersection points are used to form surface patches. In this embodiment, the surface patches are triangulated, meaning that every three adjacent intersection points form a triangular patch. The normal vector of the patch is obtained by normalizing the gradient of the sign distance function corresponding to that location, which ensures that the normal vector direction is consistent with the physical surface and that it is topologically continuous. Adjacent patches share vertices, avoiding gaps in ray tracing or flux calculation errors caused by gaps between patches.

[0035] It should be noted that surface extraction only occurs within the narrow band; grid points outside the narrow band do not participate in sign change detection. If there are no sign-changing grid points within the narrow band, it means the surface has moved out of the current narrow band, at which point the narrow band adjustment mechanism is triggered.

[0036] Step S3: In a non-sequential physical simulation environment, probe rays are emitted from the gas source region and surface patches respectively; the rays propagate in a straight line in space, and the intersection detection between the rays and the surface is accelerated by using a hierarchical bounding box tree; when the ray intersects the surface, the reflection direction is determined according to the scattering model, and the deposition absorption ratio and scattering ratio are determined according to the surface reaction coefficient; the corresponding weight of the deposition absorption ratio is added to the flux accumulator of the corresponding surface patch, and the flux accumulation result is normalized according to the area of ​​the surface patch; the ray weight of the scattered part continues to be tracked after decay until the termination condition is met; the gas flux of each surface patch is calculated according to the final weight of each ray received by each surface patch and the incident angle.

[0037] In a non-sequential physical simulation environment, the surface patch extracted in step S2 is simultaneously treated as a receiver to receive incident rays and an emitter to emit scattered rays. Probe rays are emitted from the configured gas source region and the deposited surface patch, respectively.

[0038] When rays are emitted from the gas source region, the emission start points are uniformly sampled by area on the user-defined geometry of the gas source region. The emission direction follows a cosine distribution along the normal hemisphere of the gas source region, which is consistent with the actual emission law of gas molecules in the molecular flow region of LPCVD. The scattering model adopts the Lambertian scattering model. The initial weight of each ray is set to 1, representing a representative molecular beam. Each ray is independent of each other, and the entire ray tracing process can be executed concurrently on the GPU.

[0039] Rays propagate in straight lines in space. When a ray intersects a surface, the reflection direction is determined according to the Lambertian scattering model; that is, the reflection direction follows a cosine distribution along the normal hemisphere of the irradiated surface. According to the adhesion coefficient... Determine the deposition absorption ratio and the scattering ratio, where the deposition absorption ratio is: With the same adhesion coefficient, the scattering ratio is Specifically, the adhesion coefficient characterizes the probability that gas molecules will undergo chemisorption and participate in the deposition reaction after impacting a solid surface. At that time, all colliding molecules do not react and are completely scattered back to the gas phase; when At this time, all impacting molecules react and are completely absorbed by the surface. It should be noted that the adhesion coefficient is determined by specific process conditions, such as reaction temperature, pressure, type and partial pressure of reactant gases, and type of surface material, and is a known parameter set by those skilled in the art based on the actual process scenario.

[0040] When the ray has weight When impacting the surface sheet, according to the adhesion coefficient Probabilistic flow splitting is performed, where the depositional absorption portion has a weight of... The molecules undergo deposition reactions, and their contribution is included in the flux accumulator of the surface patch; scattering component: weighted at... The molecules did not react and left the surface by scattering, continuing to fly in the gas phase.

[0041] Ray weights are updated according to the following rules: In the formula, Indicates the weight of the scattered rays; Indicates the weight of the incident ray; This represents the adhesion coefficient.

[0042] It should be noted that in the molecular flow region of LPCVD, after gas molecules undergo adsorption-desorption processes on the solid surface, the desorption direction follows a cosine law, meaning that the probability density of desorbed molecules is highest in the normal direction and tends to zero in the tangential direction. This law has been confirmed by experimental observations in the field of molecular flow, namely Knudsen's cosine law, which is a standardized description of gas-surface interactions in the molecular flow region. Therefore, the embodiments of this invention employ a Lambert scattering model that matches the physical mechanism of the LPCVD molecular flow region.

[0043] The termination conditions for ray tracing include any of the following: the ray position exceeds the simulation boundary; the ray weight is lower than a preset threshold, which is set to 10 by default in this embodiment. -6 The number of ray scattering attempts reaches the maximum scattering limit, which is set to 10 by default in this embodiment. The selection of the weight threshold needs to balance computational accuracy and efficiency: a threshold that is too high will prematurely truncate the reflection chain and underestimate the contribution of multiple reflections to the flux; a threshold that is too low will increase the amount of ineffective computation. -6 In typical applications, a reasonable balance can be achieved between accuracy and efficiency. Similarly, the maximum number of scattering times can be set, as the ray weight has decayed to an extremely low level after multiple reflections, and its contribution to the flux is negligible.

[0044] For intersection detection, a hierarchical bounding box tree is used to accelerate the intersection detection of rays and surfaces. The hierarchical bounding box tree is a spatial partitioning-based acceleration structure that organizes surface patches into hierarchical bounding box levels, quickly eliminating impossible intersection patches and significantly reducing the computational cost of intersection detection. Since surface deformation is usually small between adjacent time steps, the position and topological relationships of most patches remain unchanged. Only the bounding box nodes corresponding to the changed patches need to be updated, and the update time accounts for no more than 5%, which can improve the overall simulation efficiency. Therefore, for the dynamic surface updated after each evolution step, incremental hierarchical bounding box tree updates are used instead of full reconstruction.

[0045] To reduce statistical noise in flux calculations, a low-discrepancy Halton sequence is used to generate ray emission directions, replacing purely random sampling. The Halton sequence is a low-discrepancy sequence that can achieve more uniform spatial coverage with the same number of samples, thereby reducing statistical fluctuations in flux calculations. Purely random sampling exhibits clustering and gaps with a limited number of samples, leading to differences in flux statistics across different runs; the Halton sequence, through deterministic low-discrepancy sampling, can obtain more stable statistical results with the same number of rays.

[0046] After tracing all rays, the gas flux of each surface patch is calculated based on the final weight of each ray received by each surface patch and the incident angle.

[0047] Preferably, in one embodiment of the present invention, the method for obtaining gas flux includes: The gas flux is obtained using the gas flux calculation formula, which is shown below: In the formula, Indicates the first Gas flux per surface patch; Indicates the first The area of ​​each surface patch; Indicates hitting the first The number of rays on each surface patch; Indicates hitting the first The first surface patch The weight carried by a ray at the moment of impact; Indicates the first The cosine of the angle of incidence of the ray; This represents the cosine function.

[0048] In the formula, the gas flux equals the sum of the weights carried by all incident rays across the surface area multiplied by the cosine of the incident angle. The larger the incident angle, the smaller the effective flux received across the surface area, which perfectly matches the actual physical law of gas molecules colliding with the surface per unit time and per unit area in the molecular flow region. This is then divided by the surface area. Normalization is performed to make the gas flux a physical quantity independent of the size of the surface, which facilitates the comparison of gas flux between surface patches of different sizes and the continuous mapping of the subsequent velocity field.

[0049] Step S4: Multiply the gas flux by the deposition rate coefficient to map it into the normal velocity field of the level set evolution; extrapolate the normal velocity field along the gradient direction of the sign distance function to the grid points in the narrow band, solve the level set convection equation in the narrow band, update the sign distance function value, and repeat the operation of traversing the grid points to update the sign distance function value until the target deposition time is reached.

[0050] First, the gas flux is multiplied by the deposition rate coefficient to map it to the normal velocity field of the level set evolution. The formula for calculating the normal velocity field is as follows: The formula for calculating the deposition rate constant is as follows: In the formula, Indicates surface position Normal velocity field at the location; Represents the deposition rate constant; Indicates surface position Gas flux at the location; Indicates the molar mass of the deposited material; Indicates the density of the deposited material; This represents the adhesion coefficient.

[0051] In the formula, the unit of gas flux is... multiplied by The velocity dimension was then obtained. This formula characterizes the deposition rate of the rays. It establishes a one-to-one correspondence between measurable values ​​in engineering and calculable values ​​in simulation; that is, the simulation input is the gas flux, and the output is the normal velocity field, used to characterize the surface position. The absolute deposition rate. Personnel can verify the simulation results using actual process experiments, such as SEM / TEM testing of the film thickness. If the simulated values ​​match the experimental values, the simulation model is accurate; if they do not match, the simulation model can be optimized by adjusting parameters such as the adhesion coefficient.

[0052] Since the convection equations for the level set require calculating velocity values ​​at all grid points within the narrow band, and the above steps only define the normal velocity field on the surface patch, the surface velocity must be extended along the surface normal into the narrow band. In this embodiment of the invention, a fast sweep method is used to extrapolate along the gradient direction of the signed distance function, ensuring the smoothness and continuity of the velocity field within the narrow band, and providing complete velocity field data for the subsequent solution of the level set equations. The fast sweep method achieves rapid extrapolation of the velocity field by solving the equations in alternating directions.

[0053] Solving the convection equations for level sets within a narrow band: In the equation, The symbolic distance function represents the symbolic distance to the surface; Indicates the simulation time; Represents the normal velocity field; The gradient of the signed distance function is represented. This represents the magnitude of the gradient of the sign distance function.

[0054] In the equation, The unknown in the equation is denoted by , which states that if the surface travels at a velocity ... During the advancement, the symbolic distance function Response time.

[0055] It should be noted that in this embodiment of the invention, a fifth-order WENO scheme is used for spatial discretization, and a third-order TVD Runge-Kutta scheme is used for temporal discretization to ensure numerical accuracy and stability. The WENO (Weighted Essentially Non-Oscillatory) scheme can effectively suppress numerical oscillations while maintaining high-order accuracy, making it suitable for capturing sharp morphologies that may appear on the deposition surface; the TVD (Total Variation Diminishing) Runge-Kutta scheme ensures that the total variation does not increase during time progression, avoiding non-physical oscillations in the numerical solution.

[0056] Meanwhile, regarding time progression, since the convection equation of the level set is a hyperbolic partial differential equation, the numerical solution is stable only when the time step is less than a certain critical value. Therefore, each time step is dynamically constrained by the CFL conditions (Courant–Friedrichs–Lewy conditions). In the formula, The basic unit representing the time discretization of the level set convection equations is the time step. This represents the dimensionless safety factor in the CFL condition, with a value range of [value missing]. ; The basic unit representing the spatial discretization of the level set convection equation, namely the grid size, limits the maximum safe displacement of the surface within a single time step; This represents the maximum normal velocity at all grid points within the narrow band at the current moment.

[0057] In the formula, the time step This represents the time interval during which the surface topography progresses from the current state to the next state. CFL is the dimensionless safety factor in the CFL condition, ranging from (0,1), with a default value of 0.5 in this embodiment. This parameter ensures that the distance the surface advances within a time step does not exceed half the mesh size, preventing numerical distortion caused by surface propagation across meshes. A value that is too large will cause the calculation to reach an unstable critical state; a value that is too small, while ensuring stability, increases the total number of steps and computation time. Mesh size. It is the side length of each voxel in the 3D mesh, representing the spatial resolution. The finer the mesh, the higher the resolution. The smaller the value, the more the maximum safe displacement of the surface within a single time step under the CFL condition. The finer the mesh, the smaller the allowed time step, because at a finer mesh scale, the physical space allowed for the surface to advance per step is smaller. This represents the velocity value of the region with the fastest propagation on the current surface, determined by the gas flux calculation results mentioned above. Regions with higher gas flux exhibit higher deposition rates. The larger the step size, the greater the numerical stability of the entire surface. To ensure the numerical stability of the entire surface, the fastest propagation speed must be used as the constraint benchmark. If the step size is determined with a smaller speed, the fastest region will violate the CFL condition, leading to local divergence.

[0058] It should be noted that the 3D mesh is a uniform mesh, and the mesh size is... It is set by the user according to the simulation accuracy requirements and computing resources. The smaller the mesh size, the higher the spatial resolution and the more accurate the surface positioning, but the increased mesh size leads to increased memory usage and computational load. The larger the value, the lower the value. Those skilled in the art can determine the trade-off between accuracy and efficiency. The specific value to be taken.

[0059] After 5-10 time steps of evolution, the symbol distance function within the narrowband is reinitialized using a fast sweep method, i.e., the SDF is corrected to satisfy the given conditions. During the evolution of the level set, due to numerical dissipation, the symbolic distance function gradually deviates from the ideal symbolic distance property, leading to a decrease in the accuracy of surface normal vector calculation, which in turn affects the intersection detection accuracy of ray tracing. Reinitialization corrects the SDF to an accurate symbolic distance function by solving the reinitialization equation, while keeping the position of the zero isosurface unchanged, thus ensuring the long-term accuracy of surface extraction.

[0060] If the surface approaches the narrowband boundary within a single time step, it may shift out of the current narrowband in the next time step, resulting in a loss of computational range. Early detection and expansion of the narrowband ensures that the zero isosurface remains within the central region of the narrowband, avoiding computational interruptions or accuracy degradation due to insufficient narrowband range. Therefore, before the evolution of each time step, the displacement of the zero isosurface along the normal direction is calculated, i.e. If the displacement of the zero isosurface along the normal exceeds 50% of the distance from the current narrowband inner boundary to the zero isosurface, then narrowband adjustment is triggered, and the narrowband is expanded bidirectionally by 3 to 5 meshes along the surface normal in the next iteration.

[0061] After completing the evolution of each time step, return to step S2, the surface patch extraction process, and start the cycle for the next time step until the cumulative evolution time reaches the target total deposition time, and output the updated 3D deposition morphology.

[0062] In summary, a 3D simulation scene is constructed, and the gas source area is configured. The initial base structure is mapped to a 3D mesh using a signed distance function, and the level set evolution narrowband and narrowband width are set. The mesh points of the 3D mesh are traversed, and the sign of the signed distance function value of adjacent mesh points is checked. If it does, it is determined that the zero isosurface passes through the mesh edge between adjacent mesh points, and the coordinates of the intersection point between the zero isosurface and the mesh edge are calculated using linear interpolation. The surface patch formed by all intersection points is extracted, and the gradient of the signed distance function is normalized as the normal vector of the surface patch. In a non-sequential physical simulation environment, probe rays are emitted from the gas source area and the surface patch respectively. The rays propagate in a straight line in space, and the intersection detection between the rays and the surface is accelerated using a hierarchical bounding box tree. When the ray intersects the surface, the reflection direction is determined according to the scattering model. The deposition absorption and scattering ratios are determined according to the surface reaction coefficient. The corresponding weights of the deposition absorption ratios are added to the flux accumulator of the corresponding surface patches, and the flux accumulation results are normalized according to the area of ​​the surface patches. The weights of the scattered rays are attenuated and the tracing continues until the termination condition is met. The gas flux of each surface patch is calculated based on the final weights of each ray received by each surface patch and the incident angle. The gas flux is multiplied by the deposition rate coefficient and mapped to the normal velocity field of the level set evolution. The normal velocity field is extrapolated to the grid points in the narrow band along the gradient direction of the sign distance function. The level set convection equation is solved in the narrow band, the sign distance function value is updated, and the operation of traversing the grid points and updating the sign distance function value is repeated until the target deposition time is reached.

[0063] Please see Figure 2 The second objective of this invention is to provide a 3D deposition simulation acceleration system based on a hybrid algorithm of ray tracing and level set, which is achieved through the above steps, specifically including: The scene construction and initialization module 101 is used to construct a three-dimensional simulation scene, configure the gas source area, map the initial base structure to the three-dimensional mesh through the signed distance function, and set the level set evolution narrowband and narrowband width. The surface patch extraction module 102 is used to traverse the grid points of the three-dimensional grid, detect whether the sign distance function value of adjacent grid points changes sign; if so, it determines that the zero isosurface passes through the grid edge between the adjacent grid points, calculates the coordinates of the intersection point of the zero isosurface and the grid edge using linear interpolation; extracts the surface patch formed by all intersection points, and normalizes the gradient of the sign distance function as the normal vector of the surface patch; The flux calculation module 103 is used to emit probe rays from the gas source area and the surface patches respectively in a non-sequential physical simulation environment. The rays propagate in a straight line in space, and the intersection detection between the rays and the surface is accelerated by using a hierarchical bounding box tree. When the rays intersect the surface, the reflection direction is determined according to the scattering model, and the deposition absorption ratio and scattering ratio are determined according to the surface reaction coefficient. The corresponding weight of the deposition absorption ratio is added to the flux accumulator of the corresponding surface patch, and the flux accumulation result is normalized according to the area of ​​the surface patch. The ray weight of the scattered part continues to be tracked after decay until the termination condition is met. The gas flux of each surface patch is calculated according to the final weight of each ray received by each surface patch and the incident angle. The morphology evolution driving module 104 is used to multiply the gas flux by the deposition rate coefficient and map it into a normal velocity field of the horizontal set evolution; to extrapolate the normal velocity field along the gradient direction of the sign distance function to the grid points in the narrow band; to solve the horizontal set convection equation in the narrow band; to update the sign distance function value; and to repeat the operation of traversing the grid points to update the sign distance function value until the target deposition time is reached.

[0064] Program code for performing the operations of this invention can be written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Java and C++, and conventional procedural programming languages ​​such as C or similar languages. The program code can execute entirely on the user's computing device, partially on the user's device, as a standalone software package, partially on the user's computing device and partially on a remote computing device, or entirely on a remote computing device or server. In cases involving remote computing devices, the remote computing device can be connected to the user's computing device via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computing device (e.g., via the Internet using an Internet service provider).

[0065] In summary, this invention can be implemented in hardware, or as software modules running on one or more processors, or a combination thereof. Those skilled in the art will understand that in practice, general-purpose data processing devices such as microprocessors or digital signal processors (DSPs) can be used to implement some or all of the functions of some or all of the components according to the embodiments of the invention. This invention can also be implemented as a device or apparatus program (e.g., a computer program and computer program product) for performing part or all of the methods described herein. Such programs implementing the invention can be stored on a computer-readable medium or can take the form of one or more signals. Such signals can be downloaded from an Internet website, provided on a carrier signal, or provided in any other form.

[0066] The above specific embodiments further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the present invention is not inherently related to any specific computer, virtual device, or electronic device, and various general-purpose devices can also implement the present invention. The above are merely specific embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0067] The various embodiments in this specification are described in a progressive manner. The same or similar parts between the various embodiments can be referred to each other. Each embodiment focuses on describing the differences from other embodiments.

[0068] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A method for accelerating 3D deposition simulation based on a hybrid algorithm of ray tracing and level set, characterized in that, include: Construct a 3D simulation scene and configure the gas source area; The initial basis structure is mapped to a 3D mesh using a signed distance function, and the level set evolution narrowband and narrowband width are set. Traverse the grid points of the three-dimensional grid and check whether the sign distance function value of adjacent grid points has changed sign; If so, determine that the zero isosurface passes through the grid edge between the adjacent grid points, and use linear interpolation to calculate the coordinates of the intersection point between the zero isosurface and the grid edge; extract the surface patch formed by all intersection points, and normalize the gradient of the signed distance function as the normal vector of the surface patch; In a non-sequential physical simulation environment, probe rays are emitted from the gas source region and the surface patches respectively; the rays propagate in a straight line in space, and the intersection detection between the rays and the surface is accelerated using a hierarchical bounding box tree; when the rays intersect the surface, the reflection direction is determined according to the scattering model, and the deposition absorption ratio and scattering ratio are determined according to the surface reaction coefficient; the corresponding weight of the deposition absorption ratio is added to the flux accumulator of the corresponding surface patch, and the flux accumulation result is normalized according to the area of ​​the surface patch; the ray weight of the scattered part is attenuated and the tracing continues until the termination condition is met; the gas flux of each surface patch is calculated according to the final weight of each ray received by each surface patch and the incident angle; The gas flux is multiplied by the deposition rate coefficient and mapped to the normal velocity field of the level set evolution. The normal velocity field is extrapolated to the grid points in the narrow band along the gradient direction of the sign distance function. The level set convection equation is solved in the narrow band, the sign distance function value is updated, and the operation of traversing the grid points to update the sign distance function value is repeated until the target deposition time is reached.

2. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 1, characterized in that, The gas source region is a surface parallel to the top of the substrate. The starting point of the ray emission is uniformly sampled on the gas source region, and the ray emission direction is cosine distributed along the normal hemispherical space of the gas source region. The scattering model adopts the Lambert scattering model, and the reflection direction of the Lambert scattering model is cosine distributed along the normal hemispherical space of the irradiated surface.

3. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 1, characterized in that, The sign change refers to the sign distance function value changing from positive to negative or from negative to positive. Adjacent sign-changing grid points are x1 and x2, with corresponding sign distance function values ​​Φ1>0 and Φ2<0, respectively. The intersection point x0 of the zero isosurface and the grid edge is calculated using the following formula: In the formula, Indicates the sign change grid point The sign distance function value; Indicates the sign change grid point The sign distance function value; Represents the absolute value function; The surface patches are triangulated, and adjacent surface patches share vertices.

4. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 1, characterized in that, The surface reaction coefficient is the adhesion coefficient, and the value range of the adhesion coefficient is [0,1]; the deposition absorption ratio is the adhesion coefficient, and the scattering ratio is 1 minus the adhesion coefficient.

5. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 1, characterized in that, The ray weights are updated according to the following rules: In the formula, Indicates the weight of the scattered rays; Indicates the weight of the incident ray; This represents the adhesion coefficient.

6. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 1, characterized in that, The termination conditions include: The ray position exceeds the simulation boundary; the ray weight is lower than the preset threshold; the number of ray scatterings reaches the maximum scattering limit.

7. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 1, characterized in that, The method for obtaining the gas flux includes: The gas flux is obtained according to the gas flux calculation formula, which is as follows: In the formula, Indicates the first Gas flux per surface patch; Indicates the first The area of ​​each surface patch; Indicates hitting the first The number of rays on each surface patch; Indicates hitting the first The first surface patch The weight carried by a ray at the moment of impact; Indicates the first The angle of incidence of the ray; Indicates the first The cosine of the angle of incidence of the ray; This represents the cosine function.

8. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 4, characterized in that, The method for obtaining the normal velocity field includes: The normal velocity field is obtained according to the formula for calculating the normal velocity field, which is shown below: In the formula, Indicates surface position Normal velocity field at the location; Represents the deposition rate constant; Indicates surface position Gas flux at the location; Indicates the molar mass of the deposited material; Indicates the density of the deposited material; This represents the adhesion coefficient.

9. The 3D deposition simulation acceleration method based on a hybrid algorithm of ray tracing and level set as described in claim 1, characterized in that, The hierarchical bounding box tree uses incremental updates for the dynamic surface updated after each evolution step; the level set convection equations are spatially discretized using the fifth-order WENO scheme and temporally discretized using the third-order TVD Runge-Kutta scheme. The normal velocity field is extrapolated to the grid points within the narrow band along the gradient direction of the sign distance function using a fast sweep method. Every preset number of steps, the symbolic distance function within the narrow band is reinitialized using a fast sweep method.

10. A 3D deposition simulation acceleration system based on a hybrid algorithm of ray tracing and level set, characterized in that, include: The scene construction and initialization module is used to build a 3D simulation scene and configure the gas source area; The initial basis structure is mapped to a 3D mesh using a signed distance function, and the level set evolution narrowband and narrowband width are set. The surface patch extraction module is used to traverse the grid points of the three-dimensional grid and detect whether the sign distance function value of adjacent grid points has changed sign; If so, determine that the zero isosurface passes through the grid edge between the adjacent grid points, and use linear interpolation to calculate the coordinates of the intersection point between the zero isosurface and the grid edge; extract the surface patch formed by all intersection points, and normalize the gradient of the signed distance function as the normal vector of the surface patch; The flux calculation module is used to emit probe rays from the gas source region and the surface patch respectively in a non-sequential physical simulation environment. The rays propagate in a straight line in space, and the intersection detection between the rays and the surface is accelerated by using a hierarchical bounding box tree. When the ray intersects the surface, the reflection direction is determined according to the scattering model, and the deposition absorption ratio and scattering ratio are determined according to the surface reaction coefficient. The corresponding weight of the deposition absorption ratio is added to the flux accumulator of the corresponding surface patch, and the flux accumulation result is normalized according to the area of ​​the surface patch. The ray weight of the scattered part is decayed and the tracing continues until the termination condition is met. The gas flux of each surface patch is calculated based on the final weight of each ray received by each surface patch and the incident angle. The morphology evolution driving module is used to multiply the gas flux by the deposition rate coefficient and map it into a normal velocity field of the level set evolution; the normal velocity field is extrapolated to the grid points in the narrow band along the gradient direction of the sign distance function; the level set convection equation is solved in the narrow band; the sign distance function value is updated; and the operation of traversing the grid points to update the sign distance function value is repeated until the target deposition time is reached.

Images