Film coating processing technology simulation and optimization system based on digital twinning

Abstract

The invention belongs to the technical field, and discloses a film coating processing technology simulation and optimization system based on digital twinning. The system is composed of a data acquisition module, a feature recognition and analysis module, a behavioral digital twinborn simulation module, a coating uniformity evaluation module, a process parameter optimization module and a coating distribution prediction module. The method comprises the following steps: acquiring a three-dimensional roughness digital model of a base material, accurately identifying a micron-sized local rough region and extracting geometrical characteristic parameters of the micron-sized local rough region; constructing a coating growth behavior digital twinborn model based on the geometrical characteristic parameters and the process parameters, and accurately simulating the coating growth process of the local rough area; and establishing a global coating uniformity quantitative evaluation system, and optimizing a process parameter combination by applying machine learning and a genetic algorithm. The system breaks through the limitation of traditional static shadow hypothesis, real-time tracking and accurate simulation of microstructure evolution in the deposition process are achieved, and the coating growth prediction accuracy is remarkably improved.

Description

Technical Field

[0001] The present invention relates to the technical field, and more particularly, to a thin film coating processing simulation and optimization system based on digital twins. Background Art

[0002] As a key process in modern manufacturing, thin-film coating technology plays an irreplaceable role in fields such as microelectronics, optics, aerospace, and biomedicine. As high-end manufacturing demands ever-increasing product performance, coating uniformity, consistency, and reliability have become core indicators of product quality. However, in actual production, the impact of substrate surface microroughness on thin-film coating quality is becoming increasingly prominent.

[0003] However, the existing thin film coating processing simulation technology faces the "statistical thinking limitation" and cannot accurately reflect the dynamic evolution of the deposition process. In actual production, the deposition behavior of the high aspect ratio groove structure of microelectronic devices, the complex curved surface of optical components, and the precision surface of aerospace coatings all show obvious time-varying characteristics, while the traditional model performs static shadow calculations based on the initial morphology, ignoring the changes in local geometric shape caused by the gradual filling of the coating. This simplified method of "one-time calculation and permanent use" leads to significant deviations between the predicted results and the actual coating distribution in practical applications, especially at the bottom of deep grooves, sharp turning areas, and micro-depressions. In industrial practice, weak bottom areas frequently appear in the deposition of integrated circuit interconnect layers and barrier layers, discontinuities appear in the edge areas of OLED display cathode coatings, and uneven reflectivity of precision optical lens reflective films in areas with changing curvatures. These problems can all be traced back to the prediction deviations caused by the model's inability to track morphological evolution. In addition, existing optimization algorithms rely on empirical parameters and simplified models, and lack the ability to accurately describe dynamic processes at the microscale. This results in strong dependence and poor portability of optimization results. Process engineers have to repeatedly conduct a large number of experimental verifications on each type of new product, which greatly delays the product development cycle and increases R&D costs.

[0004] In view of this, the present invention proposes a thin film coating processing simulation and optimization system based on digital twin to solve the above problems. Summary of the Invention

[0005] In order to overcome the above-mentioned defects of the prior art and achieve the above-mentioned objectives, the present invention provides the following technical solution: a thin film coating processing simulation and optimization system based on digital twin, comprising: 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 is used to identify local rough areas with micron-scale unevenness on the surface of the substrate based on the three-dimensional roughness digital model, and to 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.

[0006] Preferably, 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.

[0007] Preferably, 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, based on 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.

[0008] Preferably, the method for calculating 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.

[0009] Preferably, 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.

[0010] Preferably, the method for calculating 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.

[0011] Preferably, 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.

[0012] Preferably, constructing 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 morphology evolution trend of the local rough area during the thin film coating process; Constructing an initial three-dimensional geometric digital twin model of the local rough area according to the dynamic change characteristics, wherein the initial three-dimensional geometric digital twin model is generated by parameterizing geometric characteristic parameters of the local rough area; 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.

[0013] Preferably, 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.

[0014] Preferably, 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.

[0015] The technical effects and advantages of the thin film coating processing simulation and optimization system based on digital twins of the present invention are as follows: The present invention breaks through the limitations of the traditional static shadow assumption by constructing a dynamically updated digital twin model, and realizes real-time tracking and precise simulation of the evolution of microscopic morphology during the deposition process. The dynamic simulation capability of the present invention 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 of the shadow effect as the morphology changes, the present invention not only optimizes the coating uniformity, but also greatly reduces the generation of local weak areas, reducing the coating structure defect rate to a lower level. At the same time, the efficiency of process parameter optimization is improved, reducing the number of tests and material consumption. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 Schematic diagram of the thin film coating processing simulation and optimization system based on digital twins of the present invention. DETAILED DESCRIPTION

[0017] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0018] This application example provides a thin-film coating process simulation and optimization system based on digital twins. The system's execution entities include: data acquisition equipment, a coating simulation platform, a process optimization server, network transmission equipment, etc., which can be considered as general computing nodes in this application. The coating simulation platform includes but is not limited to: a coating process monitoring system, a process parameter optimization platform, and at least one of the digital twin simulation systems.

[0019] The present invention provides a thin film coating processing simulation and optimization system based on digital twin, comprising the following steps: 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 is used to identify local rough areas with micron-level unevenness on the substrate surface based on a three-dimensional roughness digital model and extract geometric feature parameters of the local rough areas; The behavioral digital twin simulation module is used to construct a digital twin model of coating growth behavior in the local rough area based on the geometric characteristic parameters and process parameter digital model of the local rough area, and obtain the predicted value of the local coating thickness distribution in the local rough area during the thin film coating process; A coating uniformity evaluation module is used to calculate the 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 the flat area of ​​the substrate surface; Process parameter optimization module, which is used to optimize the process parameter digital model based on the global coating uniformity index and the preset uniformity threshold to obtain the optimized process parameter combination; The coating distribution prediction module is used to update the digital twin model of coating growth behavior based on the optimized process parameter combination and obtain the 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.

[0020] The present invention obtains the microstructural 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 to accurately simulate the coating growth process in local rough areas; calculates the global coating uniformity index to deeply evaluate the coating quality; optimizes the process parameter combination based on the uniformity index; and uses the optimized parameters to update the digital twin model, achieving high-precision and high-reliability thin-film coating process optimization.

[0021] In the embodiment of the present invention, see Figure 1 , which is a schematic diagram of a thin film coating process simulation and optimization system based on digital twins of the present invention. In this example, the thin film coating process simulation and optimization system based on digital twins includes: 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; In this embodiment, advanced surface morphology measurement equipment is first deployed, including atomic force microscope (AFM), laser confocal microscope, white light interferometer and other high-precision surface measurement systems to accurately scan the surface of the substrate.

[0022] The scanning range of the substrate surface is usually from 10μm×10μm to 100mm×100mm, and the resolution can reach the nanometer level. The surface height data obtained by scanning is subjected to noise filtering, and Gaussian filtering or median filtering is used to eliminate the measurement noise to obtain high-quality surface morphology raw data. The filtered surface morphology data is three-dimensionally reconstructed to generate a three-dimensional roughness digital model of the substrate surface, which represents the microscopic morphology characteristics of the substrate surface in the form of a 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-10nm / s), substrate temperature (such as room temperature to 800℃), deposition angle (such as 0° to 90°), background pressure (such as 10 -6 to 10 -3Pa), reaction gas flow rate (such as 0-100sccm) and other key process parameters.

[0023] These process parameters are digitized to construct a digital process parameter model. This model includes the numerical ranges of each process parameter, their interrelationships, and their impact on coating growth. The 3D roughness model and the process parameter model are time-stamped and correlated to ensure temporal and spatial consistency, providing fundamental data support for subsequent analysis.

[0024] A feature recognition and analysis module is used to identify local rough areas with micron-level unevenness on the substrate surface based on a three-dimensional roughness digital model and extract geometric feature parameters of the local rough areas; In this embodiment, the three-dimensional roughness digital model is meshed, and the grid unit 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 grid unit is calculated. The local surface normal vector is obtained by calculating the spatial position of the adjacent points, reflecting the local tilt state of the surface. The global normal vector of the substrate surface is determined. The global normal vector usually takes the normal direction of the main plane of the substrate surface and represents the direction of the ideal smooth surface. The angle between the local surface normal vector and the global normal vector of each grid unit is calculated to obtain the local normal deviation angle. The larger the local normal deviation angle, the more serious the surface tilt is and the higher the roughness is. A preset deviation angle threshold is set, usually 5° to 15°, and the specific value is determined according to the coating process and substrate type. Grid units with local normal deviation angles greater than the preset deviation angle threshold are identified, and these units are clustered to form a continuous area, which is a local rough area. Geometric parameters are extracted from identified local roughness regions, including their average depth (the average distance between points within the roughness region and the fitted plane), maximum width (the maximum extension of the roughness region within the surface plane), and the distribution variance of the local normal deflection angle (reflecting the consistency of surface inclination within the roughness region). These geometric parameters comprehensively describe the three-dimensional morphology of the local roughness region and provide key input for subsequent simulations of coating growth behavior.

[0025] The behavioral digital twin simulation module is used to construct a digital twin model of coating growth behavior in the local rough area based on the geometric characteristic parameters and process parameter digital model of the local rough area, and obtain the predicted value of the local coating thickness distribution in the local rough area during the thin film coating process; In this embodiment, a 3D geometric digital twin sub-model of a local roughness region is constructed based on the geometric characteristic parameters of the region. This sub-model accurately replicates the microscopic topography of the local roughness region, including surface undulations, slope variations, and boundary shapes.

[0026] Key parameters of the deposition process are extracted from the digital process parameter model, including the incident angle distribution of the deposited particles and the deposition rate distribution. The incident angle distribution describes the angle range and probability distribution of the deposited particles reaching the substrate surface, which is generally related to the geometry of the deposition equipment and process conditions. The deposition rate distribution describes the number of deposited particles reaching each location on the substrate surface per unit time, reflecting the spatial heterogeneity of the deposition process.

[0027] Based on the three-dimensional geometric digital twin 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 due to the uneven surface, some areas are blocked and cannot receive deposited particles. The shadow effect factor is defined as the effective area ratio of the grid cell to receive deposited particles, and the value range is 0-1. The smaller the value, the stronger the shadow effect. Combining 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 affected by the shadow effect and the local deposition rate, reflecting the uneven growth of the coating on the rough surface. Based on the local coating thickness growth rate and the preset processing time, the local coating thickness distribution prediction value of the local rough area is obtained by cumulative calculation. This prediction value shows the growth of the coating in the local rough area in the form of a three-dimensional distribution, intuitively reflecting the impact of surface roughness on coating uniformity.

[0028] A coating uniformity evaluation module is used to calculate the 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 the flat area of ​​the substrate surface; In this embodiment, the same digital twin model of coating growth behavior is applied to calculate the predicted coating thickness distribution value for the flat area of ​​the substrate surface. Flat areas refer to areas where the local normal deflection angle is less than a preset deflection 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 areas is statistically analyzed, and its average value is calculated, which is recorded as the mean coating thickness of the flat areas. The predicted coating thickness distribution values ​​of the locally rough areas are statistically analyzed, and their average value is calculated, which is recorded as the mean coating thickness of the locally rough areas. The mean coating thickness of the locally rough areas is compared with the mean coating thickness of the flat areas, and the ratio of the two is calculated to obtain the local uniformity deviation. The closer the local uniformity deviation value is to 1, the smaller the difference in coating thickness between the locally rough areas and the flat areas, and the better the coating uniformity; the further it deviates from 1, the greater the difference and the worse the uniformity. Considering the distribution of the locally rough areas on the entire substrate surface, the area proportion of the locally rough areas is calculated. The larger the area proportion, the more significant the impact of the locally rough areas on the global coating uniformity.

[0029] Taking into account the local uniformity deviation and the area ratio of the local rough area, the global coating uniformity index is calculated. This index comprehensively reflects the overall uniformity of the coating thickness distribution on the substrate surface. The specific calculation formula is: Global coating uniformity index Where QU is the global coating uniformity index, JU is the local uniformity deviation, and CU is the area ratio of the local rough area. This index ranges from 0 to 1, with values ​​closer to 1 indicating a more uniform coating and values ​​closer to 0 indicating a more severe coating unevenness.

[0030] Process parameter optimization module, which is used to optimize the process parameter digital model based on the global coating uniformity index and the preset uniformity threshold to obtain the optimized process parameter combination; In this embodiment, a preset uniformity threshold is set, typically between 0.85 and 0.95, with the specific value determined based on the application requirements of the thin film coating. The global coating uniformity index is compared with the preset uniformity threshold. If the global coating uniformity index falls below the preset uniformity threshold, process parameters need to be optimized to improve coating uniformity. Based on a digital twin model of coating growth behavior, a mapping relationship model between process parameters and the global coating uniformity index is constructed. This mapping relationship model is trained using machine learning methods (such as neural networks and support vector machines) and can predict the global coating uniformity index under different process parameter combinations. A genetic algorithm is used to perform a global optimization search for process parameters. In the genetic algorithm, each individual represents a set of process parameter combinations, and the fitness function is the global coating uniformity index. The goal is to maximize the global coating uniformity index or exceed the preset uniformity threshold. An initial population is randomly generated with several process parameter combinations, typically with a population size of 50-200 individuals. Iterative evolution is performed through selection, crossover, and mutation operations, increasing the average fitness of the population generation by generation. The number of iterations is typically 50-500 generations, or until the optimal fitness no longer significantly improves. From the final population, process parameter combinations that meet the preset uniformity threshold are selected as candidate process parameter combinations. Feasibility verification is performed on these candidate process parameter combinations to assess the feasibility of each parameter set on the actual equipment. This feasibility verification includes checking parameter ranges, analyzing the rationality of parameter combinations, and verifying equipment constraints. Infeasible process parameter combinations are eliminated, and the one with the highest global coating uniformity index is selected from the remaining feasible solutions as the optimized process parameter combination.

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

[0032] In this embodiment, the optimized process parameter combination is applied to the digital twin model of coating growth behavior, and the relevant parameters in the model are updated. Based on the optimized process parameters, the incident angle distribution and deposition rate distribution of the deposited particles are recalculated. The optimized process parameters typically change the dynamic 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 shadow effect factor and local coating thickness growth rate of the local rough area are recalculated. The shadow effect factor varies with the incident angle distribution, thus 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. Simultaneously, the coating thickness distribution prediction value of the flat area is recalculated using the same updated parameters. The updated coating thickness distribution prediction value of the local rough area and the coating thickness distribution prediction value of the flat area are integrated to form the optimized global coating thickness distribution prediction value of the substrate surface. The global coating thickness distribution prediction value intuitively 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 expected coating thickness distribution results are output as process guidance for actual production.

[0033] In this embodiment, the method for extracting geometric characteristic parameters of a local rough area includes: Mesh the three-dimensional roughness digital model to obtain the height distribution data of each grid unit on the substrate surface; According to the height distribution data, the angle between the local surface normal vector of each grid cell and the global normal vector of the substrate surface is calculated and recorded as the local normal deflection angle; According to the local normal deflection angle and the preset deflection angle threshold, the grid cells on the substrate surface whose local normal deflection angle is greater than the preset deflection angle threshold are identified to form a local rough area; The geometric characteristic parameters of the local rough area are extracted, and the geometric characteristic parameters include the average depth, maximum width and distribution variance of the local normal deviation angle of the local rough area.

[0034] In this embodiment, a tetrahedron, hexahedron or hybrid meshing algorithm is applied to the three-dimensional roughness digital model to generate a grid structure suitable for surface morphology analysis. The grid unit size is adaptively adjusted according to the complexity of the surface features, with a higher grid density in complex areas and a lower grid density in flat areas. The typical grid unit size range is 50nm-5μm, ensuring that micron-level surface features can be accurately captured. For each grid unit, its node coordinates and height values ​​are obtained to form the height distribution data H(x,y) of the unit. The height data is expressed as a distance perpendicular to the main plane of the substrate, with a positive value indicating a convexity and a negative value indicating a concaveness. The least squares method is used to fit the local surface to the height distribution data within each grid unit and its neighborhood, usually using a second-order or third-order polynomial fitting.

[0035] Based on the fitted surface, calculate the local surface normal vector n_l at the center point of the mesh cell. The local surface normal vector is perpendicular to the tangent plane of the fitted surface and reflects the tilt direction of the surface at that point. Determine the global normal vector n_g of the substrate surface, usually taking the normal direction of the substrate's principal plane (e.g., [0,0,1]). Calculate the angle θ between the local surface normal vector and the global normal vector, i.e., the local normal deviation angle: . The local normal deviation angle ranges from 0° to 90°, and the larger the value, the more severe the surface tilt. Set a preset deviation angle threshold θ_th, which is used to distinguish between flat areas and rough areas. The selection of the threshold needs to take into account material properties, coating process and quality requirements, and is usually in the range of 5° to 15°. Scan the substrate surface to identify grid cells whose local normal deviation angle θ is greater than the preset deviation angle threshold θ_th. Using the connected region analysis algorithm, adjacent grid cells that meet the conditions are clustered to form continuous local rough areas. For each identified local rough area, its geometric characteristic parameters are extracted: the average depth d_avg, calculated as the average distance from all points in the rough area to the fitting plane; the maximum width w_max, which measures the maximum extension size of the rough area in the surface plane; the distribution variance of the local normal deviation angle, which calculates the statistical variance of the local normal deviation angles of all grid cells in the rough area, reflecting the consistency of the surface tilt within the rough area. These geometric characteristic parameters constitute the characteristic vector of the local rough area, providing a quantitative basis for the subsequent simulation of coating growth behavior.

[0036] In this embodiment, the method for constructing a digital twin model of coating growth behavior includes: According to the geometric characteristic parameters of the local rough area, a three-dimensional geometric digital twin model of the local rough area is constructed; According to the digital model of process parameters, the incident angle distribution and deposition rate distribution of deposited particles during the thin film coating process are obtained; Based on the 3D geometric digital twin 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 factor is the effective area ratio of the grid cell receiving the deposited particles. According to the shadow effect factor and deposition rate distribution, the local coating thickness growth rate of each grid cell in the local rough area is calculated; Based on the local coating thickness growth rate and processing time, the local coating thickness distribution prediction value of the local roughness area is obtained.

[0037] In this embodiment, a parametric geometric model of a local roughness region is constructed based on its geometric characteristic parameters (average depth, maximum width, and local normal deviation angle distribution variance). This parametric model, represented by B-spline or NURBS surfaces, accurately reconstructs the three-dimensional topography of the local roughness region. This parametric geometric model is finely meshed, with the mesh density adaptively adjusted based on the local curvature to ensure sufficient mesh accuracy in areas with significant curvature variations.

[0038] The fine mesh model is compared and verified with the original three-dimensional roughness digital model to ensure geometric reconstruction accuracy, typically requiring a root mean square error of less than 10% of the original data resolution. Once verified, a three-dimensional geometric digital twin model of the local roughness region is generated, fully preserving the microscopic morphology of the local roughness region. Based on the digital model of process parameters, the geometry and process conditions of the thin-film coating processing equipment are analyzed to extract the incident characteristics of the deposited particles. A physical model is established to calculate the incident angle distribution P(θ, φ) of the deposited particles, where θ is the polar angle (the angle with the surface normal) and φ is the azimuthal angle. The incident angle distribution is typically related to parameters such as the equipment type, target-substrate distance, and operating pressure.

[0039] Using Monte Carlo methods or analytical models, the deposition rate distribution R(x,y) of the deposited particles at each location on the substrate surface is calculated. This deposition rate distribution takes into account factors such as device geometry, target material size, and particle scattering, and is typically expressed in nm / s. Based on the 3D geometric digital twin model and the incident angle distribution of the deposited particles, a ray tracing analysis is performed on each grid cell within the local roughness region. Virtual rays are emitted at different angles to simulate the path of the deposited particles and detect whether the rays are blocked by surrounding surfaces.

[0040] Count the proportion of incident light that each grid cell can receive and calculate the shadow effect factor SF(x,y): Where P(θ_i,φ_i) is the incidence angle distribution probability, representing the probability of a deposited particle incident from the direction (θ_i,φ_i); SB_i is the shadow occlusion ratio, representing the degree to which particles from the direction (θ_i,φ_i) are occluded; and α_i is the incidence deflection angle, i.e., the angle between the incident direction and the local surface normal. The shadow effect factor ranges from 0 to 1, where 0 represents complete occlusion, 1 represents no occlusion, and intermediate values ​​represent partial occlusion.

[0041] Combine the shadow effect factor SF(x,y) and the deposition rate distribution R(x,y) to calculate the local coating thickness growth rate of each grid cell in the local rough area The local coating thickness growth rate directly reflects the effect of surface roughness on coating growth and reflects the spatial heterogeneity 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 surface morphology does not change significantly during coating growth, the equation can be simplified to T(x,y)=G(x,y)×t. This ultimately yields the predicted coating thickness distribution T(x,y) for the roughened area. This prediction represents the coating's growth in the roughened area in a three-dimensional distribution.

[0042] In this embodiment, 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; Calculate the mean coating thickness of the local rough area based on the predicted value of the local coating thickness distribution; Calculate the ratio of the mean coating thickness in the local rough area to the mean coating thickness in the flat area, which is recorded as the local uniformity deviation; The global coating uniformity index is calculated based on the local uniformity deviation and the area ratio of the local rough area on the substrate surface.

[0043] In this example, a digital twin model of coating growth behavior is applied to identified flat areas of the substrate surface (areas where the local normal deflection angle is less than a preset deflection angle threshold) to calculate the coating thickness distribution T_f(x,y) in the flat areas. Since the surface inclination in the flat areas is relatively small, the shadow effect is weaker, and the coating growth is relatively uniform, but it is still affected by the deposition rate distribution. A statistical analysis is performed on the coating thickness distribution T_f(x,y) in the flat areas to calculate its mean μ_f, standard deviation σ_f, and uniformity index (the ratio of the standard deviation to the mean). The mean coating thickness μ_f in the flat area is used as the reference value, representing the coating thickness on an ideal flat surface. The predicted local coating thickness distribution T_r(x,y) in the local rough area is statistically analyzed to calculate its mean μ_r, standard deviation σ_r, and uniformity index. The mean coating thickness μ_r in the local rough area reflects the overall effect of surface roughness on coating thickness. The ratio of the mean coating thickness in the local rough area to the mean coating thickness in the flat area is calculated to obtain the local uniformity deviation. The local uniformity deviation UD reflects the difference in coating thickness between the rough area and the flat area. UD=1 means the thickness of the two areas is the same, UD<1 means the thickness of the rough area is less than that of the flat area (common in depressions with strong shadow effects), and UD>1 means the thickness of the rough area is greater than that of the flat area (common in local protrusions or edge areas). Statistically analyze the distribution of local rough areas on the entire substrate surface, calculate the total area A_r of the local rough areas and the total area A_t of the substrate surface, and obtain the area ratio. . The area ratio AR reflects the macroscopic distribution characteristics of the surface roughness of the substrate. The larger the AR, the higher the overall surface roughness. Taking into account the local uniformity deviation UD and the area ratio AR, the global coating uniformity index GU is calculated: GU=1-|1-UD|×AR. The value range of the global coating uniformity index GU is 0-1. The closer the value is to 1, the more uniform the coating is, and the closer the value is to 0, the more severe the unevenness is. When the local uniformity deviation UD is close to 1 (the thickness of the rough area is similar to that of the flat area) or the area ratio AR is close to 0 (there are few rough areas), the global uniformity index GU is close to 1, indicating that the global coating uniformity is good.

[0044] In this embodiment, the method for optimizing the process parameter digital model and obtaining the optimized process parameter combination includes: Construct 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 digital model of process parameters to obtain a candidate process parameter combination that makes the global coating uniformity index meet the preset uniformity threshold; The feasibility of candidate process parameter combinations is verified, and infeasible process parameter combinations are eliminated to obtain the optimized process parameter combinations.

[0045] In this embodiment, based on a digital twin model of coating growth behavior, a large number of sample data are collected, representing process parameter combinations and corresponding global coating uniformity indicators. The sampling range covers all dimensions of the process parameter space, ensuring representativeness and diversity of the samples. Common process parameters include deposition rate, substrate temperature, deposition angle, background pressure, and reaction gas flow rate. A machine learning approach is used to construct a mapping relationship model f(P)→GU between process parameters and global coating uniformity indicators, where P is the process parameter vector and GU is the global coating uniformity indicator. Possible machine learning methods include artificial neural networks (ANN), support vector machines (SVM), random forests (RF), or Gaussian process regression (GPR). The trained mapping relationship model is cross-validated to ensure that the model accuracy and generalization ability meet the requirements, typically requiring a root mean square error of less than 5%. A genetic algorithm is used to perform a global optimization search for 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 to make GU exceed a preset uniformity threshold GU_th. An initial population is randomly generated, containing N process parameter combinations. The population size N is typically 50-200. The process parameter values ​​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 index, is calculated.

[0046] A roulette wheel or tournament selection method is used to select individuals for reproduction based on their fitness ratio. Individuals with high fitness have a greater probability of being selected. A crossover operation is performed on the selected individuals to generate new offspring individuals. Crossover methods include single-point crossover, multi-point crossover, or uniform crossover. Mutation operations are performed on some individuals to randomly change the values ​​of certain process parameters to increase population diversity and search capabilities. Elite individuals (a small number of individuals with the highest fitness) are retained and directly enter the next generation to ensure that excellent genes are not lost. The selection, crossover, mutation, and elite retention steps are repeated, and evolution is iterated for multiple generations (usually 50-500 generations) until the termination condition is met. The termination condition can be that the maximum number of iterations is reached, or that the optimal fitness no longer increases significantly over multiple consecutive generations. From the final population, process parameter combinations whose global coating uniformity index GU exceeds the preset uniformity threshold GU_th are screened to form a set of candidate process parameter combinations.

[0047] Conduct feasibility verification on candidate process parameter combinations to assess their feasibility on actual equipment. Feasibility verification includes: parameter range checks to ensure that each parameter is within the permitted range of the equipment; parameter combination rationality analysis to check whether there are conflicts or incompatibilities between parameters; equipment constraint test, such as power limit, temperature upper limit, gas flow range, etc.; process stability assessment to analyze whether the parameter combination will lead to process instability or difficulty in control. Eliminate infeasible process parameter combinations and select one or more combinations with the highest global coating uniformity index from the remaining feasible solutions as the optimized process parameter combination. If there are multiple parameter combinations with similar uniformity indicators, other factors (such as energy consumption, processing efficiency, cost, etc.) can be further considered for the final selection.

[0048] In this embodiment, the calculation method of the shadow effect factor includes: According to the three-dimensional geometric digital twin model of the local rough area, the local surface normal vector of each grid cell in the local rough area is obtained; According to the incident angle distribution of the deposited particles, the angle between the incident direction of the deposited particles and the local surface normal vector is calculated and recorded as the incident deflection angle; Calculate the shadow occlusion ratio of the grid cell based on the incident deflection angle and the height difference between adjacent grid cells in the local rough area; The shadow effect factor of the grid cell is calculated based on the shadow occlusion ratio and the incident deflection angle.

[0049] In this embodiment, the geometric information of each grid cell is extracted from the three-dimensional geometric digital twin model of the local rough area, including node coordinates, cell connection relationships and surface topology. Local surface fitting (such as the least squares method) is applied to each grid cell to obtain the local surface normal vector n_lo(x,y) at the center point of the cell. The local surface normal vector is perpendicular to the tangent plane of the fitted surface, and its direction follows the right-hand rule and points to the outside of the substrate. According to the deposition source characteristics in the digital model of the process parameters, the incident angle distribution P(θ,φ) of the deposited particles is determined. The incident angle distribution describes the probability distribution of deposited particles reaching the substrate surface from all directions, where θ is the polar angle (the 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 based on 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 cell, calculate the angle between the two and obtain the incident deflection angle corresponding to the i-th incident direction . The range of the incident deflection angle α_i is 0° to 180°. When α_i>90°, it means that the incident direction enters from the back, and deposition cannot be performed in this direction. For directions with incident deflection angles α_i≤90°, continue to calculate the shadow occlusion effect. Emit rays from each grid cell in the opposite direction of the incident direction vector d_i to detect whether the ray intersects with the surrounding surface. If the ray intersects with the surrounding surface, it means that the deposited particles in this direction are blocked. For each incident direction d_i, calculate the shadow occlusion ratio SB_i of the grid cell in this direction. If the ray does not intersect with any surface, SB_i=0 (no occlusion); if the ray intersects with the surface, SB_i=1 (complete occlusion).

[0050] Combine the incident deflection angle α_i and the shadow occlusion ratio SB_i to calculate the effective deposition factor of the grid cell in the incident direction The effective deposition factor takes into account the effect of the incident angle on the deposition efficiency (which conforms to the cosine law) and the shadowing effect.

[0051] Calculate the comprehensive shadow effect factor of the grid cell based on the incident angle distribution P(θ_i,φ_i) and the effective deposition factor EF_i The comprehensive shadow effect factor (SF) is normalized to a range of 0-1 and represents the average effective receiving area ratio after accounting for the incident angle distribution. SF = 1 indicates no shadow effect, SF = 0 indicates complete shadowing, and intermediate values ​​indicate partial shadowing. The shadow effect factor (SF) directly affects the local coating growth rate and is a key parameter in the digital twin model of coating growth behavior.

[0052] In this embodiment, the three-dimensional roughness digital model is meshed to obtain the height distribution data of each mesh unit on the substrate surface, including: Scanning the substrate surface with an atomic force microscope or a laser scanning microscope to obtain three-dimensional topography data of the substrate surface; Based on the three-dimensional topography data, a feature significance distribution map of the substrate surface is constructed. The feature significance distribution map is obtained by calculating the local curvature and height change rate of each measurement point in the three-dimensional topography data, where 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 sum of the squares of the height differences between the measurement point and its neighboring measurement points; Identifying multi-scale feature regions on the substrate surface based on the feature significance distribution map, wherein the multi-scale feature regions include a macro-curvature-dominated region, a micron-scale roughness-dominated region, and a nano-scale texture-dominated region. The identification method includes performing threshold segmentation on the feature significance distribution map, wherein the macro-curvature-dominated region corresponds to a low significance threshold range, the micron-scale roughness-dominated region corresponds to a medium significance threshold range, and the nano-scale texture-dominated region corresponds to a high significance threshold range. Based on the distribution of multi-scale feature regions, a dynamic mesh division scheme is generated. The dynamic mesh division scheme includes: allocating a first mesh density to the macro-curvature-dominated region, where the mesh unit size of the first mesh density is determined based on the average curvature radius of the macro-curvature-dominated region; allocating a second mesh density to the micron-scale roughness-dominated region, where the mesh unit size of the second mesh density is determined based on the characteristic length of the micron-scale roughness-dominated region, where the characteristic length is the spatial scale corresponding to the maximum value of the height change rate in the micron-scale roughness-dominated region; allocating a third mesh density to the nano-scale texture-dominated region, where the mesh unit size of the third mesh density is determined based on the texture period of the nano-scale texture-dominated region, where the texture period is obtained by Fourier transforming the height distribution data of the nano-scale texture-dominated region; According to the dynamic meshing scheme, the 3D roughness digital model is meshed to obtain the height distribution data of each mesh unit on the substrate surface. The height distribution data of the mesh units in the macro-curvature-dominated area is extracted from the 3D topography data by local quadratic surface fitting, the height distribution data of the mesh units in the micron-scale roughness-dominated area is extracted from the 3D topography data by local spline interpolation, and the height distribution data of the mesh units in the nano-scale texture-dominated area is extracted from the 3D topography data by local high-frequency filtering. The grid cell height distribution data is corrected for consistency. The consistency correction includes adjusting the grid cell 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 based on the boundary transition characteristics between multi-scale feature regions using a boundary smoothing algorithm to eliminate the influence of the discontinuity of the grid density between regions on the height distribution data.

[0053] In this embodiment, a suitable surface measurement device is selected to precisely scan the substrate surface. For micron-scale features, a laser scanning microscope (LSM) or white light interferometer (WLI) is used, with a scanning range typically ranging from 100μm × 100μm to 10mm × 10mm and a vertical resolution of 10-100nm. For nanoscale features, an atomic force microscope (AFM) is used, with a scanning range typically ranging from 1μm × 1μm to 50μm × 50μm and a vertical resolution of 0.1-1nm. Multi-scale features can be measured using multiple devices or a multi-resolution scanning strategy to ensure complete capture of features at different scales. The raw scan data is preprocessed, including noise filtering, defect repair, and benchmark 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), and the local curvature and height change rate of each measurement point (x,y) are calculated.

[0054] 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 geometric properties. The height change 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 N measurement points in its neighborhood: ; where Z(x_i,y_i) represents the height value of the i-th point in the neighborhood, that is, the height value of the i-th measurement point adjacent to (x,y), and Z(x,y) is the height value of the currently investigated point.

[0055] Combining the local curvature k(x,y) and the height change rate g(x,y), a feature significance distribution map of the substrate surface is constructed , where w1 and w2 are weight coefficients used to balance the contributions of curvature and height change rate. The feature significance distribution map S(x,y) reflects the geometric complexity of each region of the surface, and the higher the value, the more significant the geometric features.

[0056] Perform multi-threshold segmentation on the feature significance distribution map S(x,y) to identify multi-scale feature regions. Set a low significance threshold S_low, a medium significance threshold S_mid, and a high significance threshold S_high, and divide the substrate surface into three types of regions: a macro curvature-dominated region (S(x,y) < S_low), where the geometric change is slow, mainly manifested as large-scale curvature changes; a micro-roughness-dominated region (S_low ≤ S(x,y) < S_high), which has obvious micro-scale concavo-convex features; a nano-scale texture-dominated region (S(x,y) ≥ S_high), which has high-frequency nano-scale texture features.

[0057] Based on the distribution of multi-scale feature regions, design a dynamic grid division scheme. For the macro curvature-dominated region, assign the first grid density, and the grid cell size , where r_curv is the average curvature radius of this region, and k1 is a safety factor (usually 5 - 10). For the micro-roughness-dominated region, assign the second grid density, and the grid cell size , where L_feat is the characteristic length of this region. The characteristic length is defined as the spatial scale corresponding to the maximum value of the height change rate, and k2 is a safety factor (usually 5 - 10). For the nano-scale texture-dominated region, assign the third grid density, and the grid cell size , where T_texture is the texture period of this region, which is obtained by performing a two-dimensional fast Fourier transform (2D-FFT) on the height distribution data of this region and analyzing the power spectral density, and k3 is a safety factor (usually 5 - 10).

[0058] The 3D roughness digital model is meshed according to a dynamic meshing scheme. For regions dominated by macro-curvature, local quadratic surface fitting is used to extract grid cell height distribution data from the 3D topography data. For regions dominated by micron-scale roughness, local spline interpolation is used to extract grid cell height distribution data from the 3D topography data. Spline interpolation preserves the fine structure of micron-scale features while also providing good smoothness. For regions dominated by nano-scale texture, local high-frequency filtering is used to extract grid cell height distribution data from the 3D topography data. High-frequency filtering preserves nano-scale texture features and removes measurement noise. To address boundary transitions between regions with multi-scale features, a boundary smoothing algorithm is used to adjust the grid cell height distribution data at the intersection. This boundary smoothing algorithm includes identifying region boundaries and determining the width of the transition zone (typically 3-5 grid cells in the high-density mesh region); applying a weighted interpolation function within the transition zone to ensure a smooth transition between mesh density and height data; and verifying the geometric continuity of the transition zone to ensure tangent continuity or higher-order continuity. Finally, a grid division result with multi-scale adaptability is obtained, and each grid cell contains accurate height distribution data, providing a solid foundation for subsequent local rough area analysis and coating growth simulation.

[0059] In this embodiment, a three-dimensional geometric digital twin model of the local rough area is constructed based on the geometric characteristic parameters of the local rough area, including: Based on the geometric characteristic parameters of the local rough area, the dynamic change characteristics of the local rough area are extracted. The dynamic change characteristics include the morphological evolution trend of the local rough area during the thin film coating process. The morphological evolution trend is obtained by analyzing the change pattern of the geometric characteristic parameters with processing time or process conditions. Based on the dynamic change characteristics, an initial three-dimensional geometric digital twin model of the local rough area is constructed. The initial three-dimensional geometric digital twin model is generated by parametric modeling of the geometric characteristic parameters of the local rough area, wherein the parametric modeling includes fitting the boundary contour of the local rough area with a spline curve; Based on the morphology evolution trend, the initial three-dimensional geometric digital twin model is updated to obtain the three-dimensional geometric digital twin model of the local rough area. The update includes adjusting the control point position of the spline curve according to the changes in processing time or process conditions to reflect the morphology changes of the local rough area during the thin film coating process.

[0060] In this example, the geometric characteristic parameters (average depth, maximum width, and local normal angle distribution variance) of the locally roughened areas are analyzed in depth to extract their variations under different conditions. Topographic data of the locally roughened areas at different time points during the thin-film coating process are collected or simulated to establish a time series database of the geometric characteristic parameters. Analysis of the changing trends of the geometric characteristic parameters over processing time typically reveals the following: the average depth gradually decreases with coating deposition; the maximum width may slightly increase; and the distribution variance of the local normal angle generally decreases, indicating a gradual smoothing of the surface. The sensitivity of the geometric characteristic parameters to process conditions is analyzed to identify key influencing factors, such as the effect of temperature on surface diffusion and the influence of the angle of incidence on shadowing.

[0061] Based on time series data and sensitivity analysis, a dynamic evolution model of geometric characteristic parameters is constructed: , where F0 is the initial feature parameter vector, ΔF is the parameter change, t is the processing time, and P is the process parameter vector. The dynamic evolution model can be a physics-based analytical model or a data-based machine learning model (such as Gaussian process regression or neural network). Through the dynamic evolution model, the morphology 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 parametric modeling method is used to construct the initial three-dimensional geometric digital twin model of the local rough area. First, the boundary contour 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 boundary contour is fitted with a B-spline curve to obtain a continuous and smooth boundary representation. The expression of the B-spline curve is: , where P_i is the i-th control point, N_i,k(u) is the k-th order B-spline basis function, u is the parameter value, and n is the total number of control points. Next, a 3D height distribution function h(x,y) is constructed based on the average depth and distribution characteristics of the local roughness region. The height distribution function can be an analytical function (such as a Gaussian function or a polynomial function) or an interpolation function (such as a radial basis function (RBF) or a thin plate spline (TPS). The height distribution function is applied to generate the surface topography of the initial 3D geometric digital twin model. 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, and usually take values ​​in a certain interval (such as [0,1]×[0,1]); x(u,v) is the x-coordinate of the point on the surface, which is a function of the parameters u and v, y(u,v) is the y-coordinate of the point on the surface, which is a function of the parameters u and v, z(u,v) is the z-coordinate of the point on the surface, which represents the height value, which is a function of the parameters u and v, h(x,y) is the height distribution function, which maps the plane coordinates (x,y) to the corresponding height values, and h(x(u,v),y(u,v)) means substituting the plane coordinates generated by the parameters u and v into the height function to obtain the corresponding height value.

[0062] Based on the extracted morphology evolution trends, the initial 3D geometric digital twin model is dynamically updated. Key control parameters of the model are adjusted based on predicted changes in geometric feature parameters. The control points of the boundary contour B-spline curve are repositioned to reflect the evolution of the boundary shape. The parameters of the height distribution function h(x,y,t) are updated to change over time, such as: h(x,y,t)=h(x,y)×f(t), where f(t) is a time-decay function that reflects the smoothing process of the surface morphology. The adjusted 3D geometric digital twin model can reflect the dynamic changes in local roughness during thin-film coating processing, providing a geometric foundation for accurate simulation of coating growth behavior. The dynamically updated 3D geometric digital twin model is meshed to generate a discrete representation suitable for numerical calculations. The meshing process maintains the geometric accuracy of the model while meeting computational efficiency requirements. The resulting three-dimensional geometric digital twin model of the local rough area not only accurately expresses the static geometric characteristics of the rough area, but also reflects its dynamic evolution behavior during the coating processing, providing a high-fidelity geometric basis for subsequent shadow effect calculations and coating growth simulations.

[0063] In this embodiment, the method for obtaining the predicted value of the global coating thickness distribution of the optimized substrate surface includes: Based on the optimized process parameter combination, the incident angle distribution and deposition rate distribution of the deposited particles are updated; Recalculating the predicted value of the local coating thickness distribution in the rough area and the predicted value of the coating thickness distribution in 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 optimized global coating thickness distribution prediction value of the substrate surface.

[0064] In this embodiment, the optimized process parameter combination P * Substitute the process parameters into the digital model and update the key parameters in the model. Based on the updated process parameters, recalculate the incident angle distribution P of the deposited particles * (θ, φ). The update of the incident angle distribution takes into account the impact of process parameters on the characteristics of the deposition source, such as changes in emission characteristics caused by power changes, changes in scattering characteristics caused by pressure changes, etc. Generally, the optimized incident angle distribution will be more concentrated or directional, which is conducive to reducing the shadow effect. Based on the updated process parameters, the deposition rate distribution R is recalculated. * (x,y). The update of the deposition rate distribution takes into account the influence of process parameters on deposition dynamics, such as changes in adhesion coefficient caused by temperature changes, changes in reaction rate caused by changes in gas flow, etc. The optimized deposition rate distribution is usually more uniform, which is conducive to improving coating uniformity. Using the updated incident angle distribution P * (θ,φ) and the deposition rate distribution R * (x,y), recalculate the shadow effect of the local rough area. Due to the change of the incident angle distribution, the shadow effect factor SF * (x,y) will change accordingly, usually the shadow effect will be weakened, so that the effective deposition area of ​​the local rough area will increase. Combined with the updated shadow effect factor SF * (x,y) and the deposition rate distribution R * (x,y), recalculate the local coating thickness growth rate for each grid cell within the local roughness area According to the optimized local coating thickness growth rate and processing time t, the local coating thickness distribution prediction value of the local rough area is recalculated. Similarly, using the updated deposition rate distribution R * (x,y), recalculate the coating thickness distribution prediction value T in the flat area * _f(x,y). 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 of the local rough area is * _r(x,y) and the coating thickness distribution prediction value T in the flat area * _f(x,y) is integrated according to the spatial position to form a complete prediction value T of the global coating thickness distribution on the substrate surface* (x,y). The integration process needs to deal with the transition problem of the regional boundary to ensure the continuity of the thickness distribution. The predicted value of the global coating thickness distribution T * (x,y) performs statistical analysis and calculates the mean μ * , standard deviation σ * , maximum value T * _max, minimum value T * _min and uniformity index These statistics intuitively reflect the uniformity of the optimized coating. The optimized global coating thickness distribution prediction value T * 3D visualization of the (x,y) coordinates provides an intuitive display of the coating's spatial distribution. Visualization options include pseudo-color, contour plots, and 3D surface plots, facilitating intuitive understanding and process decision-making. Compare global coating thickness distribution predictions and statistics before and after optimization to quantitatively evaluate optimization results.

[0065] Evaluation indicators include uniformity improvement rate , thickness range reduction rate etc., 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) is output as the optimization result to provide scientific guidance for actual production.

[0066] In this embodiment, calculating the local coating thickness growth rate of each grid cell in the local rough area includes: Based on the digital model of process parameters, the spatial heterogeneity characteristics of the deposition rate distribution during the thin film coating process are obtained. The spatial heterogeneity characteristics include the deposition rate gradient distribution and spatial correlation parameters in different areas of the substrate surface. Determine the spatial position coordinates of each grid cell on the substrate surface within the local rough area based on the three-dimensional geometric digital twin model; Based on the spatial heterogeneity of spatial position coordinates and sedimentation rate distribution, the local sedimentation rate correction value corresponding to each grid cell in the local rough area is calculated. The local sedimentation rate correction value is the sedimentation rate adjustment value after considering the sedimentation rate gradient distribution and spatial correlation parameters. According to the local deposition rate correction value and the shadow effect factor, the local coating thickness growth rate of each grid unit in the local rough area is calculated. The local coating thickness growth rate is the product of the local deposition rate correction value and the shadow effect factor.

[0067] In this embodiment, based on the digital model of process parameters, the deposition source characteristics and cavity geometry of the thin film coating processing equipment are analyzed to construct a theoretical distribution model of the deposition rate. The theoretical model is usually based on the classical cosine law (applicable to thermal evaporation sources), the modified cosine law (applicable to sputtering sources) or other special models (such as PECVD, ALD and other processes). The baseline deposition rate distribution R0(x,y) is obtained from the theoretical model, which reflects the spatial distribution of the deposition rate under ideal conditions. The actual deposition rate distribution R(x,y) is obtained through experimental measurement or high-precision simulation. The actual distribution is compared with the baseline distribution to extract the spatial heterogeneity characteristics of the deposition rate distribution. The spatial gradient vector field ∇R(x,y) of the deposition rate is calculated. The gradient vector field describes the rate of change of the deposition rate in each direction and reflects the directional characteristics of the deposition uniformity. Based on the gradient vector field, a deposition rate gradient distribution map is constructed. The gradient distribution diagram displays the severity of the deposition rate change in the form of a scalar field. High-gradient areas indicate rapid deposition rate changes, while low-gradient areas indicate more uniform deposition. Spatial statistical methods are used to analyze the spatial autocorrelation of the deposition rate distribution and calculate the spatial correlation parameters. The spatial correlation parameters include the correlation length L_corr (describing the characteristic scale of the spatial variation of the deposition rate), the anisotropy factor A_aniso (describing the difference in the deposition rate changes in different directions), and the spatial autocorrelation function C(r) (describing the degree of correlation between the deposition rates at different distances). These parameters comprehensively characterize the spatial structural characteristics of the deposition rate distribution. From the three-dimensional geometric digital twin model, the center point coordinates (x_i, y_i, z_i) of each grid cell in the local rough area are extracted. The three-dimensional coordinates are projected onto the substrate surface plane to obtain the two-dimensional spatial position coordinates (x_i, y_i). Considering the height variation of the local rough area, a height correction is introduced into the spatial position coordinates to obtain the corrected 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 corrected 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; calculating the rate change ΔR_grad caused by the position offset based on the deposition rate gradient distribution G(x, y); and calculating the local correlation correction factor C_corr based on the spatial correlation parameter; Specifically, ; Where L_corr is the correlation length, which describes the characteristic scale of the spatial variation of the sedimentation 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, which describes the difference in the variation of the sedimentation rate in different directions, θ(x,y) is the azimuth from the reference point to the position (x,y), and θ_max is the direction angle with the strongest correlation of the sedimentation rate.

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

[0069] Combining the above factors, the local deposition rate correction value 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. Combining 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 of each grid cell is calculated. , it should be explained that i here is the index of the corresponding grid cell. The local coating thickness growth rate comprehensively considers the spatial heterogeneity of the deposition rate and the shadow effect caused by the 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 film growth dynamics model (such as island growth, layered growth or mixed growth mode), the time evolution of the coating thickness is predicted. For simple cases, it can be approximated by a linear integral: ; For complex situations, it is necessary to consider the dynamic changes of surface morphology during coating growth, and use a recursive iterative algorithm to calculate the coating thickness.

[0070] The foregoing description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art will be able to modify the technical solutions described in the foregoing embodiments or to substitute equivalents for some of the technical features. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention shall be included within the scope of protection of the present invention.

[0071] The formulas in this manual are all dimensionless and calculated using numerical values. The formulas are obtained by collecting a large amount of data and performing software simulation to obtain the most recent real situation. The preset parameters and thresholds in the formulas are set by technicians in this field based on actual conditions.

[0072] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to the embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the 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