Detailed Description
The technical solutions of the present invention will be clearly and completely described in connection with the embodiments, and it is obvious that the described embodiments are only some embodiments of the present invention, not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
As shown in fig. 1, the control method of the composite copper foil thickness control apparatus of the present embodiment specifically may include:
Step S101, constructing a multi-physical field coupling model describing the production process of the composite copper foil according to the mutual coupling relation of the electric field, the flow field and the mass transfer process of the electroplating process, integrating the interface layer morphology and the thickness parameter of the composite copper foil into the multi-physical field coupling model, and acquiring electric field distribution, electrolyte flow state and copper ion concentration distribution data on interface layers with different shapes through simulation of the multi-physical field coupling model.
Dividing a quadrilateral grid unit on the surface of a composite copper foil substrate by adopting a finite element numerical calculation method, calculating according to three-dimensional coordinates of boundary points of the grid unit to obtain an interface layer thickness distribution function, establishing a composite copper foil three-dimensional geometric structure digital model according to the interface layer thickness distribution function, carrying out discrete solution on a Navier-Stokes equation by adopting a finite volume method to obtain electrolyte velocity field and pressure field distribution data near the interface layer, solving a Poisson equation according to the electrolyte velocity field and pressure field distribution data by adopting a finite element method to obtain potential distribution, calculating by combining an electrolyte conductivity distribution function to obtain current density at the interface layer, establishing a copper ion concentration control equation according to the current density at the interface layer, and solving a mass transfer equation by adopting a finite difference method to obtain ion concentration distribution data.
Specifically, a quadrilateral grid unit is divided on the surface of the composite copper foil substrate by adopting finite element numerical calculation, an interface layer thickness distribution function is calculated according to the three-dimensional coordinates of boundary points of the grid unit, and a composite copper foil three-dimensional geometric structure digital model is constructed through substrate surface roughness parameters and interface layer thickness data. And introducing dynamic parameters such as electrolyte density, viscosity, temperature and the like, performing discrete solution on a Navier-Stokes equation by adopting a finite volume method, and acquiring distribution data of an electrolyte speed field and a pressure field near an interface layer by taking a composite copper foil three-dimensional geometric structure digital model as a boundary condition. Based on the three-dimensional geometry of the composite copper foil and the electrolyte flow field data, a Poisson equation is solved by adopting a finite element method to calculate potential distribution, an electric field intensity distribution is calculated by combining an electrolyte conductivity distribution function, and a current density calculation expression at an interface layer is established. And (3) adopting current density distribution at the interface layer as an electrochemical reaction source term, establishing a copper ion concentration control equation, and solving a mass transfer equation by using a finite difference method by combining a speed field transport term and a diffusion coefficient diffusion term. And solving a multi-physical field coupling equation set by adopting a Newton iteration method according to the electric field intensity, the fluid speed and the ion concentration distribution data, and obtaining a steady-state solution as a multi-physical field coupling model in the electroplating process of the composite copper foil. Substituting the multi-physical field coupling model into different interface layer morphology parameters, and adopting a Longg-Kutta method to carry out time-driven solution on the coupling equation set so as to obtain evolution data of each physical field along with time. In the production process of the composite copper foil, key factors influencing the electroplating quality are interface layer structure and multi-physical field distribution. Taking the interface layer thickness distribution function h (x) as An example, when the surface of the composite copper foil substrate presents a periodic coarse structure, the fourier series expansion can be used for describing the height distribution of h (x) =h0+Σ (Ancos (2pi nx/L) + Bnsin (2pi nx/L)), wherein h0 is the average thickness, an and Bn are fourier coefficients, L is the period length, and x is the position parameter of different interface layer morphologies. And (3) performing least square fitting on the actually measured contour data to obtain a first 5 th order Fourier coefficient, and establishing an accurate three-dimensional geometric model. After the geometric model is obtained, the calculation of the electrolyte flow field is particularly important. For a typical laminar flow regime in a plating bath, the reynolds number re=ρvd/μ, where ρ is the electrolyte density, v is the characteristic velocity, D is the hydraulic diameter, and μ is the dynamic viscosity. The Navie-Stokes equation is discretized by a finite volume method, the mesh size is taken to be 1/10 of the thickness of the interface layer, namely 1 micron, and the pressure-speed coupling equation is solved by using a SIMPLE algorithm to obtain the three-dimensional speed field distribution near the interface layer. The electric field distribution is calculated based on poisson equation for potential: wherein the conductivity kappa of the electrolyte changes with the ion concentration c, wherein kappa=k0 (1+alpha (c-c 0)), wherein kappa 0 is the standard state conductivity of 5.8S/m, c0 is the standard ion concentration of 0.6mol/L, and alpha is the concentration coefficient of 0.25L/mol. Applying potential-0.3V on the surface of the cathode and potential-0V on the surface of the anode, solving to obtain potential distribution phi by a finite element method, and further calculating the electric field intensity The mass transfer process is controlled by the convection-diffusion equation: Wherein the diffusion coefficient d=6.5×10 -10m2/s. At the interface layer, the current density and the ion flux satisfy the relation that j= nFN, N is the valence number 2,F, the Faraday constant 96485C/mol, and N is the ion flux mol/(m 2. S). And adopting a finite difference format of an interlaced grid, taking the time step length as 0.001s, keeping the space step length consistent with the flow field calculation, and solving to obtain the concentration field evolution. When solving the multi-physical field coupling equation set, selecting the relative error margin as 10 -6, and adopting a Newton iteration method to solve. For a given interface layer topography, such as a periodic trench structure with a 30 degree pitch, 75 micron pitch, calculations indicate a 25% decrease in current density at the trench bottom and a 35% increase in current density at the edge, consistent with the experimentally observed trend of the coating thickness distribution. The dynamic evolution characteristics of each physical field can be obtained by performing time propulsion by a fourth-order Dragon-Kutta method, wherein the time step is 0.01s and the calculation period is 10 s.
Step S102, judging whether an electric field local distortion phenomenon exists according to electric field distribution data obtained through simulation, if so, determining the position coordinates of a distortion area, marking the electric field data of the distortion area as abnormal data, meanwhile, judging whether a flow disturbance phenomenon exists according to electrolyte flow state data, if so, determining the position coordinates of a disturbance area, and marking the flow field data of the disturbance area as abnormal data.
According to grid point electric field distribution data, calculating an electric field intensity gradient value, if the gradient value exceeds three times of a local average value, obtaining a distorted area boundary point coordinate set by adopting a three-dimensional space clustering method, calculating a speed gradient tensor and a rotation field from a grid point speed vector field, if the rotation value is larger than twice of a laminar flow reference value, obtaining a disturbance area boundary point coordinate set by adopting a boundary tracking algorithm, carrying out grid point matching on the distorted area boundary point coordinate set and the disturbance area boundary point coordinate set to obtain an overlapped area volume fraction, if the volume fraction is larger than a critical value, judging the area as a strong coupling abnormal area, obtaining strong coupling abnormal area temperature field distribution by adopting a thermocouple array, and obtaining local conductivity according to the temperature field distribution and the conductivity temperature coefficient.
Specifically, the gradient value g (x, y, z) of the electric field intensity of each grid point is calculated according to the electric field distribution data, if the gradient value g (x, y, z) exceeds three times of the standard deviation of the local average value, a three-dimensional space clustering method is adopted to determine the coordinate set P (x, y, z) of the boundary point of the distorted region, and a distorted region position marking matrix M1 is generated. And calculating a velocity gradient tensor and a rotation field omega (x, y, z) from the velocity vector field v (x, y, z) of the fluid grid point, and if the rotation value is more than twice the laminar flow reference value or the velocity gradient exceeds the Reynolds stress threshold value, determining a disturbance zone boundary point coordinate set Q (x, y, z) by adopting a boundary tracking algorithm, and generating a disturbance zone position mark matrix M2. And carrying out grid point matching on the electric field distortion region marking matrix M1 and the flow field disturbance region marking matrix M2, calculating the volume fraction alpha (x, y, z) of the overlapped region, and marking the region as a strong coupling abnormal region if the volume fraction alpha (x, y, z) is larger than a critical value. And acquiring abnormal region temperature field distribution T (x, y, z) by adopting a thermocouple array, calculating local conductivity kappa (x, y, z) by combining the conductivity temperature coefficient, and generating temperature compensated electric field intensity distribution E (x, y, z). And extracting the electric field intensity E (x, y, z), the fluid velocity v (x, y, z) and the temperature T (x, y, z) data of the abnormal region according to the distortion region position marking matrix M1 and the disturbance region position marking matrix M2, and calculating the instantaneous frequency characteristic of the abnormal data by adopting Hilbert transformation. Wavelet decomposition is carried out on the abnormal region data, abnormal characteristic coefficients under different scales are extracted, multi-scale abnormal characteristic vectors are generated, and an abnormal region characteristic database is established. In the electric field distortion detection process, the electric field intensity gradient is calculated by adopting a central differential format:
g (x, y, z) = [ (E (x+h, y, z) -E (x-h, y, z))/2 h, (E (x, y+h, z) -E (x, y-h, z))/2 h, (E (x, y, z+h) -E (x, y, z-h))/2 h ], where h is the grid spacing value of 0.1mm. For a typical plating bath, when the gradient value of the localized region exceeds the mean value by a factor of 3 standard deviation, i.e., |g (x, y, z) | > μg+3σg, this point is labeled as the distortion point, where μg is the gradient field mean value of 0.5kV/m 2 and σg is the standard deviation of 0.15kV/m 2. In flow field disturbance judgment, velocity gradient tensor The calculation of (1) also adopts a central difference method, and the rotation fieldThe local vortex strength was characterized. Under laminar flow conditions, the reference curl value is 10s-1, and when the actual curl exceeds 20s-1, a significant disturbance is indicated. The Reynolds stress threshold was set at 0.1Pa, corresponding to a turbulence intensity of 5%. The region matching process adopts grid overlapping degree analysis, and the volume fraction alpha (x, y, z) represents the proportion of simultaneous electric field distortion and flow field disturbance in unit volume. When the α value exceeds 0.3, it is indicated that there is strong coupling abnormality in this region. In actual plating, a typical strongly coupled anomaly size is about 5mm by 2mm. The effect of temperature field distribution on conductivity is described by a linear relationship of κ (T) =κ 0[1+β(T-T0), where κ 0 is the conductivity 5.8S/m at reference temperature T 0 =298K and β is the temperature coefficient 0.02K-1. When the temperature is raised to 308K, the conductivity increases by about 20%, resulting in a corresponding decrease in the local electric field strength. The hilbert transform is used to extract the instantaneous frequency characteristics of the anomaly signal. The electric field intensity data E (t) is transformed to obtain an analysis signal z (t) =E (t) +jH [ E (t) ], wherein H [ E (t) ] is a Hilbert transformation result. The instantaneous frequency f (t) = (1/2 pi) ·d (argz (t))/dt reflects the frequency characteristics of the electric field fluctuation in the abnormal region, with typical values in the range of 0.1-10 Hz. And 4 layers of decomposition are carried out on the abnormal region data by adopting db4 wavelet basis in wavelet decomposition, so as to obtain different scale coefficients. The first coefficient reflects high frequency disturbances (> 5 Hz) and the fourth coefficient reflects low frequency variations (< 0.5 Hz). By analyzing the energy distribution of each scale factor, the dominant frequency component of the disturbance can be identified. When the high frequency component energy ratio exceeds 50%, severe local fluctuations are indicated. The multi-scale feature extraction method can comprehensively describe the dynamic characteristics of the abnormal region and provides a basis for the optimization of the subsequent process parameters.
Step S103, copper ion concentration is clustered according to copper ion concentration distribution data to obtain copper ion groups with different concentration levels, center coordinates of the copper ion groups with different concentration levels are calculated, whether the center coordinates of the copper ion groups with different concentration levels are located in an abnormal area or not is judged according to position coordinates of a distortion area and a disturbance area, and if the copper ion groups with different concentration levels are located in the abnormal area, the copper ion groups with the concentration levels are marked as abnormal concentration groups.
The method comprises the steps of constructing a three-dimensional concentration field according to copper ion concentration distribution data, obtaining three concentration level areas of high, medium and low through a density clustering algorithm, obtaining mass center coordinates and point density according to the concentration level areas, obtaining a group space distribution function and a group characteristic radius through a Gaussian kernel function according to the mass center coordinates and the point density, constructing an abnormal area space envelope surface by adopting a boundary point set of an electric field distortion area and a fluid disturbance area, obtaining a shortest distance value between the abnormal area space envelope surface and a concentration group mass center, judging whether the ratio of the shortest distance value to the group characteristic radius is smaller than a critical value, and marking the concentration group as an abnormal group and generating an abnormal group space distribution characteristic matrix if the ratio is smaller than the critical value.
Specifically, a three-dimensional concentration field is constructed according to copper ion concentration distribution data c (x, y, z), the concentration field is spatially clustered by adopting a density clustering algorithm based on Euclidean distance, and a concentration gradient threshold value is passedAnd a minimum point number threshold n to obtain three concentration level areas C 1、C2、C3 of high, medium and low. For each concentration level region C, calculating centroid coordinates r= (x, y, z), counting point densities ρ (R) in the region, and calculating a group space distribution function phi (R) and a group feature radius R by using a gaussian kernel function. Based on the boundary point set of the electric field distortion region D 1 and the fluid disturbance region D 2, an abnormal region space envelope surface S (x, y, z) is constructed, and the shortest distance D from the centroid of the concentration group to the envelope surface is calculated. And taking the ratio lambda=d/R of the shortest distance d between the mass center of the concentration group and the abnormal region and the characteristic radius R of the group as a judging parameter, and marking the concentration group as an abnormal group if lambda is smaller than a critical value lambda c. For the concentration group overlapped by the space position, calculating an overlapped volume proportion gamma j, and if gamma j is larger than an overlapping degree threshold gamma c, merging overlapping areas to form a new concentration group, and updating group characteristic parameters. And establishing coding mapping for the marked abnormal concentration group to generate an abnormal group space distribution characteristic matrix M (x, y, z), wherein matrix elements comprise characteristic quantities such as concentration values, centroid positions, group radiuses and the like. In copper ion concentration distribution analysis, a density clustering algorithm distinguishes different concentration level areas by setting a concentration gradient threshold value and a minimum point number. Taking an actual plating tank as an example, when a concentration gradient threshold is setWhen the minimum point number threshold n=100, the concentration field can be divided into a high concentration region (c >0.8 mol/L), a medium concentration region (0.4-0.8 mol/L), and a low concentration region (c <0.4 mol/L). For the divided concentration groups, the centroid coordinate calculation adopts a weighted average method: r=Σ (c jrj)/Σcj, where c j is the concentration value at point j, R j is its position vector: in the abnormal region judgment, if the shortest distance d between the mass center of a certain concentration group and the boundary of an abnormal region is 2mm, and the characteristic radius R is 4mm, a judgment parameter lambda=0.5 is set, a critical value lambdac=0.8, and the concentration group is marked as an abnormal group. When gamma j is more than 0.3, the new concentration group characteristic parameters after combination are obtained through weighted average, the weight is in direct proportion to the number of the original group, the abnormal group characteristic matrix M (x, y, z) adopts a multidimensional data structure, each matrix element comprises information such as concentration value c (x, y, z), centroid position r= (x, y, z), group radius R, point density rho and the like, the evolution process of the abnormal concentration group can be tracked through the matrix, the distortion area of an electric field is observed, the ion concentration tends to be a local accumulation phenomenon, the concentration gradient can reach 0.2 mol/(L.mm), and in a fluid disturbance area, the ion distribution tends to be uniform, and the concentration gradient is reduced to below 0.02 mol/(L.mm). The multi-scale characteristic analysis method can comprehensively describe the concentration distribution characteristics of the abnormal region, and provides important basis for process parameter optimization.
Step S104, for the abnormal concentration group, a time sequence prediction algorithm is adopted, according to the copper ion concentration change trend at the historical moment, the concentration change at the future moment is predicted, the concentration change curves of the different Chang Nongdu groups are obtained, and meanwhile, according to the concentration change curves, a curve fitting algorithm is adopted, and the deposition rate change rule diagram of the different concentration groups is obtained.
The method comprises the steps of carrying out resampling at fixed sampling intervals according to historical data of abnormal concentration groups to obtain time sequence data with equal intervals and a concentration change rate sequence, training a long-time memory network model by adopting the concentration change rate sequence, outputting a concentration predicted value at the next moment through the network model, calculating a concentration difference value at the adjacent moment according to the concentration predicted value, calculating a theoretical deposition thickness increment sequence according to the concentration difference value and Faraday law, carrying out Fourier transformation on the theoretical deposition thickness increment sequence to extract frequency components, carrying out polynomial fitting through the frequency components to obtain a deposition rate change rule function, and if the root mean square error of the polynomial fitting is larger than a preset threshold, increasing a polynomial order and repeating the fitting process.
Specifically, a time sequence { t 1,c1 } is constructed according to the historical data of the abnormal concentration group, the concentration data is resampled by adopting a fixed sampling interval δt, the time sequence data { t, c } with equal intervals is generated, and a first-order concentration change rate sequence { t, Δc/Δt } is obtained through differential operation. Setting the width tau of a time window, modeling a concentration change rate sequence by adopting a long-short-time memory network, inputting the concentration change rate value of tau continuous moments into the network, outputting the concentration change rate value as a concentration predicted value of the next moment, and optimizing network parameters by a back propagation algorithm. And calculating concentration differences at adjacent moments according to the predicted concentration sequence, calculating theoretical deposition thickness increment delta h at each moment based on Faraday law and current density distribution, and generating a cumulative sequence { t, h } of thickness along with time. And carrying out Fourier transformation on the accumulated thickness sequence, extracting a main frequency component omega and a corresponding amplitude A thereof in the amplitude spectrum, and obtaining a deposition rate time function v (t) through inverse transformation reconstruction. And (3) carrying out least square fitting on the deposition rate v (t) by adopting an n-order polynomial function, obtaining a polynomial coefficient { a }, and calculating the root mean square error sigma of the fitting curve and the measured data. And according to the comparison result of the root mean square error sigma and the preset threshold value, adjusting the polynomial order n and repeating the fitting process until the error meets the requirement, and obtaining a final deposition rate change rule function f (t). In copper ion concentration prediction, the construction of a time series is critical to the accuracy of the prediction. Taking an actual electroplating process as an example, when the sampling interval delta t is set to 0.5 seconds, the original concentration data points {0s,0.6mol/L, 0.47s,0.58mol/L, 0.91s,0.55mol/L } are resampled to obtain an equidistant sequence {0s,0.6mol/L, 0.5s,0.59mol/L, 1.0s,0.55mol/L }, and the concentration change rate is-0.02 mol/(L.s) and-0.08 mol/(L.s) through differential calculation. The time window width tau of the long and short time memory network is set to 10 sampling points (namely 5 seconds), the input layer comprises concentration change rate values in windows corresponding to 10 neurons, the hidden layer is provided with 64 neurons, and the output layer predicts the concentration value at the next time by 1 neuron. Through 500 iterative training, the prediction error is reduced from the initial 0.15mol/L to 0.01mol/L. Based on Faraday's law, the calculation formula of the deposition thickness increment of copper ions is Δh= (M/nF.ρ). J.DELTA.t, wherein M is 63.5g/mol of molar mass of copper, n is valence number 2,F is Faraday constant 96485C/mol, ρ is copper density 8.9g/cm 3, and j is local current density. The theoretical deposition thickness increase in 1 second was 0.017 μm at a current density of 50mA/cm 2. And performing fast Fourier transform on the accumulated thickness sequence, wherein the sampling point number is 1024, and the main fluctuation period in the frequency spectrum is 20 seconds and 60 seconds, and the corresponding amplitude values are 0.05 mu m and 0.03 mu m respectively. This indicates that there are two main components of short cycle fluctuations and long cycle drifts in the deposition process. The deposition rate function v (t) =0.017+0.005sin (0.314 t) +0.003sin (0.105 t) μm/s reconstructed by inverse transformation. and 4-order polynomials are adopted to carry out least square fitting, the initial coefficient value {0.017, -0.001,0.0002, -0.00001,0.0000002} is subjected to least square iterative optimization to obtain the correction coefficient {0.0172, -0.00095,0.00018, -0.000012,0.00000018}, and the root mean square error of the fitted curve and the measured data is 0.002 mu m/s. When the polynomial order is increased to 5 th order, the error is reduced only to 0.0019 μm/s, the improvement is not significant, and therefore it is determined to use a 4 th order polynomial as a final deposition rate variation law function. The function reveals the nonlinear change characteristic of the deposition rate along with time, namely the rapid decline of the deposition rate in the initial stage (0-20 seconds), the stabilization of the deposition rate in the middle stage (20-60 seconds) and the slow rising trend in the later stage (> 60 seconds), which is closely related to the evolution of the electrode surface state and the dynamic change of the ion concentration distribution. By the multi-time scale analysis method, the rapid fluctuation characteristic is captured, and the long-term change trend is reflected.
Step S105, calculating the deposition thickness of each Chang Nongdu group at the future time according to the deposition rate change rule diagram, obtaining the copper layer thickness distribution of the abnormal region, comparing the copper layer thickness distribution of the abnormal region with the thickness distribution of the normal region, and calculating the thickness deviation value of the abnormal region.
The method comprises the steps of constructing a fourth-order Dragon lattice tower integral equation according to a deposition rate change rule function, obtaining an abnormal concentration group accumulated deposition thickness distribution function through a weighted average method, carrying out interpolation fitting on a copper layer thickness of a normal region through a radial basis function, obtaining a reference thickness distribution function through three-dimensional grid division, carrying out difference calculation on the copper layer thickness distribution function of the abnormal region and the reference thickness distribution function to obtain a thickness deviation function, carrying out probability density estimation on the thickness deviation function through a Gaussian kernel function, and obtaining a deviation distribution map according to kernel function weight calculation.
Specifically, a fourth-order Dragon-Greek tower integral equation is constructed according to a deposition rate change rule function v (t), a time step H is set to calculate a state quantity k 1、k2、k3、k4, and a weighted average method is adopted to calculate the accumulated deposition thickness H (x, y, z, t) of each abnormal concentration group. And obtaining a normal region discrete point coordinate set P (x, y, z) through three-dimensional grid division, and performing interpolation fitting on the normal region copper layer thickness H 0 (x, y, z) by adopting a radial basis function to generate a continuous reference thickness distribution function B (x, y, z). And calculating the mean value mu and the standard deviation sigma according to the standard thickness distribution function B (x, y, z) of the normal region, and establishing the thickness fluctuation range [ mu-3 sigma, mu+3 sigma ] of the normal region by adopting the three-time standard deviation principle as a standard. And calculating the difference value between the copper layer thickness distribution H (x, y, z, t) of the abnormal region and the reference thickness distribution function B (x, y, z) to generate a thickness deviation function D (x, y, z, t) =H (x, y, z, t) -B (x, y, z). And carrying out probability density estimation on the thickness deviation value by adopting a Gaussian kernel function, selecting a bandwidth parameter w to calculate kernel function weight, and generating a deviation value probability density function P (D). Dividing the deviation value interval according to the probability density function P (D), counting the spatial distribution characteristics of the deviation points in each interval, and generating a deviation distribution map G (x, y, z) by adopting a contour line method. In the process of calculating the deposition thickness of the copper layer, a fourth-order Dragon-Gregory tower method is adopted to carry out numerical integration on the deposition rate function. Taking the practical case as an example, given a deposition rate function v (t) =0.015+0.002sin (0.1 t) μm/s and a time step h=0.1 s, the state quantity calculation formula is k 1=h·v(t),k2=h·v(t+h/2),k3=h·v(t+h/2),k4 =h·v (t+h), and the final integration result h= (k 1+2k2+2k3+k4)/6. The cumulative thickness reached 0.152 μm at time t=10s. The normal region reference thickness distribution is interpolated by a multi-quadratic radial basis function in the form of B (r) = (1+ (epsilonr) 2)/(1/2), where r is the spatial distance and epsilons is the shape parameter taken to be 0.1. For typical normal area measurement data, the mean μ=0.15 μm, standard deviation σ=0.005 μm was calculated, and the reference range [0.135 μm,0.165 μm ] was established. In the thickness deviation calculation, the abnormal region appears a local thickness abrupt change phenomenon. Taking a certain abnormal concentration group as an example, the thickness of the center position reaches 0.18 mu m, and the thickness is different from the standard value by 0.03 mu m and is far beyond the normal fluctuation range. Probability density estimation was performed by gaussian kernel K (x) =exp (-x 2/2 h 2), with kernel bandwidth parameter h determined to be 0.0075 μm using SI LVERMAN criterion. The probability density distribution of the deviation value shows obvious double-peak characteristics, the main peak value is positioned at 0.02 mu m, the secondary peak value is positioned at-0.015 mu m, and the abnormal region has the thickening and thinning phenomena at the same time. On the spatial distribution, the thickening area is mainly concentrated at the edge of the electric field distortion area, the maximum deviation reaches 0.035 mu m, and the thinning area is arranged at the center of fluid disturbance, and the maximum deviation is-0.025 mu m. And drawing a deviation distribution map by a contour line method, and setting the contour line spacing to be 0.005 mu m. The thickness-increasing area is in elliptic distribution with long axis of 5mm and short axis of 3mm, and the thickness-decreasing area is in irregular circular shape with diameter of 4mm. There is a transition zone between the two regions, approximately 1mm wide, with a gradual transition from-0.01 μm to 0.01 μm in the deviation value. The z-direction profile shows that the deviation values decay exponentially with depth, tending to stabilize at 50 μm from the surface, reaching the reference thickness level. the three-dimensional distribution characteristics reflect the superposition result of the electric field effect and the fluid effect, and provide a quantitative basis for the optimization of the subsequent process.
Step S106, aiming at the thickness deviation value, adopting a deviation compensation algorithm, and adjusting electroplating process parameters of the abnormal region according to the magnitude and sign of the deviation value, wherein the electroplating process parameters comprise current density and electrolyte flow rate.
The method comprises the steps of establishing a proportional integral compensation equation according to a thickness deviation value, obtaining a current density compensation coefficient and a flow velocity compensation coefficient through deviation integral operation of a target value and an actual value, establishing a piecewise linear feedback function aiming at the current density compensation coefficient, obtaining corrected current density distribution through linear compensation if the absolute value of the compensation coefficient is smaller than a preset threshold value, transforming the flow velocity compensation coefficient through an exponential nonlinear mapping function, obtaining corrected flow velocity distribution through operation of the preset mapping coefficient and an original flow velocity value, constructing a compensation parameter evaluation function comprising a thickness deviation square term, a current density gradient term and a flow velocity gradient term, and obtaining the optimized compensation coefficient through optimization operation of a deep neural network on compensation parameters.
Specifically, a proportional integral compensation equation is established according to the thickness deviation value deltah (x, y, z), wherein u (t) =kp.e (t) +ki e (t) dt, u (t) is the output of the controller, kp and Ki are compensation coefficients, e (t) is the deviation between a target value and an actual value, and t is time, and the current density compensation coefficient alpha (x, y, z) and the flow velocity compensation coefficient beta (x, y, z) are calculated. A piecewise linear feedback function is established for the current density compensation coefficient α (x, y, z), with linear compensation if |α| is less than a threshold α 0, and saturation compensation if |α| is greater than α 0, to generate a corrected current density profile j' (x, y, z). Converting the flow rate compensation coefficient beta (x, y, z) by adopting an exponential nonlinear mapping function, wherein v '=v·exp (λβ), wherein λ is the mapping coefficient, and calculating to obtain a corrected flow rate distribution v' (x, y, z). Constructing a compensation parameter evaluation function:
Wherein J (alpha, beta) is an optimized target value of the compensation parameter, delta h is thickness deviation, J is current density distribution, v is flow velocity distribution, w 1、w2、w3 is a weight coefficient, and the optimized target value of the compensation parameter is calculated. And optimizing the compensation parameters by adopting a five-layer deep neural network, wherein the input layer is a deviation value and an original parameter, the hidden layer adopts a ReLU activation function, and the output layer is an optimized compensation coefficient. And recalculating the technological parameters according to the optimized compensation coefficient, verifying the compensation effect through numerical simulation, and generating a compensation parameter verification index R (x, y, z). In the electroplating process parameter compensation process, the deviation compensation is calculated by adopting a proportional integral algorithm. When a thickness deviation Δh=5 μm at a certain point is measured, a proportionality coefficient kp=0.2, an integral coefficient ki=0.05, and a corresponding compensation output u (t) =0.2×5+0.05×Σ5dtμm are set. The integral term accumulates over time and the compensation output gradually increases from 1 μm to 2.5 μm over 10 seconds. The current density compensation adopts piecewise function processing, a threshold value alpha 0 =0.3 is set, when the compensation coefficient alpha=0.2 is located in a linear interval, the corrected current density j '=j× (1+0.2) is obtained, when the compensation coefficient alpha=0.4 is obtained, the compensation current density exceeds the threshold value and enters a saturation region, and the saturation value j' =j×1.3 is obtained. This segmentation approach avoids overcompensation. The flow rate compensation adopts an exponential mapping relationship, and the mapping coefficient lambda=0.5. When β=0.4, the corrected flow velocity v' =v×exp (0.5×0.4) =1.22 v, exhibiting a nonlinear growth characteristic. The exponential relationship responds smoothly in the small compensation coefficient region and responds rapidly in the large compensation coefficient region. The weight coefficient in the evaluation function reflects the importance degree of different optimization targets, and the typical value is w 1=0.5,w2=0.3,w3 =0.2. Taking a certain abnormal region as an example, the evaluation value j=2.5 before compensation, wherein the thickness deviation term contributes 1.5, the current gradient term contributes 0.6, and the flow velocity gradient term contributes 0.4. The deep neural network adopts a five-layer structure, 8 neurons of the input layer receive the deviation value and the original parameter, three hidden layers respectively comprise 32, 16 and 8 neurons, and 2 neurons of the output layer output the optimized compensation coefficient. The hidden layer output is processed using a ReLU activation function f (x) =max (0, x), avoiding the gradient vanishing problem. During the verification phase, the effect of the compensation parameters is evaluated by the index R (x, y, z). The pre-compensation index R=0.85 of a certain test point is improved to R=0.92 after one round of optimization. Analytical optimization processes have found that current density compensation dominates the improvement in thickness uniformity, while flow rate compensation significantly reduces local disturbances. The synergy of the two compensation mechanisms is shown by the fact that the sum of squares of thickness deviations is reduced by 45%, the sum of squares of current gradients is reduced by 35%, and the sum of squares of flow gradients is reduced by 30% on multiple components of the evaluation function. This multi-objective optimization strategy not only improves deposition uniformity, but also maintains a smooth distribution of process parameters.
And step S107, setting and adjusting parameters of a power supply system and an electrolyte circulation system of the composite copper foil thickness control equipment according to the optimized electroplating process parameters, so that the composite copper foil thickness control equipment is produced according to the optimized electroplating process parameters.
Calculating the theoretical position coordinates of an anode plate in a target area according to a current density distribution function, measuring the actual position of the anode plate through a space distance sensor array, obtaining a position deviation compensation quantity through a proportional integral controller, obtaining a distance change rate matrix between the anode plate and a cathode plate through a laser displacement sensor array, calculating the inclination angle of the anode plate according to the distance change rate matrix, outputting an angle adjustment signal through an angle closed-loop controller, driving an anode plate position compensation mechanism to adjust if a current density uniformity index is smaller than a preset threshold value until the uniformity index is larger than the preset threshold value, calculating an electrolyte flow rate set value according to a flow field distribution function, measuring the actual flow rate value through a flow sensor, adjusting the opening of a valve through a pressure closed-loop controller, and adjusting the flow ratio of a liquid inlet pipeline to a liquid return pipeline according to the opening of the valve.
Specifically, the theoretical position coordinates (x 0,y0,z0) of the anode plate in the target area are calculated according to the optimized current density distribution function j (x, y, z), the actual position (x, y, z) of the anode plate is measured through the space distance sensor array, and the position deviation compensation quantity deltar is calculated through the proportional integral controller. The laser displacement sensor array is adopted to measure the distance change rate between the anode plate and the cathode plate in real time, the optimal inclination angle theta (x, y) of the anode plate is calculated according to the distance change rate matrix D (x, y, z), and an adjusting signal is output through the angle closed-loop controller. And calculating anode plate fine tuning parameters according to the current density uniformity index, and driving an anode plate position compensation mechanism to carry out fine adjustment if the uniformity index is smaller than a preset threshold value eta 0 until the uniformity index is larger than the threshold value eta 0. And calculating an electrolyte flow rate set value v 0 by adopting the optimized flow field distribution function v (x, y, z), measuring an actual flow rate v by a flow sensor, and regulating the rotation speed n of the circulating pump according to the deviation Deltav=v 0 -v. The electrolyte pipeline pressure p (t) is collected in real time, a third-order Butterworth low-pass filter is adopted to filter pressure signals, and the valve opening compensation quantity phi is calculated through a pressure closed-loop controller. And adjusting the flow ratio gamma of the liquid inlet pipeline and the liquid return pipeline according to the valve opening compensation quantity phi, monitoring the liquid level h (t) of the electroplating bath through a liquid level sensor, and establishing a liquid level stable control loop. In the parameter adjustment process of electroplating equipment, a multistage closed-loop feedback mechanism is adopted for anode plate position control. When the theoretically calculated optimal position coordinates of the anode plate are (150 mm,200mm,80 mm), the actual positions measured by the space distance sensor are (152 mm,198mm,82 mm), and the position deviation vector Deltar= (2, -2, 2) mm. And setting a proportional coefficient Kp=0.8 and an integral coefficient Ki=0.2 by adopting a proportional-integral controller, calculating to obtain a position compensation amount, and driving a servo motor to carry out fine adjustment. In the inclination angle adjustment of the anode plate, a laser displacement sensor array is arranged around the cathode plate, and the sampling frequency is 100Hz. When the measured distance change rate matrix shows that the right distance increase rate is 0.5mm/s and the left distance decrease rate is-0.4 mm/s, the anode plate has a rightward inclination tendency. And obtaining the optimal inclination angle theta=2.5 degrees through matrix operation, and outputting a corresponding stepping motor pulse signal by the angle controller. The current density uniformity index η is calculated by using the ratio of the standard deviation to the average value, and a preset threshold η 0 =0.05. the measured data shows that the uniformity index η=0.08 of a certain region is out of the threshold range. At this time, the precise position compensation mechanism is activated, fine adjustment is carried out by taking 0.1mm as a stepping amount, and the uniformity index is reduced to 0.045 after 3 iterations. In terms of flow field control, the optimized flow velocity distribution function gives a set value v 0 =0.15 m/s, and the flow sensor measures the actual flow velocity v=0.12 m/s. The rotation speed increment delta n=300 rpm of the circulating pump is calculated by the PI controller, and the frequency converter adjusts the output frequency accordingly. Meanwhile, the pressure sensor monitors pipeline pressure fluctuation and contains 50Hz power frequency interference. And processing the pressure signal by adopting a third-order Butterworth filter with the cut-off frequency of 10Hz, and filtering high-frequency noise to obtain stable pressure reading p=0.2 MPa. The liquid level control adopts a double-pipeline balance strategy, and the flow ratio gamma of the liquid inlet pipeline and the liquid return pipeline is realized by adjusting the opening of a valve. When the opening of the liquid inlet valve is 60%, the opening of the liquid return valve is adjusted to 55%, and the liquid level is kept at the target height h 0 =500 mm. If the liquid level fluctuates by +/-10 mm, the liquid level controller automatically adjusts the opening of the valve, and the compensation quantity phi is determined by the deviation value and the change rate. The multiparameter coupling control scheme realizes dynamic balance of technological parameters, and the control loops cooperate with each other to maintain stable operation of electroplating process.
The preferred embodiments of the invention disclosed above are intended only to assist in the explanation of the invention. The preferred embodiments are not intended to be exhaustive or to limit the invention to the precise form disclosed. Obviously, many modifications and variations are possible in light of the above teaching. The embodiments were chosen and described in order to best explain the principles of the invention and the practical application, to thereby enable others skilled in the art to best understand and utilize the invention. The invention is limited only by the claims and the full scope and equivalents thereof.