Electro-thermal co-verification method and system for heterogeneous process 3D stacked interconnection structure

Abstract

The application provides an electro-thermal cooperative verification method and system for a heterogeneous chemical process 3D stacked interconnection structure, relates to the technical field of integrated circuit manufacturing, and comprises the following steps: collecting physical layout information, calculating carrier migration trajectory and actual conductivity distribution, performing grid division and subdivision, calculating current density, joule heat loss and temperature distribution, setting temperature monitoring points, marking hot spot regions, dynamically adjusting functional module working parameters, and realizing temperature control. The application can improve electro-thermal verification precision, reduce the temperature of hot spot regions, prolong the service life of a chip, and optimize the performance of a stacked structure.

Description

Technical Field

[0001] This invention relates to the field of integrated circuit manufacturing technology, and in particular to an electrothermal co-verification method and system for heterogeneous process 3D stacked interconnect structures. Background Technology

[0002] With the continuous development of integrated circuit technology, in order to overcome the performance bottlenecks and chip area limitations faced by traditional planar processes, heterogeneous 3D stacked interconnect structures have become an important development direction for the semiconductor industry. This technology vertically stacks chips with different process technologies and uses interconnection methods such as microbumps and vias to achieve signal transmission between different functional units, effectively improving the chip's integration and performance.

[0003] In heterogeneous 3D stacked structures, each chip layer may employ different process nodes and materials, resulting in significant differences in their physical and electrical properties. This heterogeneity leads to complex thermoelectric coupling effects during chip operation, and the heat accumulation and transfer paths between layers become more complex. With increasing integration density and operating frequency, hotspots within the chip become increasingly prominent. Uneven temperature distribution not only affects chip performance and reliability but may also lead to premature failure.

[0004] Existing verification technologies for heterogeneous 3D stacked interconnect structures have several defects and shortcomings. They typically consider electrical and thermal characteristics separately, ignoring the coupling effect between the two. The lack of a co-verification mechanism for electrothermal processes makes it difficult for existing methods to accurately evaluate the performance and reliability of chips under actual operating conditions. Traditional verification methods use a uniform mesh partitioning strategy to model and analyze the chip structure, which is computationally inefficient when dealing with complex 3D stacked structures. For areas with high current density, uniform meshes cannot provide sufficient accuracy to capture local temperature changes and hotspot formation processes. In heterogeneous 3D stacked structures, the thermal coupling effect between chip layers is significant, and heat accumulation is difficult to dissipate in a timely manner. Without targeted thermal management strategies, this can easily lead to excessively high local temperatures, accelerating chip aging and even causing permanent damage. Summary of the Invention

[0005] This invention provides a method and system for electrothermal co-verification of heterogeneous 3D stacked interconnect structures, which can solve the problems in the prior art.

[0006] A first aspect of this invention provides a method for electrothermal co-verification of heterogeneous 3D stacked interconnect structures, comprising:

[0007] Collect physical layout information of the 3D stacked interconnect structure of the heterogeneous process to generate initial stress distribution data;

[0008] Based on the initial stress distribution data, the migration path of charge carriers in each chip cell is calculated to obtain the charge carrier migration trajectory matrix, and the actual conductivity distribution of each chip cell is calculated.

[0009] The heterogeneous 3D stacked interconnect structure is divided into multiple grid cells, and a grid subdivision operation is performed in the region where the current density is greater than a preset density threshold to generate grid cell data.

[0010] Based on the grid cell data and the actual conductivity distribution, the current density and Joule heat loss of each grid cell are calculated to obtain the initial heat distribution data, and the initial temperature value of each grid cell is calculated.

[0011] Temperature monitoring points are set in the grid cells to collect temperature data in real time and compare it with the initial temperature value. Grid cells with a temperature change rate exceeding the preset temperature change threshold and their corresponding adjacent cells are marked as hotspot areas.

[0012] Calculate the temperature gradient in the hot spot area, combine it with the power consumption sensitivity parameters of the functional modules to generate power adjustment commands, and dynamically adjust the operating parameters of the functional modules.

[0013] Based on the adjusted operating parameters, the temperature distribution of each grid cell is recalculated, and temperature monitoring and power adjustment are repeated until the temperature at all monitoring points is lower than the preset temperature threshold, and the verification results are output.

[0014] In one optional embodiment, based on the initial stress distribution data, the migration paths of charge carriers within each chip cell are calculated to obtain the charge carrier migration trajectory matrix. The calculation of the actual conductivity distribution of each chip cell includes:

[0015] The effective mass tensor and deformation potential tensor of the charge carriers are calculated based on the initial stress distribution data. A modified drift-diffusion equation is constructed, and the stress-induced potential is calculated. The carrier motion equation is established based on the modified drift-diffusion equation and the stress-induced potential, and spatiotemporal discretization parameters are generated.

[0016] A scattering effect model is constructed based on the spatiotemporal discretization parameters and the carrier motion equation; the migration path of carriers in each chip cell is determined based on the scattering effect model, a carrier migration trajectory matrix is ​​generated, and the local conductivity parameters of each chip cell, including electron mobility and hole mobility, are calculated.

[0017] The local conductivity parameters are interpolated using a spatial interpolation algorithm, and anisotropy correction is performed to obtain the actual conductivity distribution of each chip cell layer.

[0018] When the calculation error of the actual conductivity distribution is greater than the preset error threshold, the time step parameter and grid size parameter are adjusted to generate spatiotemporal discretization update parameters. The scattering effect model is then reconstructed based on the spatiotemporal discretization update parameters and the carrier motion equation until the calculation error of the actual conductivity distribution is less than the preset error threshold, and the actual conductivity distribution is output.

[0019] In one optional embodiment, constructing a scattering effect model based on spatiotemporal discretization parameters and carrier motion equations includes:

[0020] Phonon scattering potential energy and impurity scattering potential field are obtained based on spatiotemporal discretization parameters and carrier motion equations; phonon occupancy number is calculated based on the phonon scattering potential energy, and phonon scattering rate and impurity scattering rate are calculated based on the phonon occupancy number and the impurity scattering potential field, a scattering intensity evaluation matrix is ​​constructed, and a scattering intensity threshold is set.

[0021] The scattering intensity in the scattering intensity evaluation matrix is ​​compared with the scattering intensity threshold to determine the dominant scattering mechanism. The mean free path and characteristic scattering length are calculated based on the dominant scattering mechanism to generate grid scale parameters.

[0022] An electronic wave function is constructed based on the grid scale parameters, and the corresponding eigenstates are calculated. A scattering path sequence is constructed based on the eigenstates, the scattering probability of each scattering path is calculated, the scattering probability of the scattering path is compared with a preset probability threshold, effective scattering paths are selected, the scattering time of the effective scattering paths is calculated, and a scattering time sequence is generated.

[0023] Calculate the relative error between two adjacent scattering time series. When the relative error is greater than a preset convergence threshold, update the grid scale parameter and return to the construction step. Output the final scattering effect model when the relative error is less than the preset convergence threshold.

[0024] In one optional embodiment, based on the grid cell data and the actual conductivity distribution, the current density and Joule heat loss of each grid cell are calculated to obtain initial heat distribution data. The initial temperature value of each grid cell is then calculated, including:

[0025] The potential gradient is calculated based on the grid cell data and the actual conductivity distribution, generating the grid cell electric field strength. The grid cell current density is calculated based on the grid cell electric field strength, and the Joule heat power density is calculated based on the grid cell current density and the grid cell electric field strength, generating grid cell heat source distribution data.

[0026] A heat conduction equation is constructed based on the heat source distribution data of the grid cells. A heat conduction coefficient matrix is ​​generated according to the specific heat capacity and density of the material, and the heat flux distribution of the grid cells is calculated.

[0027] The heat diffusion rate is calculated based on the heat flux distribution of the grid cells, and the grid size parameters and temperature field iteration parameters are determined based on the heat diffusion rate.

[0028] A Hilbert-Huang transform model is constructed for the heat flux distribution of the grid cells to obtain the frequency-energy density spectrum. An adaptive filter bank is constructed, and the heat flux is divided into a first frequency band and a second frequency band through singular value decomposition. The temperature field is solved separately to obtain the first temperature distribution and the second temperature distribution.

[0029] A temperature field reconstruction weight function is constructed based on the first and second temperature distributions, and the initial temperature value of the grid cell is calculated.

[0030] When the ratio of the amplitude of the first temperature distribution to the amplitude of the second temperature distribution exceeds a preset ratio threshold, the grid size parameter is reduced by a preset reduction value and the temperature field iteration parameter is updated, and the Fourier transform step is returned to be executed.

[0031] In one optional embodiment, a Hilbert-Huang transform model is constructed to obtain the frequency-energy density spectrum of the heat flux distribution of the grid cells. An adaptive filter bank is constructed, and the heat flux is divided into a first frequency band and a second frequency band through singular value decomposition. The temperature field is then solved separately to obtain the first temperature distribution and the second temperature distribution, which include:

[0032] A Hilbert-Huang transform model is constructed for the heat flux distribution of grid cells, the eigenmode functions of the heat flux are extracted, the instantaneous frequency characteristics are calculated, and the frequency-energy density spectrum is generated.

[0033] An adaptive filter bank is constructed using the frequency-energy density spectrum, and the center frequency of the adaptive filter bank is updated according to the energy density distribution.

[0034] The adaptive filter bank is applied to the heat flux distribution to obtain multiple inherent frequency band components. Singular value decomposition is performed on the multiple inherent frequency band components to extract the frequency band feature vector.

[0035] Calculate the singular values ​​of the frequency band feature vector, and divide the singular values ​​into the first frequency band heat flux and the second frequency band heat flux according to the preset dividing boundary value;

[0036] The temperature field equations for the heat flux in the first frequency band and the heat flux in the second frequency band are constructed using a variational multi-scale method, respectively.

[0037] Calculate the temperature gradient value of the temperature field grid cell. When the temperature gradient value exceeds the preset gradient threshold, divide the corresponding grid cell into four sub-grid cells and update the grid cell.

[0038] Solve the temperature field equations of the updated mesh cells to obtain the first and second temperature distributions.

[0039] In one optional embodiment, the temperature gradient of the hotspot region is calculated, and a power adjustment command is generated in combination with the power consumption sensitivity parameter of the functional module. The operating parameters of the functional module are dynamically adjusted, including:

[0040] Extract temperature data of each marked grid cell within the hotspot area, calculate the temperature gradient vector between adjacent grid cells, calculate the magnitude of the temperature gradient vector, and label the temperature gradient distribution of the hotspot area;

[0041] Acquire real-time power consumption data and operating parameter data of each functional module, and establish a correlation model between power consumption and operating parameters; based on the correlation model, calculate the power consumption sensitivity parameters of each functional module, including the temperature sensitivity coefficient to power consumption and the power consumption sensitivity coefficient to operating parameters; construct a power adjustment objective function based on the temperature gradient vector of the hot spot area and the power consumption sensitivity parameters, including a temperature gradient term and a power consumption deviation term.

[0042] The gradient descent method is used to optimize and solve the power regulation objective function to obtain the target operating parameters of the functional module, and a power regulation command is generated based on the target operating parameters.

[0043] The system monitors the temperature distribution in the adjusted hotspot area in real time, calculates the temperature control error, and generates a comprehensive control quantity based on the integral, proportional, and derivative terms of the temperature control error. The system then corrects the target operating parameters based on the comprehensive control quantity and updates the power adjustment command.

[0044] A second aspect of the present invention provides an electrothermal co-verification system for heterogeneous 3D stacked interconnect structures, comprising:

[0045] The first unit is used to collect physical layout information of the heterogeneous process 3D stacked interconnect structure and generate initial stress distribution data;

[0046] The second unit is used to calculate the migration path of charge carriers in each chip cell based on the initial stress distribution data, obtain the charge carrier migration trajectory matrix, and calculate the actual conductivity distribution of each chip cell.

[0047] The third unit is used to divide the heterogeneous process 3D stacked interconnect structure into multiple grid units, and to perform grid subdivision operation in areas where the current density is greater than a preset density threshold to generate grid unit data.

[0048] The fourth unit is used to calculate the current density and Joule heat loss of each grid cell based on the grid cell data and the actual conductivity distribution, obtain the initial heat distribution data, and calculate the initial temperature value of each grid cell.

[0049] The fifth unit is used to set temperature monitoring points in the grid cells to collect temperature data in real time and compare it with the initial temperature value. Grid cells with a temperature change rate exceeding a preset temperature change threshold and their corresponding adjacent cells are marked as hotspot areas.

[0050] The sixth unit is used to calculate the temperature gradient of the hot spot area, and generate power adjustment instructions by combining the power consumption sensitivity parameters of the functional modules, so as to dynamically adjust the working parameters of the functional modules.

[0051] The seventh unit is used to recalculate the temperature distribution of each grid cell based on the adjusted operating parameters, repeat temperature monitoring and power adjustment until the temperature of all monitoring points is lower than the preset temperature threshold, and output the verification results.

[0052] A third aspect of the present invention provides an electronic device, comprising:

[0053] processor;

[0054] Memory used to store processor-executable instructions;

[0055] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0056] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0057] In this embodiment of the invention, the electrothermal co-verification method can accurately simulate the electrothermal coupling effect in the heterogeneous 3D stacked interconnect structure. By analyzing the initial stress distribution and carrier migration trajectory, it calculates a conductivity distribution that better reflects the actual working state, thus improving the accuracy of electrothermal analysis. A smart grid partitioning strategy is used to subdivide areas with high current density, optimizing the allocation of computing resources while ensuring computational accuracy, significantly improving verification efficiency, especially suitable for verification scenarios of complex multilayer heterogeneous integrated circuits. A closed-loop control mechanism for real-time temperature monitoring and dynamic power adjustment is implemented, enabling timely identification of hotspot areas and adaptive adjustment of operating parameters based on power consumption sensitivity parameters, effectively preventing local overheating, improving the reliability and lifespan of the heterogeneous 3D stacked structure, and providing reliable technical support for the design optimization of heterogeneous integrated circuits. Attached Figure Description

[0058] Figure 1 This is a flowchart illustrating the electrothermal co-verification method for a heterogeneous 3D stacked interconnect structure according to an embodiment of the present invention.

[0059] Figure 2A comparison chart of scattering path selection results under different electric field intensities;

[0060] Figure 3 This is a three-dimensional distribution diagram of the total temperature field superimposed from the two frequency bands. Detailed Implementation

[0061] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0062] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0063] Figure 1 This is a flowchart illustrating the electrothermal co-verification method for heterogeneous 3D stacked interconnect structures according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0064] Collect physical layout information of the 3D stacked interconnect structure of the heterogeneous process to generate initial stress distribution data;

[0065] Based on the initial stress distribution data, the migration path of charge carriers in each chip cell is calculated to obtain the charge carrier migration trajectory matrix, and the actual conductivity distribution of each chip cell is calculated.

[0066] The heterogeneous 3D stacked interconnect structure is divided into multiple grid cells, and a grid subdivision operation is performed in the region where the current density is greater than a preset density threshold to generate grid cell data.

[0067] Based on the grid cell data and the actual conductivity distribution, the current density and Joule heat loss of each grid cell are calculated to obtain the initial heat distribution data, and the initial temperature value of each grid cell is calculated.

[0068] Temperature monitoring points are set in the grid cells to collect temperature data in real time and compare it with the initial temperature value. Grid cells with a temperature change rate exceeding the preset temperature change threshold and their corresponding adjacent cells are marked as hotspot areas.

[0069] Calculate the temperature gradient in the hot spot area, combine it with the power consumption sensitivity parameters of the functional modules to generate power adjustment commands, and dynamically adjust the operating parameters of the functional modules.

[0070] Based on the adjusted operating parameters, the temperature distribution of each grid cell is recalculated, and temperature monitoring and power adjustment are repeated until the temperature at all monitoring points is lower than the preset temperature threshold, and the verification results are output.

[0071] In one optional implementation, based on the initial stress distribution data, the migration paths of charge carriers within each chip cell are calculated to obtain the charge carrier migration trajectory matrix. The calculation of the actual conductivity distribution of each chip cell includes:

[0072] The effective mass tensor and deformation potential tensor of the charge carriers are calculated based on the initial stress distribution data. A modified drift-diffusion equation is constructed, and the stress-induced potential is calculated. The carrier motion equation is established based on the modified drift-diffusion equation and the stress-induced potential, and spatiotemporal discretization parameters are generated.

[0073] A scattering effect model is constructed based on the spatiotemporal discretization parameters and the carrier motion equation; the migration path of carriers in each chip cell is determined based on the scattering effect model, a carrier migration trajectory matrix is ​​generated, and the local conductivity parameters of each chip cell, including electron mobility and hole mobility, are calculated.

[0074] The local conductivity parameters are interpolated using a spatial interpolation algorithm, and anisotropy correction is performed to obtain the actual conductivity distribution of each chip cell layer.

[0075] When the calculation error of the actual conductivity distribution is greater than the preset error threshold, the time step parameter and grid size parameter are adjusted to generate spatiotemporal discretization update parameters. The scattering effect model is then reconstructed based on the spatiotemporal discretization update parameters and the carrier motion equation until the calculation error of the actual conductivity distribution is less than the preset error threshold, and the actual conductivity distribution is output.

[0076] In one specific implementation, the effective mass tensor and deformation potential tensor of charge carriers are calculated based on given initial stress distribution data. For silicon-based chip cells, the effective mass tensor of electrons is determined by considering the band structure changes under stress, with typical values ​​of 0.19m0 for the transverse mass and 0.916m0 for the longitudinal mass (m0 being the free electron mass). The deformation potential tensor is calculated using the strain-energy level relationship; for the conduction band bottom, the deformation potential components include a volumetric deformation potential of approximately 9.5 eV and a shear deformation potential of approximately 8.7 eV. These parameters are applied to modify the classical drift-diffusion equation, introducing a stress effect term to obtain the modified drift-diffusion equation. The modified term considers the anisotropic change in the effective mass of charge carriers under stress and is expressed through a stress-induced potential.

[0077] The stress-induced potential is calculated centered on each heterogeneous interface. For interlayer interconnect structures, its value is obtained by the concatenation of the local stress tensor and the deformation potential tensor. In a specific implementation, for a region with an interface stress of 200 MPa, the calculated electron stress-induced potential is approximately 35 meV, and the hole stress-induced potential is approximately 48 meV.

[0078] Based on the modified drift-diffusion equation and stress-induced potential, a carrier motion equation is established. This equation describes the dynamic behavior of carriers under stress, including drift, diffusion, and stress-induced terms. When generating the spatiotemporal discretization parameters, the time step is set to the picosecond level (typically 0.5 ps), the spatial grid size is set to the nanometer level (typically 5 nm), and an adaptive meshing technique is used to refine the mesh in regions with large stress gradients.

[0079] A scattering effect model is constructed based on spatiotemporal discretization parameters and carrier motion equations. This model comprehensively considers multiple scattering mechanisms, including phonon scattering, Coulomb scattering, and surface scattering, and combines the Monte Carlo method to simulate the scattering behavior of carriers in a stress field. For the silicon-silicon oxide interface, the surface scattering rate parameter is set to 2.3 × 10⁵ cm / s; for the metal-semiconductor contact interface, the interface scattering enhancement factor is set to 1.7.

[0080] The migration paths of charge carriers within each chip cell are determined based on a scattering effect model, generating a carrier migration trajectory matrix. For a specific chip cell (e.g., a 5μm×5μm×2μm silicon substrate region), the migration trajectories of 1000 sample charge carriers are tracked, and their position coordinates and velocity vectors in three-dimensional space are recorded to form a migration trajectory matrix. During the calculation, the effects of accelerating electric field, scattering events, and boundary conditions are considered. The migration trajectory data is used to calculate local conductivity parameters.

[0081] Based on the carrier migration trajectory matrix, the local conductivity parameters of each chip cell are calculated. Electron mobility is determined by the ratio of the average drift velocity of the migration trajectory to the electric field strength, and hole mobility is calculated similarly. In regions with large stress gradients, such as around interconnects, electron mobility decreases by approximately 15%, typically from 1350 cm² / V·s to 1150 cm² / V·s; hole mobility decreases by approximately 22%, from 450 cm² / V·s to 350 cm² / V·s.

[0082] Spatial interpolation algorithms were used to interpolate local conductivity parameters, including trilinear interpolation and radial basis function interpolation. Trilinear interpolation is suitable for regions with gentle stress gradients, while radial basis function interpolation is suitable for regions with steep stress gradients. After interpolation, anisotropy correction was performed, taking into account the angle between the lattice direction and the stress principal axis direction to correct the anisotropic properties of the conductivity tensor. For stress in the

[110] crystal direction, the principal axis ratio of the conductivity tensor is 1:0.78:0.65.

[0083] When the calculation error of the actual conductivity distribution exceeds a preset error threshold (usually set to 3%), the time step and mesh size parameters are adjusted. In the high stress gradient region, the time step is reduced to 50% of its original size (e.g., from 0.5 ps to 0.25 ps), and the mesh size is reduced to 40% of its original size (e.g., from 5 nm to 2 nm). After generating the spatiotemporal discretization update parameters, the scattering effect model is reconstructed until the calculation error of the actual conductivity distribution is less than the preset error threshold.

[0084] In a practical application, the conductivity distribution of a 3D stacked structure manufactured using a 10nm process was verified. The initial mesh size was set to 5nm, and the time step was 0.5ps. After three iterations of optimization, the calculation error was reduced to 2.4%, below the preset threshold of 3%. The final conductivity distribution data showed that in the interlayer interconnect structure, the conductivity in the region most affected by stress decreased by 17.5% compared to the nominal value. This is of significant importance for circuit performance analysis.

[0085] This verification method can also be extended to different process nodes and material systems. For example, for GaN-based chips, by adjusting the corresponding effective mass tensor (typical values ​​of lateral mass 0.20m0 and longitudinal mass 0.98m0) and deformation potential tensor (conduction band bottom deformation potential of about 8.3eV), the conductivity distribution of III-V semiconductor heterostructures can be verified.

[0086] The output actual conductivity distribution data can be used for subsequent electrothermal co-analysis, providing a scientific basis for the reliability assessment and optimization design of 3D stacked interconnect structures. The conductivity accuracy verified by this method is significantly better than that of traditional methods, and the computational efficiency is improved by about 40%, providing strong support for the design of 3D stacked structures with heterogeneous processes.

[0087] In one alternative implementation, constructing a scattering effect model based on spatiotemporal discretization parameters and carrier motion equations includes:

[0088] Phonon scattering potential energy and impurity scattering potential field are obtained based on spatiotemporal discretization parameters and carrier motion equations; phonon occupancy number is calculated based on the phonon scattering potential energy, and phonon scattering rate and impurity scattering rate are calculated based on the phonon occupancy number and the impurity scattering potential field, a scattering intensity evaluation matrix is ​​constructed, and a scattering intensity threshold is set.

[0089] The scattering intensity in the scattering intensity evaluation matrix is ​​compared with the scattering intensity threshold to determine the dominant scattering mechanism. The mean free path and characteristic scattering length are calculated based on the dominant scattering mechanism to generate grid scale parameters.

[0090] An electronic wave function is constructed based on the grid scale parameters, and the corresponding eigenstates are calculated. A scattering path sequence is constructed based on the eigenstates, the scattering probability of each scattering path is calculated, the scattering probability of the scattering path is compared with a preset probability threshold, effective scattering paths are selected, the scattering time of the effective scattering paths is calculated, and a scattering time sequence is generated.

[0091] Calculate the relative error between two adjacent scattering time series. When the relative error is greater than a preset convergence threshold, update the grid scale parameter and return to the construction step. Output the final scattering effect model when the relative error is less than the preset convergence threshold.

[0092] In one specific implementation, when constructing the scattering effect model, the system first obtains the phonon scattering potential energy and impurity scattering potential field based on the spatiotemporal discretization parameters and the carrier motion equation. Specifically, the system uses the lattice constant of the semiconductor material as the spatial discretization parameter; for silicon, the lattice constant can be set to 0.543 nm; the temporal discretization parameter can be set to 1 femtosecond. The carrier motion equation is based on the Boltzmann transport equation, considering the acceleration and scattering process of carriers under the action of an electric field. Based on this, the system calculates the phonon scattering potential energy; for silicon, the acoustic phonon scattering potential energy is approximately 9.5 eV, and the optical phonon scattering potential energy is approximately 11.0 eV. The impurity scattering potential field is calculated using the Coulomb potential; for n-type silicon with a doping concentration of 1 × 10¹⁸ / cm³, the average value of the impurity scattering potential field is approximately 0.3 eV.

[0093] The system calculates the phonon occupancy number based on the calculated phonon scattering potential energy using the Bose-Einstein distribution. At room temperature (300 K), for silicon, the occupancy number of the transverse optical phonon mode is approximately 0.18, and the occupancy number of the longitudinal optical phonon mode is approximately 0.21. Combining the phonon occupancy number and the impurity scattering potential field, the system calculates the phonon scattering rate and the impurity scattering rate. For n-type silicon with a doping concentration of 1 × 10¹⁸ / cm³, the phonon scattering rate is approximately 3.2 × 10¹³ / s, and the impurity scattering rate is approximately 1.5 × 10¹³ / s. The system then constructs a scattering intensity evaluation matrix, where the matrix elements represent the scattering intensity between different energy and momentum states, and sets the scattering intensity threshold to 2 × 10¹³ / s.

[0094] The system compares the scattering intensities in the scattering intensity assessment matrix with set thresholds to determine the dominant scattering mechanism. In moderately doped silicon materials at room temperature, phonon scattering is typically the dominant mechanism. Based on the dominant scattering mechanism, the system calculates the carrier mean free path and characteristic scattering length. Under the above conditions, the electron mean free path is approximately 15 nm, and the characteristic scattering length is approximately 22 nm. The system generates grid scale parameters based on these parameters, with the initial grid scale set to 1 / 10 of the mean free path, i.e., 1.5 nm.

[0095] Using the generated grid scale parameters, the system constructs an electronic wavefunction. The wavefunction adopts a Gaussian wave packet form, with an initial width of 5 nanometers and a central momentum corresponding to an energy of 0.2 electron volts. The system solves the Schrödinger equation to obtain the eigenstates of the wavefunction; typically, calculating the first 20 eigenstates is sufficient to represent the main quantum properties of the system. Based on the calculated eigenstates, the system constructs a scattering path sequence. The scattering path sequence contains transitions from the initial state to each possible final state, and each path records information such as the scattering type, initial and final states, and scattering angle.

[0096] The system calculates the scattering probability of each scattering path. For example, from an initial state with energy of 0.2 eV and wave vector direction

[100] to a final state with energy of 0.18 eV and wave vector direction

[110] , the probability of acoustic phonon scattering is approximately 0.025. The system compares the scattering probability of each scattering path with a preset probability threshold of 0.001 and selects effective scattering paths with probabilities greater than the threshold. For each effective scattering path, the system calculates the scattering time and generates a scattering time series. The typical phonon scattering time is approximately 0.03 picoseconds, and the impurity scattering time is approximately 0.06 picoseconds.

[0097] The system calculates the relative error between two adjacent scattering time series, defined as the difference between the average values ​​of the scattering time series between two iterations divided by the average value of the previous iteration. When the relative error exceeds a preset convergence threshold (usually set to 1%), the system updates the grid scale parameter based on the current scattering characteristics. For regions with strong scattering, the grid scale is reduced to 0.8 times its original size; for regions with weak scattering, the grid scale is increased to 1.2 times its original size. The system then returns to the step of constructing the electron wavefunction, recalculating the eigenstates and scattering paths.

[0098] The system repeats the above iterative process until the relative error is less than the preset convergence threshold, and finally outputs the scattering effect model. The model includes optimized grid parameters (typically in the range of 0.8-2.5 nanometers), scattering time series (typical values ​​in the range of 0.02-0.08 picoseconds), and the relative contributions of each scattering mechanism (e.g., phonon scattering accounts for 67%, and impurity scattering accounts for 33%).

[0099] In practical applications, the model parameters can be adjusted to suit devices with different operating temperatures and doping concentrations. For example, when the temperature rises to 400K, the phonon scattering rate increases by about 30%, and the grid scale parameter needs to be adjusted to 0.85 times its original value. When the doping concentration increases to 5×10¹⁸ / cm³, impurity scattering becomes the dominant mechanism, and the characteristic scattering length decreases to about 8 nanometers, requiring the grid scale parameter to be adjusted to 0.6 times its original value. This scattering effect model can be effectively applied to the performance prediction and optimization design of nanoscale semiconductor devices.

[0100] like Figure 2 As shown in the comparison chart of scattering path screening results under different electric field intensities, the performance of scattering path screening is comprehensively demonstrated within the range of electric field intensities from 0.5 kV / cm to 5.0 kV / cm. With increasing electric field intensity, the total number of scattering paths shows a significant upward trend, increasing from 342 to 985, indicating that the scattering process becomes more complex under high electric field conditions. The number of effective paths screened by this technical solution is significantly less than that of the traditional method. At an electric field intensity of 5.0 kV / cm, this technical solution screened 346 effective paths, while the traditional method screened 473, a reduction of approximately 27%. In terms of screening efficiency, the screening efficiency of this technical solution fluctuates between 25.44% and 38.04%, with the highest point occurring at 3.5 kV / cm, while the screening efficiency of the traditional method ranges from 41.81% to 49.68%, consistently higher than that of this technical solution. This indicates that this technical solution employs stricter screening criteria and can more effectively eliminate scattering paths that contribute little to the final result. The computational cost savings data further confirm this. Over the entire electric field strength range, this technical solution saves an average of approximately 28% of the computational cost compared to traditional methods, especially in the low electric field region (0.5 kV / cm), where the cost savings reach as high as 39.17%. This efficient path selection method not only improves computational efficiency but also enhances the accuracy of the final scattering effect model by more precisely focusing on physically significant scattering processes.

[0101] In one optional implementation, based on the grid cell data and the actual conductivity distribution, the current density and Joule heat loss of each grid cell are calculated to obtain initial heat distribution data. The initial temperature value of each grid cell is then calculated, including:

[0102] The potential gradient is calculated based on the grid cell data and the actual conductivity distribution, generating the grid cell electric field strength. The grid cell current density is calculated based on the grid cell electric field strength, and the Joule heat power density is calculated based on the grid cell current density and the grid cell electric field strength, generating grid cell heat source distribution data.

[0103] A heat conduction equation is constructed based on the heat source distribution data of the grid cells. A heat conduction coefficient matrix is ​​generated according to the specific heat capacity and density of the material, and the heat flux distribution of the grid cells is calculated.

[0104] The heat diffusion rate is calculated based on the heat flux distribution of the grid cells, and the grid size parameters and temperature field iteration parameters are determined based on the heat diffusion rate.

[0105] A Hilbert-Huang transform model is constructed for the heat flux distribution of the grid cells to obtain the frequency-energy density spectrum. An adaptive filter bank is constructed, and the heat flux is divided into a first frequency band and a second frequency band through singular value decomposition. The temperature field is solved separately to obtain the first temperature distribution and the second temperature distribution.

[0106] A temperature field reconstruction weight function is constructed based on the first and second temperature distributions, and the initial temperature value of the grid cell is calculated.

[0107] When the ratio of the amplitude of the first temperature distribution to the amplitude of the second temperature distribution exceeds a preset ratio threshold, the grid size parameter is reduced by a preset reduction value and the temperature field iteration parameter is updated, and the Fourier transform step is returned to be executed.

[0108] In one specific implementation, the potential gradient is calculated based on grid cell data and the actual conductivity distribution. In the implementation process, the stacked interconnect structure is discretized into grid cells with dimensions of 5nm × 5nm × 2nm, and the potential gradient is calculated for each grid cell using the finite difference method. The potential gradient calculation considers the ratio of the potential difference between grid cell boundaries to the grid cell size. Boundary conditions include a source potential of 1.0V, a drain potential of 0V, and other boundaries set as insulating boundary conditions. For grid cells at heterogeneous interfaces, considering conductivity discontinuities, harmonic averaging is used to process the interface conductivity.

[0109] The electric field strength of the grid cells is generated based on the calculated potential gradient. For each grid cell, the electric field strength vector is determined by the negative value of the potential gradient. Calculation results show that the electric field strength reaches the order of 10⁶ V / m around the TSV, while it drops to the order of 10⁴ V / m in regions far from the TSV. At interconnect vias and line intersections, the local electric field strength can reach 1.8 × 10⁶ V / m.

[0110] The current density of the grid cells is calculated based on the electric field strength of the grid cells. The current density is calculated using Ohm's law, which is the product of conductivity and electric field strength. Considering the anisotropic nature of conductivity under stress, tensor multiplication is used for the current density calculation. Results show that the typical current density is 5 × 10¹⁰ A / m² in the metal interconnect region, while it drops to the order of 10⁶ A / m² in the silicon substrate region.

[0111] The Joule thermal power density is calculated based on the grid cell current density and grid cell electric field strength. The Joule thermal power density is equal to the scalar product of the current density and the electric field strength. The calculation results show that in the high-stress region near the metal / semiconductor interface, the Joule thermal power density is approximately 1.2 × 10¹⁶ W / m³, while it drops to the order of 10¹⁴ W / m³ in regions far from the interface. This method generates grid cell heat source distribution data, which serves as input for subsequent thermal analysis.

[0112] A heat conduction equation is constructed based on the heat source distribution data of the grid cells. This equation describes the relationship between temperature and time and space, and includes heat source terms, thermal conduction terms, and heat capacity terms. A thermal conductivity matrix is ​​generated based on the specific heat capacity and density of the materials. For silicon substrates, the thermal conductivity is 148 W / (m·K), the specific heat capacity is 705 J / (kg·K), and the density is 2330 kg / m³. For copper interconnects, the thermal conductivity is 398 W / (m·K), the specific heat capacity is 385 J / (kg·K), and the density is 8960 kg / m³. For silicon dioxide dielectric layers, the thermal conductivity is 1.4 W / (m·K), the specific heat capacity is 730 J / (kg·K), and the density is 2200 kg / m³. At heterogeneous interfaces, the thermal conductivity is calculated using harmonic averaging to ensure heat flux continuity.

[0113] The heat flux distribution within the grid cells was calculated. The heat flux distribution was calculated using Fourier's law of heat conduction, which is the product of the heat conduction coefficient matrix and the temperature gradient. Under steady-state conditions, the temperature distribution was obtained by iteratively solving the heat conduction equation, and then the heat flux distribution was calculated. The results show that the heat flux density can reach 10⁵ W / m² in hotspot regions, while it drops to the order of 10³ W / m² in regions far from hotspots.

[0114] The heat diffusion rate was calculated based on the heat flux distribution of the grid cells. The heat diffusion rate was calculated using the thermal diffusivity coefficient, which is the thermal conductivity coefficient divided by the product of density and specific heat capacity. The calculation results show that the thermal diffusion rate of the silicon substrate is approximately 8.8 × 10⁻⁵ m² / s, the thermal diffusion rate of the copper interconnect is approximately 1.15 × 10⁻⁴ m² / s, and the thermal diffusion rate of the silicon dioxide dielectric layer is approximately 8.7 × 10⁻⁷ m² / s. Based on the heat diffusion rate, the grid size parameters and temperature field iteration parameters were determined. The grid size parameters should be less than the minimum thermal diffusion length, and the temperature field iteration parameters were set with a convergence threshold of 0.01 K and a maximum number of iterations of 1000.

[0115] A Hilbert-Huang transform model is constructed to obtain the frequency-energy density spectrum of the heat flux distribution in the grid cells. The Hilbert-Huang transform decomposes the heat flux distribution into different frequency components, yielding the frequency-energy density spectrum. An adaptive filter bank is constructed, including a low-pass filter and a high-pass filter, with the filter cutoff frequency set to 0.2 nm⁻¹. Singular value decomposition is used to divide the heat flux distribution into a first frequency band and a second frequency band, with the first band corresponding to low-frequency thermal diffusion processes and the second band corresponding to high-frequency thermal diffusion processes.

[0116] The temperature fields corresponding to the heat flux in the first and second frequency bands are solved separately to obtain the first and second temperature distributions. The first temperature distribution mainly reflects the global heat diffusion characteristics, with a gradual temperature change and a maximum temperature of approximately 353 K; the second temperature distribution mainly reflects the local heat diffusion characteristics, with a drastic temperature change and a temperature fluctuation range of approximately ±15 K.

[0117] A temperature field reconstruction weighting function was constructed based on the first and second temperature distributions. The weighting function was designed according to the heat flux spectrum characteristics, with a weight of 0.75 for the low-frequency portion and 0.25 for the high-frequency portion. The initial temperature values ​​of the grid cells were calculated by combining the first and second temperature distributions using the weighting function. The results show that the temperature rises by approximately 35 K around the TSV and by approximately 28 K at the line intersection, forming significant hotspots.

[0118] When the ratio of the amplitude of the first temperature distribution to the amplitude of the second temperature distribution exceeds a preset ratio threshold, the mesh size parameter is reduced by a preset reduction value, and the temperature field iteration parameters are updated. The preset ratio threshold is set to 3.0, and the preset reduction value is set to 0.5, meaning the mesh size is halved. For example, when the ratio is 3.5, the mesh size is reduced from 5nm×5nm×2nm to 2.5nm×2.5nm×1nm, the temperature field iteration parameters are updated to the convergence threshold of 0.005K, and the maximum number of iterations is 1500. Then, the process returns to the Hilbert-Huang transform step to resolve the temperature field.

[0119] In a practical application case, electrothermal co-validation was performed on a 3D stacked structure using a 14nm process. The initial mesh size was 5nm. After two mesh optimizations, the mesh size was reduced to 1.25nm in the hotspot region. The amplitude ratio of the first temperature distribution to the second temperature distribution decreased to 2.8, below the preset ratio threshold of 3.0. The final temperature distribution showed that the highest temperature in the stacked structure reached 375K, occurring in the interconnect cross-section region with a power density of 1.5×10¹⁶ W / m³. The maximum temperature gradient was 15K / μm, occurring at interfaces with significant differences in thermal conductivity.

[0120] This electrothermal co-verification method accurately captures the current and temperature distribution characteristics in heterogeneous 3D stacked interconnect structures, providing important basis for optimized design and reliability assessment. Compared with traditional methods, the calculation accuracy is improved by approximately 30%, especially demonstrating higher prediction accuracy at heterogeneous interfaces and high-stress regions.

[0121] In one optional implementation, a Hilbert-Huang transform model is constructed to obtain the frequency-energy density spectrum of the heat flux distribution of the grid cells. An adaptive filter bank is then constructed, and the heat flux is divided into a first frequency band and a second frequency band through singular value decomposition. The temperature fields are then solved separately to obtain the first temperature distribution and the second temperature distribution, which include:

[0122] A Hilbert-Huang transform model is constructed for the heat flux distribution of grid cells, the eigenmode functions of the heat flux are extracted, the instantaneous frequency characteristics are calculated, and the frequency-energy density spectrum is generated.

[0123] An adaptive filter bank is constructed using the frequency-energy density spectrum, and the center frequency of the adaptive filter bank is updated according to the energy density distribution.

[0124] The adaptive filter bank is applied to the heat flux distribution to obtain multiple inherent frequency band components. Singular value decomposition is performed on the multiple inherent frequency band components to extract the frequency band feature vector.

[0125] Calculate the singular values ​​of the frequency band feature vector, and divide the singular values ​​into the first frequency band heat flux and the second frequency band heat flux according to the preset dividing boundary value;

[0126] The temperature field equations for the heat flux in the first frequency band and the heat flux in the second frequency band are constructed using a variational multi-scale method, respectively.

[0127] Calculate the temperature gradient value of the temperature field grid cell. When the temperature gradient value exceeds the preset gradient threshold, divide the corresponding grid cell into four sub-grid cells and update the grid cell.

[0128] Solve the temperature field equations of the updated mesh cells to obtain the first and second temperature distributions.

[0129] This embodiment provides a method for processing the heat flux distribution of grid cells and calculating the temperature field. The method uses the Hilbert-Huang transform model to analyze the heat flux distribution, and divides the heat flux into different frequency bands through adaptive filtering and singular value decomposition, solving for the temperature field separately to obtain more accurate temperature distribution results.

[0130] In one specific implementation, when constructing a Hilbert-Huang transform model for the heat flux distribution of grid cells, the heat flux time series data is first decomposed into a series of intrinsic mode functions (IMFs). Specifically, the heat flux distribution data is set as q(t). Extrema are found, and cubic spline interpolation is used to construct upper and lower envelopes. The average value of the upper and lower envelopes is calculated as the mean function m(t). The mean function is subtracted from the original signal to obtain h(t) = q(t) - m(t). It is then checked whether h(t) satisfies the IMF conditions: the difference between the number of extrema and the number of zero-crossing points does not exceed 1, and the local mean is zero. If not, the above process is repeated until the conditions are met, yielding the first IMF1. IMF1 is subtracted from the original signal to obtain the residual signal. The residual signal is then treated as a new processing object, and the above steps are repeated until the residual signal becomes a monotonic function or a constant. In practical cases, for heat conduction problems, typically 5-8 IMFs can be extracted.

[0131] When calculating the instantaneous frequency, a Hilbert transform is applied to each intrinsic mode function (IMFi) to obtain the analytic signal z(t) = IMFi(t) + jH[IMFi(t)], where H[IMFi(t)] is the Hilbert transform of IMFi(t). The instantaneous frequency is calculated as the time derivative of the phase angle. Based on the calculated instantaneous frequency and the amplitude of the intrinsic mode functions, a frequency-energy density spectrum is generated. For example, for heat conduction problems, the low-frequency range (0-0.1Hz) typically corresponds to internal heat diffusion in the material, the mid-frequency range (0.1-1Hz) corresponds to boundary condition effects, and the high-frequency range (>1Hz) corresponds to noise or measurement errors.

[0132] When constructing an adaptive filter bank using the frequency-energy density spectrum, the most representative frequency point is selected as the center frequency of the filter based on the energy density distribution. In practice, the energy density spectrum can be divided into several equal parts, and the energy peak point in each interval can be used as the center frequency. For heat conduction problems, 3-5 filters are generally sufficient to capture the main characteristics. The filter bandwidth is adaptively adjusted according to the frequency resolution, typically set to half the difference between adjacent center frequencies. For example, if adjacent center frequencies are 0.2Hz and 0.4Hz, the bandwidth is set to 0.1Hz.

[0133] An adaptive filter bank is applied to the heat flux distribution, and singular value decomposition (SVD) is performed on the resulting multiple intrinsic frequency band components. A data matrix A is constructed, where each column represents the time series of a frequency band component. SVD is then performed on matrix A: A = USVT, where U and V are the left and right singular vector matrices, respectively, and S is the singular value diagonal matrix. The column vectors of the left singular vector matrix U are the frequency band feature vectors.

[0134] The singular values ​​of the frequency band eigenvectors are calculated, and the frequency bands are divided into two categories based on the magnitude of the singular values. A preset division boundary value is set to 80% of the total energy. The singular values ​​are arranged in descending order, and the eigenvectors with a cumulative contribution rate of 80% are assigned to the first frequency band heat flux, while the rest are assigned to the second frequency band heat flux. For example, if the first two singular values ​​are 8.5 and 6.2, and the third singular value is 2.1, with a sum of 16.8, then the first two singular values ​​account for 87.5% of the total energy, exceeding the preset boundary value of 80%. Therefore, the frequency band components corresponding to the first two eigenvectors are assigned to the first frequency band heat flux, and the remaining frequency band components are assigned to the second frequency band heat flux.

[0135] The temperature field equations for the heat flux in the first and second frequency bands were constructed using a variational multiscale method. The temperature field function was decomposed into coarse-scale and fine-scale components, which were then substituted into the heat conduction governing equations to obtain two coupled equations. These two equations were solved iteratively to obtain the temperature fields for the first and second frequency bands, respectively.

[0136] The temperature gradient values ​​of the temperature field grid cells are calculated. When the temperature gradient value exceeds a preset gradient threshold, the grid is refined. The preset gradient threshold is typically set to twice the average temperature gradient. For example, if the average temperature gradient is 5℃ / mm, the preset gradient threshold is 10℃ / mm. For grid cells exceeding the threshold, they are divided into four equal sub-grid cells, and the grid topology and related calculation parameters are updated.

[0137] The temperature field equations for the first and second frequency bands are solved separately on the updated mesh. The governing equations are discretized using the finite element method, and the stiffness matrix and load vector are constructed. The linear equation system is then solved iteratively. For example, a temperature field calculation with an initial mesh of 5000 elements may increase to 8000-10000 elements after adaptive mesh refinement, with approximately 60% of the refinement occurring in regions with large temperature gradients, such as material interfaces or near heat sources. The final first and second temperature distributions reflect the influence of heat flux in different frequency bands on the temperature field.

[0138] like Figure 3As shown, the total temperature field formed by the superposition of two frequency bands is illustrated, resulting in a more precise temperature distribution. The highest temperature reaches 143.1°C, which is 16.7°C higher than the traditional method, representing a 13.2% improvement in accuracy. The table data shows that this technical solution achieves a maximum temperature gradient capture capability of 128.5°C / m, a 39.2% improvement over the traditional method. Through adaptive meshing and frequency band decomposition, the total number of meshes is reduced from 6400 in the traditional method to 2145, a reduction of 66.5%. The total computation time is reduced from 68.7 seconds to 23.4 seconds, a speedup of 65.9%. Overall, this technical solution, through precise frequency band division and adaptive mesh refinement, significantly reduces computational resource requirements while substantially improving the accuracy of temperature field calculations, especially in critical regions with large temperature gradients.

[0139] In one optional implementation, the temperature gradient of the hotspot region is calculated, and a power adjustment command is generated based on the power consumption sensitivity parameters of the functional module. The operating parameters of the functional module are dynamically adjusted, including:

[0140] Extract temperature data of each marked grid cell within the hotspot area, calculate the temperature gradient vector between adjacent grid cells, calculate the magnitude of the temperature gradient vector, and label the temperature gradient distribution of the hotspot area;

[0141] Acquire real-time power consumption data and operating parameter data of each functional module, and establish a correlation model between power consumption and operating parameters; based on the correlation model, calculate the power consumption sensitivity parameters of each functional module, including the temperature sensitivity coefficient to power consumption and the power consumption sensitivity coefficient to operating parameters; construct a power adjustment objective function based on the temperature gradient vector of the hot spot area and the power consumption sensitivity parameters, including a temperature gradient term and a power consumption deviation term.

[0142] The gradient descent method is used to optimize and solve the power regulation objective function to obtain the target operating parameters of the functional module, and a power regulation command is generated based on the target operating parameters.

[0143] The system monitors the temperature distribution in the adjusted hotspot area in real time, calculates the temperature control error, and generates a comprehensive control quantity based on the integral, proportional, and derivative terms of the temperature control error. The system then corrects the target operating parameters based on the comprehensive control quantity and updates the power adjustment command.

[0144] In one specific embodiment of the present invention, the method for calculating the temperature gradient of the hot spot area and generating a power adjustment command in combination with the power consumption sensitivity parameters of the functional module is described in detail below.

[0145] To extract temperature data from each marked grid cell within the hotspot region, the system first monitors the chip surface temperature in real time using a temperature sensor array. The temperature sensors are distributed at a density of four per square millimeter, with a sampling frequency of 10 Hz. The temperature data is then converted from analog to digital and stored in the system's temperature data buffer. The system divides the entire chip area into 0.5 mm × 0.5 mm grid cells, forming a 100 × 100 grid matrix. The system marks grid cells with temperatures exceeding 80°C within three consecutive sampling periods as part of the hotspot region. For example, in a single detection, the system identifies multiple grid cells within the region at coordinates (25, 30) to (35, 40) whose temperatures exceed the threshold, thus marking this region as a hotspot.

[0146] When calculating the temperature gradient vector between adjacent grid cells, the system calculates the temperature difference between each marked grid cell within the hotspot region and its four adjacent grid cells (upper, lower, left, and right), and divides this difference by the grid spacing to obtain the temperature gradient components in the x and y directions. For example, for the grid cell at coordinates (26, 31) with a temperature of 85℃, the grid cell to its right has a temperature of 82℃, and the grid cell below it has a temperature of 87℃. The gradient component in the x direction is -6℃ / mm, and the gradient component in the y direction is 4℃ / mm. These two components are combined to form the temperature gradient vector for that grid cell. The system calculates the magnitude of this vector, which is the temperature gradient at that grid cell; in this example, it is 7.21℃ / mm. Repeating this process, the system completes the temperature gradient distribution calibration for the entire hotspot region and stores it as a heatmap, with areas having a gradient value greater than 10℃ / mm marked as high gradient areas.

[0147] When acquiring real-time power consumption data for each functional module, the system uses a built-in power consumption monitoring circuit to monitor the power consumption of each module in the chip in real time. For example, for the graphics processing module, the system records its current power consumption as 3.2W; for the memory controller module, the current power consumption is 1.8W. Simultaneously, the system records the operating parameters of each module, including the clock frequency of the graphics processing module as 1.2GHz and the voltage as 0.95V; and the clock frequency of the memory controller as 800MHz and the voltage as 0.85V, etc.

[0148] The system establishes a correlation model between power consumption and operating parameters through regression analysis. Taking the graphics processing module as an example, the system collects 150 power consumption data points under different operating parameters over the past 30 minutes, including power consumption values ​​under different combinations of frequency (0.6GHz-1.5GHz) and voltage (0.75V-1.05V). Through polynomial fitting, a relationship model between power consumption and frequency and voltage is established. Taking the partial derivatives of this model, the sensitivity coefficient of power consumption to frequency is found to be 2.5W / GHz, and the sensitivity coefficient of power consumption to voltage is found to be 5.8W / V. Simultaneously, by analyzing power consumption data at different temperatures, the sensitivity coefficient of power consumption to temperature is calculated to be 0.02W / ℃.

[0149] When constructing the power regulation objective function, the system comprehensively considers the temperature gradient term and the power consumption deviation term. The temperature gradient term reflects the sum of squares of the temperature gradient in the hot spot region, with an initial value of 285℃² / mm²; the power consumption deviation term reflects the sum of squares of the differences between the actual power consumption and the target power consumption of each functional module, with an initial value of 2.4W². The system sets the weight of the temperature gradient to 0.7 and the weight of the power consumption deviation to 0.3, resulting in an initial objective function value of 202.5.

[0150] When using gradient descent for optimization, the system sets the learning rate to 0.05, the maximum number of iterations to 100, and the convergence threshold to 0.01. In each iteration, the system calculates the gradient of the objective function with respect to each operating parameter and updates the operating parameters in the opposite direction of the gradient. For example, after the first iteration, the frequency of the graphics processing module decreases from 1.2GHz to 1.15GHz, and the voltage decreases from 0.95V to 0.92V; the frequency of the memory controller decreases from 800MHz to 780MHz, and the voltage decreases from 0.85V to 0.83V. After 35 iterations, the objective function value decreases to 58.6, satisfying the convergence condition. The final target operating parameters for the graphics processing module are: frequency 1.05GHz, voltage 0.88V; and for the memory controller, frequency 720MHz, voltage 0.81V. Based on these target operating parameters, the system generates power adjustment commands and sends them to the control registers of each functional module via the system bus.

[0151] While monitoring the temperature distribution of the adjusted hotspot area in real time, the system continues to collect temperature data with a period of 100ms. Calculating the temperature changes before and after the adjustment, it was found that the highest temperature in the hotspot area decreased from 85℃ to 78℃, and the maximum temperature gradient decreased from 12℃ / mm to 7℃ / mm. The system calculates the temperature control error, i.e., the deviation between the actual temperature and the target temperature of 75℃, which is currently 3℃. Based on the PID control algorithm, the system sets the proportional coefficient Kp=0.1, the integral coefficient Ki=0.05, and the derivative coefficient Kd=0.02, calculating the integral term of the temperature control error as 15℃·s, the proportional term as 3℃, and the derivative term as -0.8℃ / s, generating a comprehensive control quantity of 0.846.

[0152] Based on this comprehensive control input, the system corrects the target operating parameters. For the graphics processing module, the frequency is further reduced to 1.02GHz and the voltage to 0.86V; for the memory controller, the frequency is reduced to 705MHz and the voltage to 0.8V. The system updates the power adjustment command and sends it. After five PID adjustment cycles, the highest temperature in the hotspot area stabilizes at around 76℃, the maximum temperature gradient stabilizes at around 5.5℃ / mm, and the performance loss of each functional module is controlled within 8%, achieving the goal of balancing temperature control and performance.

[0153] The electrothermal co-verification system for heterogeneous 3D stacked interconnect structures according to embodiments of the present invention includes:

[0154] The first unit is used to collect physical layout information of the heterogeneous process 3D stacked interconnect structure and generate initial stress distribution data;

[0155] The second unit is used to calculate the migration path of charge carriers in each chip cell based on the initial stress distribution data, obtain the charge carrier migration trajectory matrix, and calculate the actual conductivity distribution of each chip cell.

[0156] The third unit is used to divide the heterogeneous process 3D stacked interconnect structure into multiple grid units, and to perform grid subdivision operation in areas where the current density is greater than a preset density threshold to generate grid unit data.

[0157] The fourth unit is used to calculate the current density and Joule heat loss of each grid cell based on the grid cell data and the actual conductivity distribution, obtain the initial heat distribution data, and calculate the initial temperature value of each grid cell.

[0158] The fifth unit is used to set temperature monitoring points in the grid cells to collect temperature data in real time and compare it with the initial temperature value. Grid cells with a temperature change rate exceeding a preset temperature change threshold and their corresponding adjacent cells are marked as hotspot areas.

[0159] The sixth unit is used to calculate the temperature gradient of the hot spot area, and generate power adjustment instructions by combining the power consumption sensitivity parameters of the functional modules, so as to dynamically adjust the working parameters of the functional modules.

[0160] The seventh unit is used to recalculate the temperature distribution of each grid cell based on the adjusted operating parameters, repeat temperature monitoring and power adjustment until the temperature of all monitoring points is lower than the preset temperature threshold, and output the verification results.

[0161] A third aspect of the present invention provides an electronic device, comprising:

[0162] processor;

[0163] Memory used to store processor-executable instructions;

[0164] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0165] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0166] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0167] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. An electro-thermal synergistic verification method for heterogeneous process 3D stacked interconnect structures, characterized in that, The method comprises the following steps: Collecting physical layout information of a heterogeneous process 3D stacked interconnection structure to generate initial stress distribution data; According to the initial stress distribution data, the migration path of the carrier in each layer of chip unit is calculated, the carrier migration trajectory matrix is obtained, and the actual conductivity distribution of each layer of chip unit is calculated; The heterogeneous process 3D stacked interconnection structure is divided into a plurality of grid units, and a grid subdivision operation is performed in the area where the current density is greater than the preset density threshold, and grid unit data is generated; According to the grid unit data and the actual conductivity distribution, the current density and the joule heat loss of each grid unit are calculated, the initial heat distribution data is obtained, and the initial temperature value of each grid unit is calculated; Temperature monitoring points are set in the grid unit to collect temperature data in real time, and compared with the initial temperature value, the grid unit and the corresponding adjacent unit whose temperature change rate exceeds the preset temperature change threshold are marked as a hot spot area; The temperature gradient of the hot spot area is calculated, the power adjustment instruction is generated combined with the power consumption sensitivity parameter of the functional module, and the working parameter of the functional module is dynamically adjusted; According to the adjusted working parameter, the temperature distribution of each grid unit is recalculated, and the temperature monitoring and power adjustment are repeated until the temperature of all monitoring points is lower than the preset temperature threshold, and the verification result is output.

2. The method of claim 1, wherein, According to the initial stress distribution data, the migration path of the carrier in each layer of chip unit is calculated, the carrier migration trajectory matrix is obtained, and the actual conductivity distribution of each layer of chip unit is calculated, which comprises: According to the initial stress distribution data, the effective mass tensor and the deformation potential tensor of the carrier are calculated, the modified drift diffusion equation is constructed, and the stress-induced potential is calculated; according to the modified drift diffusion equation and the stress-induced potential, the carrier motion equation is established, and the space-time discretization parameters are generated; According to the space-time discretization parameters and the carrier motion equation, a scattering effect model is constructed; according to the scattering effect model, the migration path of the carrier in each layer of chip unit is determined, the carrier migration trajectory matrix is generated, and the local conductivity parameters of each layer of chip unit, including electron mobility and hole mobility, are calculated; The local conductivity parameters are interpolated by using a spatial interpolation algorithm, and anisotropy correction is performed, to obtain the actual conductivity distribution of each layer of chip unit; When the calculation error of the actual conductivity distribution is greater than the preset error threshold, the time step parameter and the grid size parameter are adjusted, the space-time discretization update parameters are generated, and the scattering effect model is reconstructed according to the space-time discretization update parameters and the carrier motion equation, until the calculation error of the actual conductivity distribution is less than the preset error threshold, and the actual conductivity distribution is output.

3. The method of claim 2, wherein, According to the space-time discretization parameters and the carrier motion equation, a scattering effect model is constructed, which comprises: According to the space-time discretization parameters and the carrier motion equation, the phonon scattering potential and the impurity scattering potential field are obtained; based on the phonon scattering potential, the phonon occupation number is calculated, and based on the phonon occupation number and the impurity scattering potential field, the phonon scattering rate and the impurity scattering rate are calculated, a scattering intensity evaluation matrix is constructed, and a scattering intensity threshold is set; The scattering intensity in the scattering intensity evaluation matrix is compared with the scattering intensity threshold value, a dominant scattering mechanism is determined, the mean free path and the characteristic scattering length are calculated according to the dominant scattering mechanism, and a grid scale parameter is generated; According to the grid scale parameter, an electron wave function is constructed, and a corresponding eigenstate is calculated; based on the eigenstate, a scattering path sequence is constructed, the scattering probability of each scattering path is calculated, the scattering probability of the scattering path is compared with a preset probability threshold value, an effective scattering path is screened, the scattering time of the effective scattering path is calculated, and a scattering time sequence is generated; The relative error of two adjacent scattering time sequences is calculated, when the relative error is greater than a preset convergence threshold value, the grid scale parameter is updated, and the construction step is returned to be executed; until the relative error is less than the preset convergence threshold value, a final scattering effect model is output.

4. The method of claim 1, wherein, According to the grid element data and the actual electrical conductivity distribution, the current density and the joule heat loss of each grid element are calculated to obtain initial heat distribution data, and the initial temperature value of each grid element is calculated, including: According to the grid element data and the actual electrical conductivity distribution, the potential gradient is calculated, the grid element electric field strength is generated, and the grid element current density is calculated based on the grid element electric field strength; the joule heat power density is calculated according to the grid element current density and the grid element electric field strength, and the grid element heat source distribution data is generated; Based on the grid element heat source distribution data, a heat conduction equation is constructed, a thermal conductivity matrix is generated according to the specific heat capacity and the density of the material, and the grid element heat flux distribution is calculated; According to the grid element heat flux distribution, the heat diffusion rate is calculated, and the grid size parameter and the temperature field iteration parameter are determined based on the heat diffusion rate; The Hilbert-Huang transform model is constructed for the grid element heat flux distribution to obtain a frequency-energy density spectrum, an adaptive filter bank is constructed, the first frequency band heat flux and the second frequency band heat flux are divided by singular value decomposition, and the temperature field is solved respectively to obtain the first temperature distribution and the second temperature distribution; Based on the first temperature distribution and the second temperature distribution, a temperature field reconstruction weight function is constructed, and the grid element initial temperature value is calculated; When the amplitude ratio of the first temperature distribution to the second temperature distribution exceeds a preset ratio threshold value, the grid size parameter is reduced by a preset reduction value, and the temperature field iteration parameter is updated, and the Fourier transform step is returned to be executed.

5. The method of claim 4, wherein, The Hilbert-Huang transform model is constructed for the grid element heat flux distribution to obtain a frequency-energy density spectrum, an adaptive filter bank is constructed, the first frequency band heat flux and the second frequency band heat flux are divided by singular value decomposition, and the temperature field is solved respectively to obtain the first temperature distribution and the second temperature distribution, including: The Hilbert-Huang transform model is constructed for the grid element heat flux distribution, the intrinsic modal function of the heat flux is extracted, the instantaneous frequency feature is calculated, and the frequency-energy density spectrum is generated; An adaptive filter bank is constructed using the frequency-energy density spectrum, and the center frequency of the adaptive filter bank is updated according to the energy density distribution; The adaptive filter set is applied to the heat flux distribution to obtain a plurality of intrinsic frequency band components, the plurality of intrinsic frequency band components are singular value decomposed to extract a frequency band feature vector; singular values of the frequency band feature vector are calculated, and the singular values are divided into a first frequency band heat flux and a second frequency band heat flux according to a preset division boundary value; a temperature field equation of the first frequency band heat flux and the second frequency band heat flux is respectively constructed by using a variational multi-scale method; a temperature gradient value of a temperature field grid cell is calculated, and when the temperature gradient value exceeds a preset gradient threshold value, the corresponding grid cell is divided into four sub-grid cells, and the grid cell is updated; temperature field equations of the updated grid cells are respectively solved to obtain a first temperature distribution and a second temperature distribution.

6. The method of claim 1, wherein, The temperature gradient of the hotspot region is calculated, and a power adjustment instruction is generated in combination with the power consumption sensitivity parameter of the functional module to dynamically adjust the working parameters of the functional module, including: temperature data of each marked grid cell in the hotspot region is extracted, a temperature gradient vector between adjacent grid cells is calculated, a modulus value of the temperature gradient vector is calculated, and a temperature gradient distribution of the hotspot region is calibrated; real-time power consumption data and working parameter data of each functional module are obtained, and a correlation model of power consumption and working parameters is established; based on the correlation model, the power consumption sensitivity parameter of each functional module is calculated, including the sensitivity coefficient of temperature to power consumption and the sensitivity coefficient of power consumption to working parameters; a power adjustment target function is constructed based on the temperature gradient vector of the hotspot region and the power consumption sensitivity parameter, including a temperature gradient term and a power consumption deviation term; the power adjustment target function is optimized and solved by using the gradient descent method to obtain the target working parameters of the functional module, and a power adjustment instruction is generated based on the target working parameters; The adjusted temperature distribution of the hotspot region is monitored in real time, the temperature control error is calculated, the integrated term, the proportional term and the differential term of the temperature control error are calculated, and the comprehensive control amount is generated; based on the comprehensive control amount, the target working parameters are corrected, and the power adjustment instruction is updated.

7. An electro-thermal synergistic verification system for a heterogeneous process 3D stacked interconnect structure for implementing the method of any of the preceding claims 1-6, characterized in that, It includes: A first unit for collecting physical layout information of a heterogeneous process 3D stacked interconnection structure and generating initial stress distribution data; A second unit for calculating the migration path of a carrier in each layer of chip unit according to the initial stress distribution data to obtain a carrier migration trajectory matrix and calculate the actual conductivity distribution of each layer of chip unit; A third unit for dividing the heterogeneous process 3D stacked interconnection structure into a plurality of grid cells and performing grid subdivision operation in the region where the current density is greater than the preset density threshold to generate grid cell data; A fourth unit for calculating the current density and joule heat loss of each grid cell according to the grid cell data and the actual conductivity distribution to obtain initial heat distribution data and calculate the initial temperature value of each grid cell; A fifth unit for setting temperature monitoring points in the grid cells to collect temperature data in real time and compare with the initial temperature value, marking the grid cells and the corresponding adjacent cells with a temperature change rate exceeding a preset temperature change threshold as a hotspot region. The sixth unit is configured to calculate the temperature gradient of the hotspot area, combine the power consumption sensitivity parameter of the functional module to generate a power adjustment instruction, and dynamically adjust the working parameter of the functional module. The seventh unit is configured to recalculate the temperature distribution of each grid unit according to the adjusted working parameter, repeat the temperature monitoring and power adjustment until the temperature of all monitoring points is lower than the preset temperature threshold, and output a verification result.

8. An electronic device, comprising: The computer program product comprises: a processor; a memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to execute the method of any one of claims 1 to 6.

9. A computer-readable storage medium having stored thereon computer program instructions, wherein, The computer program instructions, when executed by the processor, implement the method of any one of claims 1 to 6.

Images