Thin film coating process simulation and optimization system based on digital twinning

Abstract

The application belongs to the technical field, and discloses a thin film coating process simulation and optimization system based on digital twinning; the system is composed of a data acquisition module, a feature recognition and analysis module, a behavior digital twinning simulation module, a coating uniformity evaluation module, a process parameter optimization module and a coating distribution prediction module. By obtaining a three-dimensional roughness digital model of a substrate, a micron-level local rough area is accurately identified and its geometric feature parameters are extracted; a coating growth behavior digital twinning model is constructed based on the geometric feature parameters and process parameters, and the coating growth process of the local rough area is accurately simulated; a global coating uniformity quantitative evaluation system is established, and a machine learning and genetic algorithm are applied to optimize the process parameter combination. The system breaks through the limitations of traditional static shadow assumption, realizes real-time tracking and accurate simulation of micro-morphology evolution in the deposition process, and significantly improves the coating growth prediction accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field, more particularly, the present application relates to a thin film coating process simulation and optimization system based on digital twinning. BACKGROUND

[0002] As a key process in modern manufacturing, thin film coating technology plays an irreplaceable role in microelectronics, optics, aerospace, biomedicine and other fields. With the increasing demand for product performance in high-end manufacturing, coating uniformity, consistency and reliability have become the core indicators of product quality. However, in actual production, the influence of the micro roughness of the substrate surface on the quality of the thin film coating is increasingly prominent.

[0003] However, the existing thin film coating process simulation technology is limited by static thinking, which cannot accurately reflect the dynamic evolution law in the deposition process. In actual production, the high aspect ratio trench structure of microelectronic devices, the complex curved surface of optical elements and the precise surface of aerospace coatings all show obvious time-varying characteristics, while the traditional model is based on the initial topography for static shadow calculation, ignoring the changes in local geometry caused by gradual coating filling. This simplified method of "one-time calculation and permanent use" leads to significant deviations between predicted results and actual coating distribution in practical applications, especially in deep trench bottoms, sharp corner transition zones and micro-depressions. In industry practice, integrated circuit interconnection layers, barrier layer deposition frequently appear in the bottom weak area, OLED display cathode coating appears discontinuity in the edge area, and precise optical mirror reflector has uneven reflectivity in the curvature change area. These problems can be traced back to the prediction deviation caused by the model's inability to track topography evolution. In addition, existing optimization algorithms rely on empirical parameters and simplified models, lacking the ability to accurately describe the dynamic process at the microscale, resulting in strong dependence on optimization results and poor portability. Process engineers have to repeat a large number of experiments on each new product, greatly delaying the product development cycle and increasing research and development costs.

[0004] In view of this, the present application proposes a thin film coating process simulation and optimization system based on digital twinning to solve the above problems. SUMMARY

[0005] In order to overcome the above-mentioned defects of the prior art, in order to achieve the above-mentioned purpose, the present application provides the following technical scheme: a thin film coating process simulation and optimization system based on digital twinning, comprising:

[0006] A data acquisition module is configured to obtain a three-dimensional roughness digital model of a substrate surface and a process parameter digital model of a processing equipment during a thin film coating process.

[0007] a feature recognition and analysis module configured to identify a local rough area with micron-level unevenness on the substrate surface based on the three-dimensional roughness digital model, and extract geometric feature parameters of the local rough area;

[0008] a behavior digital twin simulation module configured to construct a coating growth behavior digital twin model of the local rough area according to the geometric feature parameters of the local rough area and the process parameter digital model, and obtain a local coating thickness distribution prediction value of the local rough area in the thin film coating processing process;

[0009] a coating uniformity evaluation module configured to calculate a global coating uniformity index of the substrate surface according to the local coating thickness distribution prediction value and a coating thickness distribution prediction value of a flat area on the substrate surface;

[0010] a process parameter optimization module configured to optimize the process parameter digital model based on the global coating uniformity index and a preset uniformity threshold value, and obtain an optimized process parameter combination;

[0011] a coating distribution prediction module configured to update the coating growth behavior digital twin model based on the optimized process parameter combination, and obtain a global coating thickness distribution prediction value of the optimized substrate surface;

[0012] The modules are connected through wired and / or wireless means to realize data transmission between the modules.

[0013] Preferably, the method for extracting the geometric feature parameters of the local rough area comprises:

[0014] performing grid division on the three-dimensional roughness digital model to obtain height distribution data of each grid cell on the substrate surface;

[0015] calculate an angle between a local surface normal vector of each grid cell and a global normal vector of the substrate surface according to the height distribution data, denoted as a local normal deviation angle;

[0016] identify, according to the local normal deviation angle and a preset deviation angle threshold value, a grid cell with a local normal deviation angle greater than the preset deviation angle threshold value on the substrate surface, to constitute the local rough area;

[0017] extract the geometric feature parameters of the local rough area, the geometric feature parameters including an average depth, a maximum width of the local rough area, and a distribution variance of the local normal deviation angle.

[0018] Preferably, the method for constructing the coating growth behavior digital twin model comprises:

[0019] constructing a three-dimensional geometric digital twin sub-model of the local rough area according to the geometric feature parameters of the local rough area;

[0020] According to the process parameter digital model, the incident angle distribution and deposition rate distribution of the deposition particles in the thin film coating process are obtained.

[0021] Based on the three-dimensional geometric digital twin sub-model and the incident angle distribution of the deposition particles, the shadow effect factor of each grid cell in the local rough area is calculated, and the shadow effect factor is the effective area ratio of the grid cell receiving the deposition particles.

[0022] According to the shadow effect factor and the deposition rate distribution, the local coating thickness growth rate of each grid cell in the local rough area is calculated.

[0023] Based on the local coating thickness growth rate and the processing time, the local coating thickness distribution prediction value of the local rough area is obtained.

[0024] Preferably, the method for calculating the global coating uniformity index comprises:

[0025] The coating thickness distribution prediction value of the flat area of the substrate surface is obtained, denoted as the flat area coating thickness mean value;

[0026] According to the local coating thickness distribution prediction value, the coating thickness mean value of the local rough area is calculated.

[0027] The ratio of the coating thickness mean value of the local rough area to the coating thickness mean value of the flat area is calculated, denoted as the local uniformity deviation.

[0028] According to the local uniformity deviation and the area ratio of the local rough area on the substrate surface, the global coating uniformity index is calculated.

[0029] Preferably, the method for optimizing the process parameter digital model to obtain the optimized process parameter combination comprises:

[0030] A mapping relationship model between the process parameter digital model and the global coating uniformity index is constructed;

[0031] Based on the mapping relationship model, a genetic algorithm is used to iteratively optimize the process parameter digital model to obtain a candidate process parameter combination that satisfies the preset uniformity threshold of the global coating uniformity index.

[0032] The candidate process parameter combination is verified for feasibility, and the unfeasible process parameter combination is removed to obtain the optimized process parameter combination.

[0033] Preferably, the method for calculating the shadowing factor comprises:

[0034] According to the three-dimensional geometric digital twin sub-model of the local rough area, a local surface normal vector of each grid cell in the local rough area is obtained;

[0035] According to the incident angle distribution of the deposited particles, the included angle between the incident direction of the deposited particles and the local surface normal vector is calculated, denoted as the incident deviation angle;

[0036] According to the incident deviation angle and the height difference of adjacent grid cells in the local rough area, the shadowing blocking ratio of the grid cell is calculated;

[0037] According to the shadowing blocking ratio and the incident deviation angle, the shadowing factor of the grid cell is calculated.

[0038] Preferably, the method for performing grid division on the three-dimensional roughness digital model to obtain the height distribution data of each grid cell on the substrate surface comprises:

[0039] Scanning the substrate surface to obtain three-dimensional topography data of the substrate surface;

[0040] Based on the three-dimensional topography data, a feature saliency distribution map of the substrate surface is constructed;

[0041] According to the feature saliency distribution map, a multi-scale feature region of the substrate surface is identified, the multi-scale feature region comprising a macro curvature dominant region, a micron-scale roughness dominant region, and a nanoscale texture dominant region; based on the distribution of the multi-scale feature region, a dynamic grid division scheme is generated, the dynamic grid division scheme comprising: assigning a first grid density to the macro curvature dominant region, the grid cell size of the first grid density being determined based on the average curvature radius of the macro curvature dominant region; assigning a second grid density to the micron-scale roughness dominant region, the grid cell size of the second grid density being determined based on the characteristic length of the micron-scale roughness dominant region, the characteristic length being the spatial scale corresponding to the maximum value of the height variation rate in the micron-scale roughness dominant region; assigning a third grid density to the nanoscale texture dominant region, the grid cell size of the third grid density being determined based on the texture period of the nanoscale texture dominant region, the texture period being obtained by Fourier transform on the height distribution data of the nanoscale texture dominant region;

[0042] According to the dynamic meshing scheme, the three-dimensional roughness digital model is meshed to obtain height distribution data of each mesh element of the substrate surface, wherein the height distribution data of the mesh element of the macroscopic curvature dominant region is extracted from the three-dimensional topography data through local quadratic surface fitting, the height distribution data of the mesh element of the micrometer scale roughness dominant region is extracted from the three-dimensional topography data through local spline interpolation, and the height distribution data of the mesh element of the nanometer scale texture dominant region is extracted from the three-dimensional topography data through local high-frequency filtering.

[0043] The height distribution data of the mesh element is corrected for consistency, and the consistency correction includes adjusting the height distribution data of the mesh element at the junction of the macroscopic curvature dominant region, the micrometer scale roughness dominant region and the nanometer scale texture dominant region based on the boundary transition characteristics between the multi-scale feature regions by using a boundary smoothing algorithm.

[0044] Preferably, the constructing a three-dimensional geometric digital twin sub-model of the local rough region according to the geometric feature parameters of the local rough region comprises:

[0045] Based on the geometric feature parameters of the local rough region, the dynamic change characteristics of the local rough region are extracted, and the dynamic change characteristics include the morphology evolution trend of the local rough region in the thin film coating processing process;

[0046] According to the dynamic change characteristics, an initial three-dimensional geometric digital twin sub-model of the local rough region is constructed, and the initial three-dimensional geometric digital twin sub-model is generated by parameterized modeling of the geometric feature parameters of the local rough region.

[0047] Based on the morphology evolution trend, the initial three-dimensional geometric digital twin sub-model is updated to obtain a three-dimensional geometric digital twin sub-model of the local rough region, and the updating includes adjusting the control point position of the spline curve according to the change of the processing time or the process condition.

[0048] Preferably, the method for obtaining the optimized global coating thickness distribution prediction value of the substrate surface comprises:

[0049] Based on the optimized process parameter combination, the incident angle distribution and the deposition rate distribution of the deposition particles are updated;

[0050] According to the updated incident angle distribution and deposition rate distribution of the deposition particles, the local coating thickness distribution prediction value of the local rough region and the coating thickness distribution prediction value of the flat region are recalculated;

[0051] The updated local coating thickness distribution prediction value is integrated with the coating thickness distribution prediction value of the flat region to obtain the optimized global coating thickness distribution prediction value of the substrate surface.

[0052] Preferably, the calculation of the local coating thickness growth rate of each grid cell in the local rough region comprises:

[0053] Based on the process parameter digital model, the spatial heterogeneity characteristics of the deposition rate distribution in the thin film coating process are obtained, and the spatial heterogeneity characteristics include the deposition rate gradient distribution and spatial correlation parameters of different regions of the substrate surface.

[0054] According to the three-dimensional geometric digital twin sub-model, the spatial position coordinates of each grid cell in the local rough region on the substrate surface are determined.

[0055] Based on the spatial position coordinates and the spatial heterogeneity characteristics of the deposition rate distribution, the local deposition rate correction value corresponding to each grid cell in the local rough region is calculated.

[0056] According to the local deposition rate correction value and the shadow effect factor, the local coating thickness growth rate of each grid cell in the local rough region is calculated, and the local coating thickness growth rate is the product of the local deposition rate correction value and the shadow effect factor.

[0057] The technical effects and advantages of the thin film coating process simulation and optimization system based on digital twinning of the present application are as follows:

[0058] The present application breaks through the limitations of traditional static shadow assumption by constructing a dynamically updated digital twin model, realizes real-time tracking and accurate simulation of micro-topography evolution in the deposition process. The dynamic simulation capability of the present application improves the accuracy of coating growth prediction, especially on high aspect ratio structures and complex curved substrates. For substrates containing deep groove structures, the prediction accuracy is also improved, providing a more reliable basis for thin film performance design. By accurately capturing the dynamic evolution law of the shadow effect with the change of topography, the present application not only optimizes the coating uniformity, but also greatly reduces the generation of local weak areas, so that the coating structure defect rate is reduced to a low level. At the same time, the process parameter optimization efficiency is improved, and the number of tests and material consumption are reduced. BRIEF DESCRIPTION OF DRAWINGS

[0059] Figure 1 The schematic diagram of the thin film coating process simulation and optimization system based on digital twinning of the present application. DETAILED DESCRIPTION

[0060] With reference to the accompanying drawings: the technical solutions in the embodiments of the present application will be described clearly and completely, obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative labor belong to the scope of protection of the present application.

[0061] The present application provides a thin film coating process simulation and optimization system based on digital twinning. The execution body of the system includes: data acquisition equipment, coating simulation platform, process optimization server, network transmission equipment, etc. which can be regarded as general computing nodes of the present application, and the coating simulation platform includes but is not limited to: at least one of coating process monitoring system, process parameter optimization platform, digital twinning simulation system.

[0062] The present application provides a thin film coating process simulation and optimization system based on digital twinning, comprising the following steps:

[0063] A data acquisition module is used to obtain a three-dimensional roughness digital model of a substrate surface and a process parameter digital model of a processing equipment in a thin film coating process.

[0064] A feature recognition and analysis module is used to identify a local rough area with micron-level unevenness on the substrate surface based on the three-dimensional roughness digital model, and extract geometric feature parameters of the local rough area.

[0065] A behavior digital twinning simulation module is used to construct a coating growth behavior digital twinning model of the local rough area according to the geometric feature parameters of the local rough area and the process parameter digital model, and obtain a local coating thickness distribution prediction value of the local rough area in the thin film coating process.

[0066] A coating uniformity evaluation module is used to calculate a global coating uniformity index of the substrate surface according to the local coating thickness distribution prediction value and a coating thickness distribution prediction value of a flat area of the substrate surface.

[0067] A process parameter optimization module is used to optimize the process parameter digital model based on the global coating uniformity index and a preset uniformity threshold, and obtain an optimized process parameter combination.

[0068] A coating distribution prediction module is used to update the coating growth behavior digital twinning model based on the optimized process parameter combination, and obtain a global coating thickness distribution prediction value of the optimized substrate surface.

[0069] Each module is connected through wired and / or wireless means to realize data transmission between modules.

[0070] The present application obtains the microstructure characteristics of the substrate surface through a three-dimensional roughness digital model; identifies local rough areas and accurately extracts geometric feature parameters; constructs a digital twin model of coating growth behavior, accurately simulates the coating growth process of local rough areas; calculates the global coating uniformity index, deeply evaluates the coating quality; optimizes the process parameter combination based on the uniformity index; and updates the digital twin model using the optimized parameters, realizing high-precision and high-reliability thin film coating process optimization.

[0071] In the embodiments of the present application, referring to Figure 1 , the schematic diagram of the thin film coating processing process simulation and optimization system based on digital twin of the present application, in the present example, the thin film coating processing process simulation and optimization system based on digital twin, comprising:

[0072] The data acquisition module is used to obtain the three-dimensional roughness digital model of the substrate surface and the process parameter digital model of the processing equipment in the thin film coating processing process;

[0073] In the present embodiment, first, advanced surface topography measurement equipment is deployed, including atomic force microscope (AFM), laser confocal microscope, white light interferometer and other high-precision surface measurement systems, to accurately scan the substrate surface.

[0074] The scanning range of the substrate surface is usually 10 μm x 10 μm to 100 mm x 100 mm, and the resolution can reach nanometer level. The surface height data obtained by scanning is subjected to noise filtering processing, and methods such as Gaussian filtering or median filtering are used to eliminate measurement noise to obtain high-quality surface topography original data. The filtered surface topography data is subjected to three-dimensional reconstruction to generate a three-dimensional roughness digital model of the substrate surface, which represents the micro-topography characteristics of the substrate surface in the form of three-dimensional point cloud or grid. At the same time, the process parameters of the thin film coating processing equipment are collected, including deposition rate (such as 0.1-10 nm / s), substrate temperature (such as room temperature to 800℃), deposition angle (such as 0° to 90°), background pressure (such as 10 -6 to 10 -3 Pa), reaction gas flow (such as 0-100 sccm) and other key process parameters.

[0075] These process parameters are digitized to construct a process parameter digital model, which contains the numerical range of each process parameter, the mutual relationship and its influence mechanism on coating growth. The three-dimensional roughness digital model and the process parameter digital model are time-labeled and associated matched to ensure their consistency in time and space, providing basic data support for subsequent analysis.

[0076] The feature recognition and analysis module is used to identify local rough areas with micron-level unevenness in the substrate surface based on the three-dimensional roughness digital model, and extract the geometric feature parameters of the local rough areas;

[0077] In this embodiment, the three-dimensional roughness digital model is meshed, and the mesh cell size is dynamically adjusted according to the complexity of the surface features, usually 1 / 10 to 1 / 20 of the surface feature size. The local surface normal vector at each mesh cell is calculated, which is calculated by the spatial position of the adjacent points and reflects the local tilt state of the surface. The global normal vector of the substrate surface is determined, which is usually taken as the normal direction of the main plane of the substrate surface, representing the direction of the ideal smooth surface. The angle between the local surface normal vector and the global normal vector of each mesh cell is calculated to obtain the local normal deviation angle. The larger the local normal deviation angle, the more serious the surface tilt at that location, and the higher the roughness. A preset angle threshold is set, usually 5° to 15°, and the specific value is determined according to the coating process and substrate type. Identify the mesh cells with a local normal deviation angle greater than the preset angle threshold, and cluster these cells to form a continuous region, which is the local rough region. For the identified local rough region, extract the geometric feature parameters, including: the average depth of the local rough region (the average distance from the points in the rough region to the fitted plane), the maximum width (the maximum extension size of the rough region in the surface plane), and the distribution variance of the local normal deviation angle (reflecting the consistency of the surface tilt in the rough region). These geometric feature parameters comprehensively describe the three-dimensional morphological characteristics of the local rough region, providing key inputs for subsequent coating growth behavior simulation.

[0078] The behavior digital twin simulation module is used to construct a coating growth behavior digital twin model of the local rough region based on the geometric feature parameters of the local rough region and the process parameter digital model, and to obtain a local coating thickness distribution prediction value of the local rough region in the thin film coating processing process.

[0079] In this embodiment, based on the geometric feature parameters of the local rough region, a three-dimensional geometric digital twin sub-model of the region is constructed. This sub-model accurately replicates the micro-topographic features of the local rough region, including surface undulations, slope changes, and boundary shapes, etc.

[0080] Key parameters of the deposition process are extracted from the process parameter digital model, including the incident angle distribution and deposition rate distribution of the deposited particles. The incident angle distribution describes the angle range and its probability distribution of the deposited particles arriving at the substrate surface, which is usually related to the geometric structure of the deposition equipment and the process conditions. The deposition rate distribution describes the number of deposited particles arriving at each location on the substrate surface per unit time, reflecting the spatial non-uniformity of the deposition process.

[0081] Based on the three-dimensional geometric digital twin sub-model and the incident angle distribution of the deposited particles, the shadow effect factor of each grid cell in the local rough area is calculated. The shadow effect refers to the phenomenon that part of the area is blocked and cannot receive the deposited particles due to the unevenness of the surface. The shadow effect factor is defined as the proportion of the effective area of the grid cell receiving the deposited particles, and the value ranges from 0 to 1. The smaller the value, the stronger the shadow effect. Combined with the shadow effect factor and the deposition rate distribution, the local coating thickness growth rate of each grid cell in the local rough area is calculated. The local coating thickness growth rate is influenced by both the shadow effect and the local deposition rate, and reflects the non-uniformity of the coating growth on the rough surface. According to the local coating thickness growth rate and the preset processing time, the local coating thickness distribution prediction value of the local rough area is calculated cumulatively. This prediction value presents the growth of the coating in the local rough area in a three-dimensional distribution, and intuitively reflects the influence of surface roughness on the uniformity of the coating.

[0082] The coating uniformity evaluation module is used to calculate the global coating uniformity index of the substrate surface according to the local coating thickness distribution prediction value and the coating thickness distribution prediction value of the flat area of the substrate surface.

[0083] In this embodiment, the same coating growth behavior digital twin model is applied to calculate the coating thickness distribution prediction value of the flat area of the substrate surface. The flat area refers to the area with a local normal deviation angle less than a preset deviation angle threshold. The coating growth in these areas is relatively uniform and less affected by the shadow effect. The coating thickness distribution of the flat area is statistically analyzed to calculate its average value, denoted as the flat area coating thickness average. The local coating thickness distribution prediction value of the local rough area is statistically analyzed to calculate its average value, denoted as the local rough area coating thickness average. The local uniformity deviation is calculated by comparing the local rough area coating thickness average with the flat area coating thickness average, and calculating the ratio of the two. The closer the local uniformity deviation value is to 1, the smaller the difference in coating thickness between the local rough area and the flat area, and the better the coating uniformity. The farther the deviation is from 1, the greater the difference and the worse the uniformity. Considering the distribution of the local rough area on the entire substrate surface, the area ratio of the local rough area is calculated. The larger the area ratio, the more significant the influence of the local rough area on the global coating uniformity.

[0084] Considering the local uniformity deviation and the area ratio of the local rough area, the global coating uniformity index is calculated, which comprehensively reflects the overall uniformity of the coating thickness distribution on the substrate surface. The specific calculation formula is:

[0085] Global coating uniformity index ; wherein QU is the global coating uniformity index, JU is the local uniformity deviation, and CU is the local rough area ratio. The index value range is 0-1, the closer to 1 indicates that the coating is more uniform, and the closer to 0 indicates that the coating is more uneven.

[0086] The process parameter optimization module is configured to optimize the process parameter digital model based on the global coating uniformity index and the preset uniformity threshold, and obtain an optimized process parameter combination.

[0087] In this embodiment, the preset uniformity threshold is set, which is usually 0.85-0.95, and the specific value is determined according to the application requirement of the thin film coating. The global coating uniformity index is compared with the preset uniformity threshold, if the global coating uniformity index is lower than the preset uniformity threshold, the process parameters need to be optimized to improve the coating uniformity. Based on the coating growth behavior digital twin model, a mapping relationship model between the process parameters and the global coating uniformity index is constructed. The mapping relationship model is trained by machine learning methods (such as neural network, support vector machine, etc.), which can predict the global coating uniformity index under different process parameter combinations. Genetic algorithm is applied to search for the global optimization of the process parameters. In the genetic algorithm, each individual represents a set of process parameter combinations, the fitness function is the global coating uniformity index, and the goal is to maximize the global coating uniformity index or make it exceed the preset uniformity threshold. The initial population is randomly generated with several sets of process parameter combinations, and the population size is usually 50-200 individuals. Through selection, crossover and mutation operations, the average fitness of the population is improved generation by generation. The iteration number is usually 50-500 generations, or until the optimal fitness no longer improves significantly. From the final population, the process parameter combination that meets the preset uniformity threshold is selected as the candidate process parameter combination. The candidate process parameter combination is verified for feasibility, and the realizability of each parameter combination in the actual equipment is evaluated. The feasibility verification includes parameter range check, parameter combination rationality analysis, equipment limitation condition test, etc. The unfeasible process parameter combination is removed, and the group with the highest global coating uniformity index is selected from the remaining feasible solutions as the optimized process parameter combination.

[0088] The coating distribution prediction module is configured to update the coating growth behavior digital twin model based on the optimized process parameter combination, and obtain the predicted value of the global coating thickness distribution on the substrate surface after optimization.

[0089] In this embodiment, the optimized process parameter combination is applied to the coating growth behavior digital twin model to update the relevant parameters in the model. Based on the optimized process parameters, the incident angle distribution and deposition rate distribution of the deposited particles are recalculated. The optimized process parameters usually change the kinetic characteristics of the deposition process, thereby affecting the incident angle distribution and deposition rate distribution. Using the updated incident angle distribution and deposition rate distribution, the shadowing effect factor of the local rough area and the local coating thickness growth rate are recalculated. The shadowing effect factor changes with the change of the incident angle distribution, thereby affecting the growth pattern of the local coating thickness. Based on the updated local coating thickness growth rate, the local coating thickness distribution prediction value of the local rough area is recalculated. At the same time, using the same updated parameters, the coating thickness distribution prediction value of the flat area is recalculated. The updated local rough area coating thickness distribution prediction value and the flat area coating thickness distribution prediction value are integrated to form the optimized global coating thickness distribution prediction value of the substrate surface. The global coating thickness distribution prediction value visually displays the coating uniformity under the optimized process parameters in a three-dimensional visualization form, providing a scientific basis for the final determination of the process parameters. By comparing the global coating thickness distribution prediction values before and after optimization, the optimization effect is quantitatively evaluated, and the effectiveness of the optimization scheme is verified. Finally, the optimized process parameter combination and the expected coating thickness distribution result are output as a process guide for actual production.

[0090] In this embodiment, the method for extracting the geometric feature parameters of the local rough area includes:

[0091] Grid division is performed on the three-dimensional roughness digital model to obtain height distribution data of each grid element of the substrate surface;

[0092] According to the height distribution data, the included angle between the local surface normal vector of each grid element and the global normal vector of the substrate surface is calculated, denoted as the local normal deviation angle;

[0093] According to the local normal deviation angle and the preset angle threshold, the grid elements with a local normal deviation angle greater than the preset angle threshold in the substrate surface are identified to constitute the local rough area;

[0094] The geometric feature parameters of the local rough area are extracted, including the average depth, the maximum width of the local rough area, and the distribution variance of the local normal deviation angle.

[0095] In this embodiment, tetrahedral, hexahedral or hybrid meshing algorithms are applied to the three-dimensional roughness digital model to generate a mesh structure suitable for surface topography analysis. The mesh cell size is adaptively adjusted according to the complexity of the surface features, with higher mesh density in complex regions and lower mesh density in flat regions. The typical mesh cell size ranges from 50 nm to 5 μm, ensuring accurate capture of micron-scale surface features. For each mesh cell, the node coordinates and height values are obtained to form the height distribution data H(x, y) of the cell. The height data is represented as the distance perpendicular to the substrate's main plane, with positive values indicating protrusions and negative values indicating depressions. A least squares method is used to fit a local surface to the height distribution data of each mesh cell and its neighborhood, typically using a second or third order polynomial fit.

[0096] Based on the fitted surface, the local surface normal vector n_l at the center point of the mesh cell is calculated. The local surface normal vector is perpendicular to the tangent plane of the fitted surface, reflecting the inclination direction of the surface at that point. The global normal vector n_g of the substrate surface is determined, typically taking the normal direction of the substrate's main plane (e.g., [0, 0, 1]). The angle θ between the local surface normal vector and the global normal vector, i.e., the local normal deviation angle, is calculated: The local normal deviation angle ranges from 0° to 90°, with a larger value indicating a more severe surface tilt. A preset deviation angle threshold θ_th is set, which is used to distinguish between flat and rough regions. The threshold value should consider material properties, coating process, and quality requirements, typically within the range of 5° to 15°. The substrate surface is scanned to identify mesh cells with a local normal deviation angle θ greater than the preset deviation angle threshold θ_th. A connected component analysis algorithm is used to cluster adjacent mesh cells that meet the condition to form continuous local rough regions. For each identified local rough region, its geometric feature parameters are extracted: the average depth d_avg, calculated as the average distance of all points in the rough region to the fitted plane; the maximum width w_max, measuring the maximum extension size of the rough region in the surface plane; the distribution variance of the local normal deviation angle, calculated as the statistical variance of the local normal deviation angles of all mesh cells in the rough region, reflecting the consistency of surface tilt in the rough region. These geometric feature parameters constitute the feature vector of the local rough region, providing quantitative basis for subsequent coating growth behavior simulation.

[0097] In this embodiment, the method for constructing a digital twin model of coating growth behavior includes:

[0098] According to the geometric feature parameters of the local rough region, a three-dimensional geometric digital twin sub-model of the local rough region is constructed;

[0099] According to the process parameter digital model, the incident angle distribution and deposition rate distribution of the deposited particles in the thin film coating process are obtained;

[0100] Based on the three-dimensional geometric digital twin sub-model and the incident angle distribution of the deposition particles, a shadow effect factor of each grid cell in the local rough area is calculated, and the shadow effect factor is the effective area proportion of the grid cell receiving the deposition particles;

[0101] According to the shadow effect factor and the deposition rate distribution, the local coating thickness growth rate of each grid cell in the local rough area is calculated;

[0102] Based on the local coating thickness growth rate and the processing time, the local coating thickness distribution prediction value of the local rough area is obtained.

[0103] In this embodiment, based on the geometric feature parameters (average depth, maximum width, local normal deviation angle distribution variance) of the local rough area, a parameterized geometric model of the area is constructed. The parameterized model is represented by B-spline surface or NURBS surface, which can accurately reconstruct the three-dimensional topography of the local rough area. The parameterized geometric model is finely meshed, and the grid density is adaptively adjusted according to the local curvature to ensure sufficient grid accuracy in areas with large curvature changes.

[0104] The fine grid model is compared with the original three-dimensional roughness digital model to verify the geometric reconstruction accuracy, and the root mean square error is usually required to be less than 10% of the original data resolution. After verification, a three-dimensional geometric digital twin sub-model of the local rough area is generated, which completely retains the micro-topographic features of the local rough area. Based on the process parameter digital model, the geometric structure and process conditions of the thin film coating processing equipment are analyzed, and the incident characteristics of the deposition particles are extracted. A physical model is established to calculate the incident angle distribution P(θ,φ) of the deposition particles, where θ is the polar angle (the angle with the surface normal), and φ is the azimuth angle. The incident angle distribution is usually related to parameters such as device type, target-substrate distance, working pressure, etc.

[0105] Using the Monte Carlo method or analytical model, the deposition rate distribution R(x,y) of the deposition particles at each position on the substrate surface is calculated. The deposition rate distribution takes into account factors such as device geometry, target size, particle scattering, etc., and the unit is usually nm / s. Based on the three-dimensional geometric digital twin sub-model and the incident angle distribution of the deposition particles, ray tracing analysis is performed on each grid cell in the local rough area. Virtual light rays are emitted from different incident angles to simulate the incident path of the deposition particles, and it is detected whether the light rays are blocked by the surrounding surface.

[0106] The proportion of incident light rays that each grid cell can receive is counted, and the shadow effect factor SF(x,y) is calculated:

[0107] ; wherein P(θ_i, φ_i) is the incident angle distribution probability, representing the probability of a deposited particle incident from direction (θ_i, φ_i), SB_i is the shadow blocking ratio, representing the degree of blocking of particles from direction (θ_i, φ_i), and α_i is the incident angle, i.e. the angle between the incident direction and the local surface normal vector. The shadowing effect factor ranges from 0 to 1, 0 representing complete blocking and 1 representing no blocking, and intermediate values representing partial blocking states.

[0108] The local coating thickness growth rate of each grid cell in the local rough area is calculated by combining the shadowing effect factor SF(x, y) and the deposition rate distribution R(x, y) . The local coating thickness growth rate directly reflects the influence of surface roughness on coating growth and embodies the spatial non-uniformity of coating growth. According to the local coating thickness growth rate and the preset processing time t, the time integral of the coating thickness is calculated: . If the change in surface topography during coating growth is not significant, it can be simplified as T(x, y) = G(x, y) x t. Finally, the local coating thickness distribution prediction value T(x, y) of the local rough area is obtained, which presents the growth of the coating in the local rough area in the form of a three-dimensional distribution.

[0109] In this embodiment, the method for calculating the global coating uniformity index comprises:

[0110] The coating thickness distribution prediction value of the flat area of the substrate surface is obtained, denoted as the flat area coating thickness mean value;

[0111] According to the local coating thickness distribution prediction value, the coating thickness mean value of the local rough area is calculated;

[0112] The ratio of the coating thickness mean value of the local rough area to the coating thickness mean value of the flat area is calculated, denoted as the local uniformity deviation;

[0113] According to the local uniformity deviation and the area ratio of the local rough area on the substrate surface, the global coating uniformity index is calculated.

[0114] In this embodiment, the coating growth behavior digital twin model is applied to the identified flat area of the substrate surface (the area with a local normal angle less than a preset angle threshold) to calculate the coating thickness distribution T_f(x, y) of the flat area. Since the surface inclination of the flat area is small, the shadowing effect is weak, and the coating growth is relatively uniform, but it is still affected by the deposition rate distribution. Statistical analysis is performed on the flat area coating thickness distribution T_f(x, y) to calculate its mean value μ_f, standard deviation σ_f, and uniformity index (the ratio of the standard deviation to the mean value) The mean value of the coating thickness of the flat area is taken as a reference value, representing the coating thickness on an ideal flat surface. The mean value of the local coating thickness distribution of the local rough area is calculated, and the standard deviation and the uniformity index of the local rough area are calculated The mean value of the coating thickness of the local rough area reflects the overall influence of the surface roughness on the coating thickness. The ratio of the mean value of the coating thickness of the local rough area to the mean value of the coating thickness of the flat area is calculated, and the local uniformity deviation is obtained The local uniformity deviation UD reflects the difference in coating thickness between the rough area and the flat area. UD=1 indicates that the thicknesses of the two areas are the same, UD<1 indicates that the thickness of the rough area is smaller than that of the flat area (common in strong shadow effect recesses), and UD>1 indicates that the thickness of the rough area is greater than that of the flat area (common in local protrusions or edge areas). The distribution of the local rough area on the entire substrate surface is statistically analyzed, and the total area of the local rough area A_r and the total area of the substrate surface A_t are calculated, and the area ratio is obtained The area ratio AR reflects the macroscopic distribution characteristics of the surface roughness of the substrate, and the larger the AR, the higher the overall roughness of the surface. Considering the local uniformity deviation UD and the area ratio AR, the global coating uniformity index GU is calculated: GU=1-|1-UD|×AR. The global coating uniformity index GU has a value range of 0-1, and the closer the value to 1, the more uniform the coating, and the closer the value to 0, the more serious the unevenness. When the local uniformity deviation UD is close to 1 (the thicknesses of the rough area and the flat area are close) or the area ratio AR is close to 0 (the rough area is very small), the global uniformity index GU is close to 1, indicating that the global coating uniformity is good.

[0115] In this embodiment, the method for optimizing the process parameter digital model to obtain the optimized process parameter combination comprises the following steps:

[0116] A mapping relationship model between the process parameter digital model and the global coating uniformity index is constructed.

[0117] Based on the mapping relationship model, a genetic algorithm is used to iteratively optimize the process parameter digital model to obtain a candidate process parameter combination that satisfies the preset uniformity threshold of the global coating uniformity index.

[0118] The candidate process parameter combination is subjected to feasibility verification, and the unfeasible process parameter combination is removed to obtain the optimized process parameter combination.

[0119] In this embodiment, based on the coating growth behavior digital twin model, a large number of sample data of process parameter combinations and corresponding global coating uniformity indicators are collected. The sampling range covers all dimensions of the process parameter space, ensuring the representativeness and diversity of the samples. Common process parameters include deposition rate, substrate temperature, deposition angle, background pressure, reaction gas flow, etc. A mapping relationship model f(P)→GU is constructed between the process parameters and the global coating uniformity indicator using machine learning methods, where P is the process parameter vector and GU is the global coating uniformity indicator. The machine learning methods that can be selected include artificial neural network (ANN), support vector machine (SVM), random forest (RF), or Gaussian process regression (GPR), etc. The trained mapping relationship model is cross-validated to ensure that the model accuracy and generalization ability meet the requirements, usually requiring a prediction root mean square error of less than 5%. Genetic algorithm is applied to search for the global optimization of process parameters. In the genetic algorithm, each individual represents a set of process parameter combinations P_i, and the fitness function is the global coating uniformity indicator GU_i=f(P_i). The optimization goal is to maximize GU or make GU exceed the preset uniformity threshold GU_th. An initial population containing N process parameter combinations is randomly generated, and the population size N is usually 50-200. The process parameters of each individual are randomly generated within the allowed range. The fitness of each individual in the initial population, i.e. the global coating uniformity indicator, is calculated.

[0120] The roulette or tournament selection method is used to select individuals for reproduction according to the fitness ratio. Individuals with higher fitness have a higher probability of being selected. The selected individuals are subjected to crossover operation to generate new offspring individuals. The crossover methods include single-point crossover, multi-point crossover, or uniform crossover, etc. Some individuals are subjected to mutation operation to randomly change the values of some process parameters, in order to increase the population diversity and search ability. Elite individuals (a small number of individuals with the highest fitness) are directly entered into the next generation to ensure that excellent genes are not lost. The selection, crossover, mutation and elite preservation steps are repeated for multiple generations (usually 50-500 generations) until the termination condition is met. The termination condition can be to reach the maximum number of iterations, or the optimal fitness does not improve significantly for consecutive generations. From the final population, process parameter combinations with global coating uniformity indicators GU exceeding the preset uniformity threshold GU_th are selected to form a set of candidate process parameter combinations.

[0121] The candidate process parameter combinations are verified for feasibility to evaluate their implementability on the actual equipment. The feasibility verification includes: parameter range checking to ensure that each parameter is within the range allowed by the equipment; parameter combination rationality analysis to check whether there are conflicts or incompatibility between parameters; equipment limitation condition checking such as power limitation, temperature upper limit, gas flow range, etc.; process stability evaluation to analyze whether the parameter combination will cause process instability or difficulty in control. Unfeasible process parameter combinations are eliminated, and one or several groups with the highest global coating uniformity index are selected from the remaining feasible solutions as the optimized process parameter combinations. If there are multiple groups of parameter combinations with similar uniformity indexes, other factors (such as energy consumption, processing efficiency, cost, etc.) can be further considered for final selection.

[0122] In this embodiment, the calculation method of the shadow effect factor includes:

[0123] According to the three-dimensional geometric digital twin sub-model of the local rough area, the local surface normal vector of each grid element in the local rough area is obtained;

[0124] According to the incident angle distribution of the deposited particles, the included angle between the incident direction of the deposited particles and the local surface normal vector is calculated, which is denoted as the incident deflection angle;

[0125] According to the incident deflection angle and the height difference of adjacent grid elements in the local rough area, the shadow shielding ratio of the grid element is calculated;

[0126] According to the shadow shielding ratio and the incident deflection angle, the shadow effect factor of the grid element is calculated.

[0127] In this embodiment, the geometric information of each grid element is extracted from the three-dimensional geometric digital twin sub-model of the local rough area, including node coordinates, element connection relationship and surface topology. The local surface normal vector n_lo(x,y) at the center point of each grid element is obtained by applying local surface fitting (such as least squares method). The local surface normal vector is perpendicular to the tangent plane of the fitted surface, and its direction follows the right-hand rule, pointing to the outside of the substrate. According to the deposition source characteristics in the process parameter digital model, the incident angle distribution P(θ,φ) of the deposited particles is determined. The incident angle distribution describes the probability distribution of the deposited particles arriving at the substrate surface from various directions, where θ is the polar angle (the included angle with the global normal vector), and φ is the azimuth angle. The incident angle space is discretely sampled to generate a series of incident direction vectors. The number of sampling points is usually 50-200, and the sampling strategy can be uniform sampling or sampling according to the importance of the incident angle distribution. For the i-th incident direction vector d_i and the local surface normal vector n_lo(x,y) of each grid element, the included angle between the two is calculated to obtain the incident deflection angle of the i-th incident direction The incident angle a_i ranges from 0° to 180°, and when a_i > 90°, it represents that the incident direction enters from the back, at which time the direction cannot be deposited. For the direction with incident angle a_i ≤ 90°, the shadow blocking effect is continued to be calculated. A ray is emitted from each grid cell in the opposite direction of the incident direction vector d_i, and it is detected whether the ray intersects with the surrounding surface. If the ray intersects with the surrounding surface, it represents that the deposited particles in the direction are blocked. For each incident direction d_i, the shadow blocking ratio SB_i of the grid cell in the direction is calculated. If the ray does not intersect with any surface, SB_i = 0 (no blocking); if the ray intersects with the surface, SB_i = 1 (complete blocking).

[0128] Combined with the incident angle a_i and the shadow blocking ratio SB_i, the effective deposition factor EF_i of the grid cell in the incident direction is calculated The effective deposition factor considers the influence of the incident angle on the deposition efficiency (in accordance with the cosine law) and the shadow blocking effect.

[0129] According to the incident angle distribution P(θ_i, φ_i) and the effective deposition factor EF_i, the comprehensive shadow effect factor SF of the grid cell is calculated The comprehensive shadow effect factor SF is normalized to the range of 0-1, representing the average effective receiving area ratio after considering the incident angle distribution. SF = 1 represents no shadow effect, SF = 0 represents complete shadow blocking, and the intermediate value represents a partially blocked state. The shadow effect factor SF directly affects the local coating growth rate and is a key parameter of the digital twin model of the coating growth behavior.

[0130] In this embodiment, the three-dimensional roughness digital model is meshed to obtain the height distribution data of each grid cell of the substrate surface, including:

[0131] The substrate surface is scanned by an atomic force microscope or a laser scanning microscope to obtain three-dimensional topographic data of the substrate surface;

[0132] Based on the three-dimensional topographic data, a feature saliency distribution map of the substrate surface is constructed, and the feature saliency distribution map is obtained by calculating the local curvature and height change rate of each measurement point in the three-dimensional topographic data, wherein the local curvature is the trace of the surface curvature tensor at the measurement point, and the height change rate is the square root of the square sum of the height difference between the measurement point and its neighborhood measurement point;

[0133] According to the feature saliency distribution map, a multi-scale feature region of the substrate surface is identified, the multi-scale feature region including a macro curvature dominant region, a micron-scale roughness dominant region and a nanoscale texture dominant region, and the identification method includes threshold segmentation on the feature saliency distribution map, wherein the macro curvature dominant region corresponds to a low saliency threshold range, the micron-scale roughness dominant region corresponds to a medium saliency threshold range, and the nanoscale texture dominant region corresponds to a high saliency threshold range;

[0134] Based on the distribution of the multi-scale feature region, a dynamic grid division scheme is generated, including: assigning a first grid density to the macro curvature dominant region, the grid cell size of the first grid density being determined based on the average curvature radius of the macro curvature dominant region; assigning a second grid density to the micron-scale roughness dominant region, the grid cell size of the second grid density being determined based on the characteristic length of the micron-scale roughness dominant region, the characteristic length being the spatial scale corresponding to the maximum value of the height variation rate in the micron-scale roughness dominant region; and assigning a third grid density to the nanoscale texture dominant region, the grid cell size of the third grid density being determined based on the texture period of the nanoscale texture dominant region, the texture period being obtained by Fourier transform on the height distribution data of the nanoscale texture dominant region;

[0135] According to the dynamic grid division scheme, the three-dimensional roughness digital model is divided into grids to obtain the height distribution data of each grid cell of the substrate surface, wherein the grid cell height distribution data of the macro curvature dominant region is extracted from the three-dimensional topography data by local quadratic surface fitting, the grid cell height distribution data of the micron-scale roughness dominant region is extracted from the three-dimensional topography data by local spline interpolation, and the grid cell height distribution data of the nanoscale texture dominant region is extracted from the three-dimensional topography data by local high-frequency filtering.

[0136] The grid cell height distribution data is corrected for consistency, which includes adjusting the grid cell height distribution data at the junction of the macro curvature dominant region, the micron-scale roughness dominant region and the nanoscale texture dominant region based on the boundary transition characteristics between the multi-scale feature regions by using a boundary smoothing algorithm to eliminate the influence of the discontinuity of the grid density between the regions on the height distribution data.

[0137] In this embodiment, a suitable surface measurement device is selected to perform a precise scan of the substrate surface. For microscale features, a laser scanning microscope (LSM) or a white light interferometer (WLI) is used, with a typical scan range of 100 pm x 100 pm to 10 mm x 10 mm and a vertical resolution of 10-100 nm. For nanoscale features, an atomic force microscope (AFM) is used, with a typical scan range of 1 pm x 1 pm to 50 pm x 50 pm and a vertical resolution of 0.1-1 nm. For multiscale features, a multi-device joint measurement or a multi-resolution scanning strategy can be employed to ensure complete capture of features at different scales. The raw scan data is pre-processed, including noise filtering, defect repair, and reference correction, to obtain high-quality three-dimensional topography data Z(x, y) (height values). Feature analysis is performed on the three-dimensional topography data Z(x, y) to calculate the local curvature k(x, y) and the height variation rate g(x, y) for each measurement point (x, y).

[0138] The local curvature k(x, y) is calculated as the trace of the surface curvature tensor: where k1 and k2 are the principal curvatures. The curvature tensor is obtained by fitting a quadratic surface to the local surface and then solving its differential geometry properties. The height variation rate g(x, y) is calculated as the square root of the sum of the squares of the height differences between the measurement point and its N neighboring points:

[0139] where Z(x_i, y_i) represents the height value of the i-th point in the neighborhood, i.e., the height value of the i-th measurement point adjacent to (x, y), and Z(x, y) is the height value of the current point under consideration.

[0140] Based on the local curvature k(x, y) and the height variation rate g(x, y), a feature saliency distribution map S(x, y) is constructed for the substrate surface: where w1 and w2 are weight coefficients used to balance the contributions of curvature and height variation rate. The feature saliency distribution map S(x, y) reflects the geometric complexity of each region of the surface, with higher values indicating more prominent geometric features.

[0141] A multi-threshold segmentation is performed on the feature saliency distribution map S(x, y) to identify multiscale feature regions. A low saliency threshold S_low, a medium saliency threshold S_mid, and a high saliency threshold S_high are set to divide the substrate surface into three types of regions: a macro curvature dominant region (S(x, y) < S_low) with slow geometric variation, mainly exhibiting large-scale curvature changes; a microscale roughness dominant region (S_low ≤ S(x, y) < S_high) with obvious microscale concave-convex features; and a nanoscale texture dominant region (S(x, y) ≥ S_high) with high-frequency nanoscale texture features.

[0142] Based on the distribution of multi-scale feature regions, a dynamic mesh partitioning scheme is designed. For macro curvature dominant regions, a first mesh density is assigned, and the mesh cell size is where r_curv is the average radius of curvature of the region, and k1 is a safety factor (usually 5-10). For micrometer-scale roughness dominant regions, a second mesh density is assigned, and the mesh cell size is where L_feat is the feature length of the region, which is defined as the spatial scale corresponding to the maximum height variation rate, and k2 is a safety factor (usually 5-10). For nanometer-scale texture dominant regions, a third mesh density is assigned, and the mesh cell size is where T_texture is the texture period of the region, which is obtained by performing a two-dimensional fast Fourier transform (2D-FFT) on the height distribution data of the region and analyzing the power spectral density, and k3 is a safety factor (usually 5-10).

[0143] According to the dynamic mesh partitioning scheme, the three-dimensional roughness digital model is meshed. For macro curvature dominant regions, local quadratic surface fitting is used to extract the height distribution data of the mesh cells from the three-dimensional topography data. For micrometer-scale roughness dominant regions, local spline interpolation is used to extract the height distribution data of the mesh cells from the three-dimensional topography data. Spline interpolation preserves the fine structure of micrometer-scale features while having good smoothness. For nanometer-scale texture dominant regions, local high-frequency filtering is used to extract the height distribution data of the mesh cells from the three-dimensional topography data. High-frequency filtering preserves nanometer-scale texture features and filters out measurement noise. The boundary transition problem between multi-scale feature regions is handled by using a boundary smoothing algorithm to adjust the height distribution data of the mesh cells at the junction. The boundary smoothing algorithm includes: identifying the region boundary, determining the transition band width (usually 3-5 mesh cell widths of the high-density mesh region); applying a weighted interpolation function within the transition band to smoothly transition the mesh density and height data; checking the geometric continuity of the transition band to ensure tangent continuity or higher-order continuity. Finally, a mesh partitioning result with multi-scale adaptability is obtained, and each mesh cell contains accurate height distribution data, providing a solid foundation for subsequent local roughness region analysis and coating growth simulation.

[0144] In this embodiment, based on the geometric feature parameters of the local roughness region, a three-dimensional geometric digital twin sub-model of the local roughness region is constructed, including:

[0145] Based on the geometric feature parameters of the local roughness region, dynamic change features of the local roughness region are extracted, including the morphology evolution trend of the local roughness region in the thin film coating process, which is obtained by analyzing the variation law of the geometric feature parameters with processing time or process conditions;

[0146] According to the dynamic change characteristics, an initial three-dimensional geometric digital twin sub-model of the local rough area is constructed, and the initial three-dimensional geometric digital twin sub-model is generated by parameterized modeling of the geometric feature parameters of the local rough area, wherein the parameterized modeling includes fitting the boundary profile of the local rough area using a spline curve;

[0147] Based on the topography evolution trend, the initial three-dimensional geometric digital twin sub-model is updated to obtain a three-dimensional geometric digital twin sub-model of the local rough area, wherein the updating includes adjusting the control point positions of the spline curve according to the changes in processing time or process conditions to reflect the topography changes of the local rough area in the thin film coating processing.

[0148] In this embodiment, the geometric feature parameters (average depth, maximum width, local normal angle distribution variance) of the local rough area are analyzed in depth, and the variation law thereof under different conditions is extracted. The topography data of the local rough area at different time points in the thin film coating processing is collected or simulated, and a time series database of the geometric feature parameters is established. The variation trend of the geometric feature parameters with the processing time is analyzed, which generally shows that the average depth gradually decreases with the deposition of the coating; the maximum width may slightly increase; the distribution variance of the local normal angle generally decreases, indicating that the surface gradually smoothens. The sensitivity of the geometric feature parameters to the process conditions is analyzed, and the main influencing factors are determined, such as the influence of temperature on surface diffusion, the influence of incident angle on shadow effect, etc.

[0149] Based on the time series data and sensitivity analysis, a dynamic evolution model of the geometric feature parameters is constructed: wherein F0 is an initial feature parameter vector, ΔF is a parameter variation, t is a processing time, and P is a process parameter vector. The dynamic evolution model can be a physical-based analytical model or a data-based machine learning model (such as Gaussian process regression or neural network). Through the dynamic evolution model, the topography evolution trend of the local rough area is extracted, and the geometric feature parameters under different processing times or process conditions are predicted. Based on the initial geometric feature parameters and boundary conditions, a parameterized modeling method is used to construct an initial three-dimensional geometric digital twin sub-model of the local rough area. First, the boundary profile of the local rough area is extracted, and the boundary line between the rough area and the flat area is identified by an edge detection algorithm (such as Canny or Sobel algorithm). Then, the B-spline curve fitting is used on the boundary profile to obtain a continuous and smooth boundary representation. The B-spline curve expression is: where P_i is the i-th control point, N_i,k(u) is the k-th B-spline basis function, u is the parameter value, and n is the total number of control points. Then, according to the average depth and distribution characteristics of the local rough area, a three-dimensional height distribution function h(x, y) is constructed. The height distribution function can be an analytical function (such as a Gaussian function, a polynomial function) or an interpolation function (such as a radial basis function RBF, a thin plate spline TPS). By applying the height distribution function, the surface topography of the initial three-dimensional geometric digital twin sub-model is generated. The model surface can be represented as a parameterized surface S(u, v) = (x(u, v), y(u, v), z(u, v)), where z(u, v) = h(x(u, v), y(u, v)), u and v are variables in the parameter space, used to define points on the surface, usually taking values within a certain interval (such as [0, 1] x [0, 1]); x(u, v) is the x-coordinate of a point on the surface, which is a function of parameters u and v; y(u, v) is the y-coordinate of a point on the surface, which is a function of parameters u and v; z(u, v) is the z-coordinate of a point on the surface, representing the height value, which is a function of parameters u and v; h(x, y) is the height distribution function, which maps the plane coordinates (x, y) to the corresponding height value; h(x(u, v), y(u, v)) represents substituting the plane coordinates generated by parameters u and v into the height function to obtain the corresponding height value.

[0150] Based on the extracted topography evolution trend, the initial three-dimensional geometric digital twin sub-model is dynamically updated. According to the predicted changes in geometric feature parameters, the key control parameters of the model are adjusted. The positions of the control points of the boundary contour B-spline curve are adjusted to reflect the evolution of the boundary shape. The parameters of the height distribution function h(x, y, t) are updated to change with time, such as: h(x, y, t) = h(x, y) x f(t), where f(t) is a time decay function, reflecting the smoothing process of the surface topography. The adjusted three-dimensional geometric digital twin sub-model can reflect the dynamic changes of the local rough area during the coating process, providing a geometric basis for the accurate simulation of coating growth behavior. The dynamically updated three-dimensional geometric digital twin sub-model is meshed to generate a discrete representation suitable for numerical calculation. The meshing process maintains the geometric accuracy of the model while meeting the computational efficiency requirements. The final three-dimensional geometric digital twin sub-model of the local rough area not only accurately represents the static geometric features of the rough area, but also reflects its dynamic evolution behavior during the coating process, providing a high-fidelity geometric basis for subsequent shadow effect calculation and coating growth simulation.

[0151] In this embodiment, the method for obtaining the predicted value of the global coating thickness distribution of the optimized substrate surface includes:

[0152] Based on the optimized process parameter combination, the incident angle distribution and deposition rate distribution of the deposited particles are updated;

[0153] According to the updated incident angle distribution of the deposition particles and the deposition rate distribution, the local coating thickness distribution prediction value of the local rough area and the coating thickness distribution prediction value of the flat area are recalculated;

[0154] The updated local coating thickness distribution prediction value and the coating thickness distribution prediction value of the flat area are integrated to obtain the optimized global coating thickness distribution prediction value of the substrate surface.

[0155] In this embodiment, the optimized process parameter combination P * is substituted into the process parameter digital model to update the key parameters in the model. Based on the updated process parameters, the incident angle distribution P * (θ, φ) of the deposition particles is recalculated. The update of the incident angle distribution takes into account the influence of the process parameters on the deposition source characteristics, such as the change in emission characteristics caused by power changes, the change in scattering characteristics caused by pressure changes, etc. Generally, the optimized incident angle distribution will be more concentrated or directional, which is beneficial to weaken the shadow effect. Based on the updated process parameters, the deposition rate distribution R * (x, y) is recalculated. The update of the deposition rate distribution takes into account the influence of the process parameters on the deposition kinetics, such as the change in adhesion coefficient caused by temperature changes, the change in reaction rate caused by gas flow changes, etc. The optimized deposition rate distribution will generally be more uniform, which is beneficial to improve the coating uniformity. Using the updated incident angle distribution P * (θ, φ) and the deposition rate distribution R * (x, y), the shadow effect calculation for the local rough area is performed again. Due to the change in the incident angle distribution, the shadow effect factor SF * (x, y) will change accordingly, and generally the shadow effect will be weakened, increasing the effective deposition area of the local rough area. Combined with the updated shadow effect factor SF * (x, y) and the deposition rate distribution R * (x, y), the local coating thickness growth rate V of each grid cell in the local rough area is recalculated. According to the optimized local coating thickness growth rate and the processing time t, the local coating thickness distribution prediction value T of the local rough area is recalculated. Similarly, using the updated deposition rate distribution R * (x, y), the coating thickness distribution prediction value T * _f(x, y) of the flat area is recalculated. Since the shadow effect of the flat area is small, its coating thickness is mainly affected by the deposition rate distribution. The local coating thickness distribution prediction value T * _r(x, y) of the local rough area and the coating thickness distribution prediction value T *f(x, y) is integrated according to the spatial position to form a complete global coating thickness distribution prediction value T * (x, y). The integration process needs to handle the transition problem of the region boundary to ensure the continuity of the thickness distribution. The global coating thickness distribution prediction value T * (x, y) is statistically analyzed to calculate the mean μ * , the standard deviation σ * , the maximum value T * _max, the minimum value T * _min and the uniformity index . These statistics intuitively reflect the uniformity of the optimized coating. The optimized global coating thickness distribution prediction value T * (x, y) is visualized in three dimensions to intuitively display the spatial distribution characteristics of the coating. Visualization methods include pseudo-color maps, contour maps, three-dimensional surface maps, etc., which facilitate intuitive understanding and process decision-making. Compare the global coating thickness distribution prediction value and statistics before and after optimization to quantitatively evaluate the optimization effect.

[0156] Evaluation indicators include uniformity improvement rate , thickness range reduction rate , etc., and UI_o is the uniformity index before optimization. Finally, the optimized process parameter combination P * and the expected global coating thickness distribution prediction value T * (x, y) are output as the optimization result, providing scientific guidance for actual production.

[0157] In this embodiment, the local coating thickness growth rate of each grid cell in the local rough area is calculated, including:

[0158] Based on the process parameter digital model, the spatial heterogeneity characteristics of the deposition rate distribution in the thin film coating process are obtained, including the deposition rate gradient distribution and spatial correlation parameters of different regions on the substrate surface.

[0159] According to the three-dimensional geometric digital twin sub-model, the spatial position coordinates of each grid cell in the local rough area on the substrate surface are determined.

[0160] Based on the spatial position coordinates and the spatial heterogeneity characteristics of the deposition rate distribution, the local deposition rate correction value corresponding to each grid cell in the local rough area is calculated, which is the deposition rate adjustment value considering the deposition rate gradient distribution and spatial correlation parameters.

[0161] According to the local deposition rate correction value and the shadow effect factor, the local coating thickness growth rate of each grid cell in the local rough area is calculated, which is the product of the local deposition rate correction value and the shadow effect factor.

[0162] In this embodiment, based on the process parameter digital model, the deposition source characteristics and the cavity geometry of the thin film coating processing equipment are analyzed, and a theoretical distribution model of the deposition rate is constructed. The theoretical model is usually based on the classic cosine law (applicable to thermal evaporation source), modified cosine law (applicable to sputtering source) or other special models (such as PECVD, ALD, etc. process). The reference deposition rate distribution R0(x,y) is obtained from the theoretical model, which reflects the spatial distribution of deposition rate under ideal conditions. Through experimental measurement or high-precision simulation, the actual deposition rate distribution R(x,y) is obtained. By comparing the actual distribution with the reference distribution, the spatial heterogeneity characteristics of the deposition rate distribution are extracted. The spatial gradient vector field of the deposition rate is calculated ∇R(x,y). The gradient vector field describes the rate of change of the deposition rate in each direction, reflecting the directional characteristics of the deposition uniformity. Based on the gradient vector field, a deposition rate gradient distribution map is constructed . The gradient distribution map shows the intensity of the deposition rate change in the form of a scalar field, and the high gradient area represents the fast deposition rate change, and the low gradient area represents the uniform deposition. Using spatial statistical methods, the spatial autocorrelation of the deposition rate distribution is analyzed, and the spatial correlation parameters are calculated. The spatial correlation parameters include the correlation length L_corr (which describes the characteristic scale of the spatial variation of the deposition rate), the anisotropy factor A_aniso (which describes the difference in the change of the deposition rate in different directions), and the spatial autocorrelation function C(r) (which describes the correlation degree of the deposition rate between different distance points). These parameters comprehensively characterize the spatial structure characteristics of the deposition rate distribution. From the three-dimensional geometric digital twin sub-model, the center point coordinates (x_i, y_i, z_i) of each grid cell in the local rough area are extracted. Project the three-dimensional coordinates onto the substrate surface plane to get the two-dimensional spatial position coordinates (x_i, y_i). Considering the height variation of the local rough area, introduce height correction into the spatial position coordinates (x_i, y_i) to get the modified spatial position coordinates (x'_i, y'_i). The correction formula takes into account the geometric relationship between the local surface normal vector and the global deposition direction. Based on the modified spatial position coordinates (x'_i, y'_i) and the spatial heterogeneity characteristics of the deposition rate distribution, the local deposition rate correction value R_mod(x'_i, y'_i) is calculated. The calculation process includes: obtaining the reference deposition rate R0(x'_i, y'_i) at the position (x'_i, y'_i) from the reference deposition rate distribution; based on the deposition rate gradient distribution G(x,y), calculate the rate change ΔR_grad caused by position offset; based on the spatial correlation parameters, calculate the local correlation correction factor C_corr;

[0163] Specifically, ; wherein L corr is the correlation length, describing the characteristic scale of the spatial variation of the deposition rate, r(x, y) is the distance from the current grid cell position (x, y) to the nearest reference point; A an is the anisotropy factor, describing the difference of the deposition rate variation in different directions, θ(x, y) is the azimuth angle from the reference point to position (x, y), θ max is the direction angle of the strongest correlation of the deposition rate.

[0164] ; wherein ∇R(x, y) is the gradient vector of the deposition rate, Δp is the position offset vector.

[0165] Taking the above factors into account, the local deposition rate correction value R mod(x'_i, y'_i) is obtained The local deposition rate correction value accurately reflects the deposition rate distribution after considering spatial heterogeneity, providing accurate input for subsequent coating thickness calculation. Combined with the local deposition rate correction value R mod(x'_i, y'_i) and the shadow effect factor SF(x_i, y_i), the local coating thickness growth rate G(x_i, y_i) of each grid cell is calculated , it is necessary to explain that here i is the index of the corresponding grid cell. The local coating thickness growth rate takes into account the spatial heterogeneity of the deposition rate and the shadow effect caused by surface roughness, and can accurately describe the coating growth behavior in the local rough area. Based on the local coating thickness growth rate G(x_i, y_i), combined with the thin film growth kinetics model (such as island growth, layer growth or mixed growth mode), the time evolution process of the coating thickness is predicted. For simple cases, it can be approximated as a linear integral: ; for complex cases, the dynamic changes of the surface morphology during coating growth need to be considered, and a recursive iteration algorithm is used to calculate the coating thickness.

[0166] The above only describes the preferred embodiments of the present application and is not intended to limit the present application. Although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent replacements to some technical features. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

[0167] The formulas in the specification are dimensionless numerical calculations, the formulas are obtained by software simulation of a large amount of data to obtain a formula closest to the real situation, and the preset parameters and threshold values in the formula are set by those skilled in the art according to the actual situation.

[0168] While embodiments of the application have been shown and described, it is to be understood that the embodiments described are merely divergences of the principles and application of the present application and that numerous modifications, changes, substitutions, and alterations can be made thereto without departing from the spirit and scope of the present application, which is defined by the following claims and their equivalents.

Claims

1. Thin film coating process simulation and optimization system based on digital twin, characterized by: include: A data acquisition module is used to obtain a three-dimensional roughness digital model of the substrate surface and a digital model of the process parameters of the processing equipment during the thin film coating process; A feature recognition and analysis module, configured to identify local rough areas on the substrate surface based on the three-dimensional roughness digital model and extract geometric feature parameters of the local rough areas; a behavior digital twin simulation module, configured to construct a digital twin model of coating growth behavior in the local rough area based on the geometric characteristic parameters of the local rough area and the process parameter digital model, and obtain a predicted value of the local coating thickness distribution in the local rough area during the thin film coating process; a coating uniformity evaluation module, configured to calculate a global coating uniformity index of the substrate surface based on the local coating thickness distribution prediction value and the coating thickness distribution prediction value of a flat area on the substrate surface; a process parameter optimization module, configured to optimize the process parameter digital model based on the global coating uniformity index and a preset uniformity threshold, and obtain an optimized process parameter combination; A coating distribution prediction module is used to update the coating growth behavior digital twin model based on the optimized process parameter combination to obtain a predicted value of the global coating thickness distribution on the optimized substrate surface; The modules are connected via wired and / or wireless means to achieve data transmission between modules.

2. The thin film coating processing simulation and optimization system based on digital twin according to claim 1 is characterized in that: The method for extracting geometric characteristic parameters of the local rough area includes: Meshing the three-dimensional roughness digital model to obtain height distribution data of each mesh unit on the substrate surface; Calculate the angle between the local surface normal vector of each grid cell and the global normal vector of the substrate surface based on the height distribution data, and record it as the local normal deflection angle; According to the local normal deflection angle and a preset deflection angle threshold, identifying grid cells on the substrate surface where the local normal deflection angle is greater than the preset deflection angle threshold, forming the local rough area; The geometric characteristic parameters of the local rough area are extracted, wherein the geometric characteristic parameters include an average depth, a maximum width and a distribution variance of a local normal deviation angle of the local rough area.

3. The thin film coating processing simulation and optimization system based on digital twin according to claim 1 is characterized in that: The method for constructing the digital twin model of coating growth behavior comprises: Constructing a three-dimensional geometric digital twin model of the local rough area according to the geometric characteristic parameters of the local rough area; Obtaining, according to the process parameter digital model, an incident angle distribution and a deposition rate distribution of deposited particles during a thin film coating process; Calculating a shadow effect factor of each grid cell in the local rough area based on the three-dimensional geometric digital twin model and the incident angle distribution of the deposited particles, where the shadow effect factor is the effective area ratio of the grid cell receiving the deposited particles; Calculating a local coating thickness growth rate for each grid cell in the local rough area according to the shadow effect factor and the deposition rate distribution; Based on the local coating thickness growth rate and the processing time, a predicted value of the local coating thickness distribution of the local rough area is obtained.

4. The thin film coating processing simulation and optimization system based on digital twin according to claim 1 is characterized in that: The calculation method of the global coating uniformity index includes: Obtain the predicted value of the coating thickness distribution in the flat area of ​​the substrate surface, which is recorded as the mean coating thickness in the flat area; Calculating the mean coating thickness of the local rough area according to the predicted value of the local coating thickness distribution; Calculating the ratio of the average coating thickness of the local rough area to the average coating thickness of the flat area, and recording it as the local uniformity deviation; The global coating uniformity index is calculated according to the local uniformity deviation and the area ratio of the local rough area to the substrate surface.

5. The thin film coating processing simulation and optimization system based on digital twin according to claim 1 is characterized in that: The method of optimizing the process parameter digital model to obtain an optimized process parameter combination includes: Constructing a mapping relationship model between the process parameter digital model and the global coating uniformity index; Based on the mapping relationship model, a genetic algorithm is used to iteratively optimize the process parameter digital model to obtain a candidate process parameter combination that makes the global coating uniformity index meet the preset uniformity threshold; The feasibility of the candidate process parameter combinations is verified, and infeasible process parameter combinations are eliminated to obtain the optimized process parameter combination.

6. The thin film coating processing simulation and optimization system based on digital twin according to claim 3 is characterized in that: The calculation method of the shadow effect factor includes: Obtaining a local surface normal vector of each grid cell in the local rough area according to the three-dimensional geometric digital twin model of the local rough area; Calculating the angle between the incident direction of the deposited particles and the local surface normal vector according to the incident angle distribution of the deposited particles, and recording it as the incident deflection angle; Calculating a shadow shading ratio of the grid cell according to the incident deflection angle and a height difference between adjacent grid cells in the local rough area; The shadow effect factor of the grid unit is calculated according to the shadow shielding ratio and the incident deflection angle.

7. The thin film coating processing simulation and optimization system based on digital twin according to claim 2 is characterized in that: The step of meshing the three-dimensional roughness digital model to obtain height distribution data of each mesh unit on the substrate surface includes: Scanning the surface of the substrate to obtain three-dimensional topography data of the substrate surface; constructing a feature significance distribution map of the substrate surface based on the three-dimensional topography data; According to the feature significance distribution map, multi-scale feature areas on the surface of the substrate are identified, wherein the multi-scale feature areas include a macro-curvature-dominated area, a micron-scale roughness-dominated area, and a nano-scale texture-dominated area; based on the distribution of the multi-scale feature areas, a dynamic mesh partitioning scheme is generated, wherein the dynamic mesh partitioning scheme comprises: allocating a first mesh density to the macro-curvature-dominated area, wherein a mesh unit size of the first mesh density is determined based on an average curvature radius of the macro-curvature-dominated area; allocating a second mesh density to the micron-scale roughness-dominated area, wherein a mesh unit size of the second mesh density is determined based on a characteristic length of the micron-scale roughness-dominated area, wherein the characteristic length is a spatial scale corresponding to a maximum value of a height change rate in the micron-scale roughness-dominated area; and allocating a third mesh density to the nano-scale texture-dominated area, wherein a mesh unit size of the third mesh density is determined based on a texture period of the nano-scale texture-dominated area, wherein the texture period is obtained by performing a Fourier transform on height distribution data of the nano-scale texture-dominated area; According to the dynamic meshing scheme, meshing the three-dimensional roughness digital model is performed to obtain height distribution data of each mesh unit on the substrate surface, wherein the mesh unit height distribution data of the macro-curvature-dominated region is extracted from the three-dimensional topography data by local quadratic surface fitting, the mesh unit height distribution data of the micron-scale roughness-dominated region is extracted from the three-dimensional topography data by local spline interpolation, and the mesh unit height distribution data of the nano-scale texture-dominated region is extracted from the three-dimensional topography data by local high-frequency filtering; The grid unit height distribution data is subjected to consistency correction, and the consistency correction includes adjusting the grid unit height distribution data at the junction of the macro-curvature dominant region, the micron-level roughness dominant region and the nano-level texture dominant region using a boundary smoothing algorithm based on the boundary transition characteristics between the multi-scale feature regions.

8. The thin film coating processing simulation and optimization system based on digital twin according to claim 3 is characterized in that: The constructing of a three-dimensional geometric digital twin model of the local rough area according to the geometric characteristic parameters of the local rough area includes: Extracting dynamic change characteristics of the local rough area based on geometric characteristic parameters of the local rough area, wherein the dynamic change characteristics include a morphological evolution trend of the local rough area during the thin film coating process; the morphological evolution trend is obtained by analyzing the variation of the geometric characteristic parameters with processing time or process conditions; Constructing an initial three-dimensional geometric digital twin model of the local rough area based on the dynamic change characteristics, wherein the initial three-dimensional geometric digital twin model is generated by performing parametric modeling on the geometric characteristic parameters of the local rough area; wherein the parametric modeling includes fitting the boundary contour of the local rough area using a spline curve; Based on the morphology evolution trend, the initial three-dimensional geometric digital twin model is updated to obtain a three-dimensional geometric digital twin model of the local rough area, wherein the update includes adjusting the control point position of the spline curve according to the change of the processing time or the process conditions.

9. The thin film coating processing simulation and optimization system based on digital twin according to claim 3 is characterized in that: The method for obtaining the predicted value of the global coating thickness distribution of the optimized substrate surface comprises: Based on the optimized process parameter combination, updating the incident angle distribution and deposition rate distribution of the deposited particles; recalculating a predicted value of the local coating thickness distribution of the local rough area and a predicted value of the coating thickness distribution of the flat area according to the updated incident angle distribution and deposition rate distribution of the deposited particles; The updated local coating thickness distribution prediction value is integrated with the coating thickness distribution prediction value of the flat area to obtain the global coating thickness distribution prediction value of the optimized substrate surface.

10. The thin film coating processing simulation and optimization system based on digital twin according to claim 3 is characterized in that: Calculating the local coating thickness growth rate of each grid cell in the local rough area includes: Based on the process parameter digital model, obtaining the spatial heterogeneity characteristics of the deposition rate distribution during the thin film coating process, the spatial heterogeneity characteristics including the deposition rate gradient distribution and spatial correlation parameters in different areas of the substrate surface; determining the spatial position coordinates of each grid unit in the local rough area on the substrate surface according to the three-dimensional geometric digital twin model; Calculating a local deposition rate correction value corresponding to each grid cell in the local rough area based on the spatial position coordinates and the spatial heterogeneity characteristics of the deposition rate distribution; The local coating thickness growth rate of each grid unit in the local rough area is calculated according to the local deposition rate correction value and the shadow effect factor, and the local coating thickness growth rate is the product of the local deposition rate correction value and the shadow effect factor.

Images